Merge remote-tracking branch 'upstream/master'

Author: mateuszbaran

src/MArray.jl

@@ -73,16 +73,6 @@ end
 
 @inline MArray(a::StaticArray) = MArray{size_tuple(Size(a))}(Tuple(a))
 
-# Simplified show for the type
-#show(io::IO, ::Type{MArray{S, T, N}}) where {S, T, N} = print(io, "MArray{$S,$T,$N}")
-
-# Some more advanced constructor-like functions
-@inline one(::Type{MArray{S}}) where {S} = one(MArray{S,Float64,tuple_length(S)})
-@inline one(::Type{MArray{S,T}}) where {S,T} = one(MArray{S,T,tuple_length(S)})
-
-# MArray(I::UniformScaling) methods to replace eye
-(::Type{MA})(I::UniformScaling) where {MA<:MArray} = _eye(Size(MA), MA, I)
-
 ####################
 ## MArray methods ##
 ####################

src/MVector.jl

@@ -95,7 +95,7 @@ macro MVector(ex)
                 error("@MVector expected a 1-dimensional array expression")
             end
         else
-            error("@MVector only supports the zeros(), ones(), rand(), randn(), randexp(), and eye() functions.")
+            error("@MVector only supports the zeros(), ones(), rand(), randn(), and randexp() functions.")
         end
     else
         error("Use @MVector [a,b,c] or @MVector([a,b,c])")

src/SArray.jl

@@ -52,16 +52,6 @@ end
 
 @inline SArray(a::StaticArray) = SArray{size_tuple(Size(a))}(Tuple(a))
 
-# Simplified show for the type
-# show(io::IO, ::Type{SArray{S, T, N}}) where {S, T, N} = print(io, "SArray{$S,$T,$N}") # TODO reinstate
-
-# Some more advanced constructor-like functions
-@inline one(::Type{SArray{S}}) where {S} = one(SArray{S, Float64, tuple_length(S)})
-@inline one(::Type{SArray{S, T}}) where {S, T} = one(SArray{S, T, tuple_length(S)})
-
-# SArray(I::UniformScaling) methods to replace eye
-(::Type{SA})(I::UniformScaling) where {SA<:SArray} = _eye(Size(SA), SA, I)
-
 ####################
 ## SArray methods ##
 ####################

src/SDiagonal.jl

@@ -33,7 +33,7 @@ size(::Type{<:SDiagonal{N}}, d::Int) where {N} = d > 2 ? 1 : N
 # override to avoid copying
 diag(D::SDiagonal) = D.diag
 
-# SDiagonal(I::UniformScaling) methods to replace eye
+# SDiagonal(I::UniformScaling) methods
 (::Type{SDiagonal{N}})(I::UniformScaling) where {N} = SDiagonal{N}(ntuple(x->I.λ, Val(N)))
 (::Type{SDiagonal{N,T}})(I::UniformScaling) where {N,T} = SDiagonal{N,T}(ntuple(x->I.λ, Val(N)))
 

src/SHermitianCompact.jl

@@ -230,7 +230,7 @@ end
     end
 end
 
-@inline _eye(s::Size{S}, t::Type{SSC}) where {S, SSC <: SHermitianCompact} = _one(s, t)
+@inline _scalar_matrix(s::Size{S}, t::Type{SSC}) where {S, SSC <: SHermitianCompact} = _one(s, t)
 
 # _fill covers fill, zeros, and ones:
 @generated function _fill(val, ::Size{s}, ::Type{SSC}) where {s, SSC <: SHermitianCompact}

src/SVector.jl

@@ -108,9 +108,9 @@ macro SVector(ex)
                 error("@SVector expected a 1-dimensional array expression")
             end
         else
-            error("@SVector only supports the zeros(), ones(), rand(), randn(), randexp(), and eye() functions.")
+            error("@SVector only supports the zeros(), ones(), rand(), randn() and randexp() functions.")
         end
-    else # TODO Expr(:call, :zeros), Expr(:call, :ones), Expr(:call, :eye) ?
+    else
         error("Use @SVector [a,b,c], @SVector Type[a,b,c] or a comprehension like [f(i) for i = i_min:i_max]")
     end
 end

src/StaticArrays.jl

@@ -19,11 +19,7 @@ import LinearAlgebra: transpose, adjoint, dot, eigvals, eigen, lyap, tr,
                       kron, diag, norm, dot, diagm, lu, svd, svdvals,
                       factorize, ishermitian, issymmetric, isposdef, normalize,
                       normalize!, Eigen, det, logdet, cross, diff, qr, \
-
-# import eye for deprecation warnings
-@static if isdefined(LinearAlgebra, :eye)
-    import LinearAlgebra: eye
-end
+using LinearAlgebra: checksquare
 
 export SOneTo
 export StaticScalar, StaticArray, StaticVector, StaticMatrix
@@ -88,9 +84,15 @@ const StaticMatrixLike{T} = Union{
     StaticMatrix{<:Any, <:Any, T},
     Transpose{T, <:StaticVecOrMat{T}},
     Adjoint{T, <:StaticVecOrMat{T}},
-    Symmetric{T, <:StaticMatrix{T}},
-    Hermitian{T, <:StaticMatrix{T}},
-    Diagonal{T, <:StaticVector{<:Any, T}}
+    Symmetric{T, <:StaticMatrix{<:Any, <:Any, T}},
+    Hermitian{T, <:StaticMatrix{<:Any, <:Any, T}},
+    Diagonal{T, <:StaticVector{<:Any, T}},
+    # We specifically list *Triangular here rather than using
+    # AbstractTriangular to avoid ambiguities in size() etc.
+    UpperTriangular{T, <:StaticMatrix{<:Any, <:Any, T}},
+    LowerTriangular{T, <:StaticMatrix{<:Any, <:Any, T}},
+    UnitUpperTriangular{T, <:StaticMatrix{<:Any, <:Any, T}},
+    UnitLowerTriangular{T, <:StaticMatrix{<:Any, <:Any, T}}
 }
 const StaticVecOrMatLike{T} = Union{StaticVector{<:Any, T}, StaticMatrixLike{T}}
 const StaticArrayLike{T} = Union{StaticVecOrMatLike{T}, StaticArray{<:Tuple, T}}

src/abstractarray.jl

@@ -1,13 +1,12 @@
-length(a::SA) where {SA <: StaticArrayLike} = length(SA)
+length(a::StaticArrayLike) = prod(Size(a))
 length(a::Type{SA}) where {SA <: StaticArrayLike} = prod(Size(SA))
 
-@pure size(::Type{SA}) where {SA <: StaticArrayLike} = get(Size(SA))
+@pure size(::Type{SA}) where {SA <: StaticArrayLike} = Tuple(Size(SA))
 @inline function size(t::Type{<:StaticArrayLike}, d::Int)
     S = size(t)
     d > length(S) ? 1 : S[d]
 end
-@inline size(a::StaticArrayLike) = size(typeof(a))
-@inline size(a::StaticArrayLike, d::Int) = size(typeof(a), d)
+@inline size(a::StaticArrayLike) = Tuple(Size(a))
 
 Base.axes(s::StaticArray) = _axes(Size(s))
 @pure function _axes(::Size{sizes}) where {sizes}
@@ -93,6 +92,27 @@ mutable_similar_type(::Type{T}, s::Size{S}, ::Type{Val{D}}) where {T,S,D} = MArr
 
 sizedarray_similar_type(::Type{T},s::Size{S},::Type{Val{D}}) where {T,S,D} = SizedArray{Tuple{S...},T,D,length(s)}
 
+# Utility for computing the eltype of an array instance, type, or type
+# constructor.  For type constructors without a definite eltype, the default
+# value is returned.
+Base.@pure _eltype_or(a::AbstractArray, default) = eltype(a)
+Base.@pure _eltype_or(::Type{<:AbstractArray{T}}, default) where {T} = T
+Base.@pure _eltype_or(::Type{<:AbstractArray}, default) = default # eltype not available
+
+"""
+    _construct_similar(a, ::Size, elements::NTuple)
+
+Construct a static array of similar type to `a` with the given `elements`.
+
+When `a` is an instance or a concrete type the element type `eltype(a)` is
+used. However, when `a` is a `UnionAll` type such as `SMatrix{2,2}`, the
+promoted type of `elements` is used instead.
+"""
+@inline function _construct_similar(a, s::Size, elements::NTuple{L,ET}) where {L,ET}
+    similar_type(a, _eltype_or(a, ET), s)(elements)
+end
+
+
 # Field vectors are user controlled, and currently default to SVector, etc
 
 """

src/det.jl

@@ -43,7 +43,7 @@ end
 end
 
 @generated function _det(::Size{S}, A::StaticMatrix) where S
-    LinearAlgebra.checksquare(A)
+    checksquare(A)
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
@@ -58,7 +58,7 @@ end
 @inline logdet(A::StaticMatrix) = _logdet(Size(A), A)
 @inline _logdet(::Union{Size{(1,1)}, Size{(2,2)}, Size{(3,3)}, Size{(4,4)}}, A::StaticMatrix) = log(det(A))
 @generated function _logdet(::Size{S}, A::StaticMatrix) where S
-    LinearAlgebra.checksquare(A)
+    checksquare(A)
     if prod(S) ≤ 14*14
         quote
             @_inline_meta

src/eigen.jl

@@ -13,6 +13,19 @@ end
     end
 end
 
+function eigvals(A::StaticMatrixLike, B::StaticMatrixLike; kwargs...)
+    SA = Size(A)
+    if SA != Size(B)
+        throw(DimensionMismatch("Generalized eigenvalues can only be calculated for matrices of equal sizes: dimensions are $SA and $(Size(B))"))
+    end
+    checksquare(A)
+    return _eigvals(SA, A, B; kwargs...)
+end
+
+@inline function _eigvals(s::Size, A::StaticMatrixLike, B::StaticMatrixLike; kwargs...)
+    return SVector{s[1]}(eigvals(Array(A), Array(B); kwargs...))
+end
+
 @inline _eigvals(::Size{(1,1)}, a, permute, scale) = @inbounds return SVector(Tuple(a))
 @inline _eigvals(::Size{(1, 1)}, a::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} = @inbounds return SVector(real(parent(a).data[1]))
 
@@ -110,25 +123,33 @@ end
     return SVector{s[1], T}(vals)
 end
 
+# Utility to rewrap `Eigen` of normal `Array` into an Eigen containing `SArray`.
+@inline function _make_static(s::Size, E::Eigen{T,V}) where {T,V}
+    Eigen(similar_type(SVector, V, Size(s[1]))(E.values),
+          similar_type(SMatrix, T, s)(E.vectors))
+end
 
-@inline function _eig(s::Size, A::StaticMatrix, permute, scale)
-    # Only cover the hermitian branch, for now at least
-    # This also solves some type-stability issues that arise in Base
+@inline function _eig(s::Size, A::T, permute, scale) where {T <: StaticMatrix}
     if ishermitian(A)
-       return _eig(s, Hermitian(A), permute, scale)
+        return _eig(s, Hermitian(A), permute, scale)
     else
-       error("Only hermitian matrices are diagonalizable by *StaticArrays*. Non-Hermitian matrices should be converted to `Array` first.")
+        # For the non-hermitian branch fall back to LinearAlgebra eigen().
+        # Eigenvalues could be real or complex so a Union of concrete types is
+        # inferred. Having _make_static a separate function allows inference to
+        # preserve the union of concrete types:
+        #   Union{E{A,B},E{C,D}} -> Union{E{SA,SB},E{SC,SD}}
+        _make_static(s, eigen(Array(A); permute = permute, scale = scale))
     end
 end
 
 @inline function _eig(s::Size, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
     E = eigen(Hermitian(Array(parent(A))))
-    return (SVector{s[1], T}(E.values), SMatrix{s[1], s[2], eltype(A)}(E.vectors))
+    return Eigen(SVector{s[1], T}(E.values), SMatrix{s[1], s[2], eltype(A)}(E.vectors))
 end
 
 
 @inline function _eig(::Size{(1,1)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
-    @inbounds return (SVector{1,T}((real(A[1]),)), SMatrix{1,1,eltype(A)}(I))
+    @inbounds return Eigen(SVector{1,T}((real(A[1]),)), SMatrix{1,1,eltype(A)}(I))
 end
 
 @inline function _eig(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
@@ -157,7 +178,7 @@ end
             vecs = @SMatrix [ v11  v21 ;
                               v12  v22 ]
 
-            return (vals, vecs)
+            return Eigen(vals, vecs)
         end
     else # A.uplo == 'L'
         if !iszero(a[2]) # A is not diagonal
@@ -181,7 +202,7 @@ end
             vecs = @SMatrix [ v11  v21 ;
                               v12  v22 ]
 
-            return (vals,vecs)
+            return Eigen(vals,vecs)
         end
     end
 
@@ -197,7 +218,7 @@ end
         vecs = @SMatrix [convert(TA, 0) convert(TA, 1);
                          convert(TA, 1) convert(TA, 0)]
     end
-    return (vals,vecs)
+    return Eigen(vals,vecs)
 end
 
 # A small part of the code in the following method was inspired by works of David
@@ -230,19 +251,19 @@ end
 
         if a11 < a22
             if a22 < a33
-                return (SVector(a11, a22, a33), hcat(v1,v2,v3))
+                return Eigen(SVector((a11, a22, a33)), hcat(v1,v2,v3))
             elseif a33 < a11
-                return (SVector(a33, a11, a22), hcat(v3,v1,v2))
+                return Eigen(SVector((a33, a11, a22)), hcat(v3,v1,v2))
             else
-                return (SVector(a11, a33, a22), hcat(v1,v3,v2))
+                return Eigen(SVector((a11, a33, a22)), hcat(v1,v3,v2))
             end
         else #a22 < a11
             if a11 < a33
-                return (SVector(a22, a11, a33), hcat(v2,v1,v3))
+                return Eigen(SVector((a22, a11, a33)), hcat(v2,v1,v3))
             elseif a33 < a22
-                return (SVector(a33, a22, a11), hcat(v3,v2,v1))
+                return Eigen(SVector((a33, a22, a11)), hcat(v3,v2,v1))
             else
-                return (SVector(a22, a33, a11), hcat(v2,v3,v1))
+                return Eigen(SVector((a22, a33, a11)), hcat(v2,v3,v1))
             end
         end
     end
@@ -379,12 +400,11 @@ end
         (eigvec1, eigvec3) = (eigvec3, eigvec1)
     end
 
-    return (SVector(eig1, eig2, eig3), hcat(eigvec1, eigvec2, eigvec3))
+    return Eigen(SVector(eig1, eig2, eig3), hcat(eigvec1, eigvec2, eigvec3))
 end
 
 @inline function eigen(A::StaticMatrix; permute::Bool=true, scale::Bool=true)
-    vals, vecs = _eig(Size(A), A, permute, scale)
-    return Eigen(vals, vecs)
+    _eig(Size(A), A, permute, scale)
 end
 
 # to avoid method ambiguity with LinearAlgebra
@@ -394,8 +414,7 @@ end
 @inline eigen(A::Symmetric{<:Complex,<:StaticMatrix}; kwargs...) = _eigen(A; kwargs...)
 
 @inline function _eigen(A::LinearAlgebra.HermOrSym; permute::Bool=true, scale::Bool=true)
-    vals, vecs = _eig(Size(A), A, permute, scale)
-    return Eigen(vals, vecs)
+    _eig(Size(A), A, permute, scale)
 end
 
 # NOTE: The following Boost Software License applies to parts of the method: