Tidy one() implementation and UniformScaling constructors (#690)

* Remove all reference to the deprecated name `eye` (replace internal
  _eye with _scaler_matrix)
* Remove all the more-or-less duplicate definitions of one() and use
  _scalar_matrix as the implementation.
* Fix StaticMatrixLike definition to properly include Hermitian and
  Symmetric wrappers of static matrices.
* Import LinearAlgebra.checksquare() and use it directly.
* Remove obsolete checksquare specializations & add tests
* Introduce _construct_similar to alleviate the difficulty of computing
  eltype purely in the type domain with non-concrete UnionAlls.

Author: c42f

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,8 +84,8 @@ const StaticMatrixLike{T} = Union{
     StaticMatrix{<:Any, <:Any, T},
     Transpose{T, <:StaticVecOrMat{T}},
     Adjoint{T, <:StaticVecOrMat{T}},
-    Symmetric{T, <:StaticMatrix{T}},
-    Hermitian{T, <:StaticMatrix{T}},
+    Symmetric{T, <:StaticMatrix{<:Any, <:Any, T}},
+    Hermitian{T, <:StaticMatrix{<:Any, <:Any, T}},
     Diagonal{T, <:StaticVector{<:Any, T}}
 }
 const StaticVecOrMatLike{T} = Union{StaticVector{<:Any, T}, StaticMatrixLike{T}}

src/abstractarray.jl

@@ -93,6 +93,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/linalg.jl

@@ -1,6 +1,7 @@
 import Base: +, -, *, /, \
 
-# TODO: more operators, like AbstractArray
+#--------------------------------------------------
+# Vector space algebra
 
 # Unary ops
 @inline -(a::StaticArray) = map(-, a)
@@ -30,10 +31,7 @@ import Base: +, -, *, /, \
 @inline -(a::UniformScaling, b::StaticMatrix) = _plus_uniform(Size(b), -b, a.λ)
 
 @generated function _plus_uniform(::Size{S}, a::StaticMatrix, λ) where {S}
-    if S[1] != S[2]
-        throw(DimensionMismatch("matrix is not square: dimensions are $S"))
-    end
-    n = S[1]
+    n = checksquare(a)
     exprs = [i == j ? :(a[$(LinearIndices(S)[i, j])] + λ) : :(a[$(LinearIndices(S)[i, j])]) for i = 1:n, j = 1:n]
     return quote
         $(Expr(:meta, :inline))
@@ -46,6 +44,8 @@ end
 @inline \(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ \ b
 @inline /(a::StaticMatrix, b::UniformScaling) = a / b.λ
 
+#--------------------------------------------------
+# Matrix algebra
 
 # Transpose, conjugate, etc
 @inline conj(a::StaticArray) = map(conj, a)
@@ -85,31 +85,30 @@ end
 @inline Base.zero(a::SA) where {SA <: StaticArray} = zeros(SA)
 @inline Base.zero(a::Type{SA}) where {SA <: StaticArray} = zeros(SA)
 
-@inline one(::SM) where {SM <: StaticMatrix} = _one(Size(SM), SM)
-@inline one(::Type{SM}) where {SM <: StaticMatrix} = _one(Size(SM), SM)
-@generated function _one(::Size{S}, ::Type{SM}) where {S, SM <: StaticArray}
-    if (length(S) != 2) || (S[1] != S[2])
-        error("multiplicative identity defined only for square matrices")
-    end
-    T = eltype(SM) # should be "hyperpure"
-    if T == Any
-        T = Float64
-    end
-    exprs = [i == j ? :(one($T)) : :(zero($T)) for i ∈ 1:S[1], j ∈ 1:S[2]]
-    return quote
-        $(Expr(:meta, :inline))
-        SM(tuple($(exprs...)))
+@inline one(m::StaticMatrixLike) = _one(Size(m), m)
+@inline one(::Type{SM}) where {SM<:StaticMatrixLike}= _one(Size(SM), SM)
+function _one(s::Size, m_or_SM)
+    if (length(s) != 2) || (s[1] != s[2])
+        throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
     end
+    _scalar_matrix(s, m_or_SM, one(_eltype_or(m_or_SM, Float64)))
 end
 
-# StaticMatrix(I::UniformScaling) methods to replace eye
-(::Type{SM})(I::UniformScaling) where {N,M,SM<:StaticMatrix{N,M}} = _eye(Size(SM), SM, I)
-
-@generated function _eye(::Size{S}, ::Type{SM}, I::UniformScaling{T}) where {S, SM <: StaticArray, T}
-    exprs = [i == j ? :(I.λ) : :(zero($T)) for i ∈ 1:S[1], j ∈ 1:S[2]]
+# StaticMatrix(I::UniformScaling)
+(::Type{SM})(I::UniformScaling) where {SM<:StaticMatrix} = _scalar_matrix(Size(SM), SM, I.λ)
+# The following oddity is needed if we want `SArray{Tuple{2,3}}(I)` to work
+# because we do not have `SArray{Tuple{2,3}} <: StaticMatrix`.
+(::Type{SM})(I::UniformScaling) where {SM<:(StaticArray{Tuple{N,M}} where {N,M})} =
+    _scalar_matrix(Size(SM), SM, I.λ)
+
+# Construct a matrix with the scalar λ on the diagonal and zeros off the
+# diagonal. The matrix can be non-square.
+@generated function _scalar_matrix(s::Size{S}, m_or_SM, λ) where {S}
+    elements = Symbol[i == j ? :λ : :λzero for i in 1:S[1], j in 1:S[2]]
     return quote
         $(Expr(:meta, :inline))
-        SM(tuple($(exprs...)))
+        λzero = zero(λ)
+        _construct_similar(m_or_SM, s, tuple($(elements...)))
     end
 end
 
@@ -145,6 +144,8 @@ end
     end
 end
 
+#--------------------------------------------------
+# Vector products
 @inline cross(a::StaticVector, b::StaticVector) = _cross(same_size(a, b), a, b)
 _cross(::Size{S}, a::StaticVector, b::StaticVector) where {S} = error("Cross product not defined for $(S[1])-vectors")
 @inline function _cross(::Size{(2,)}, a::StaticVector, b::StaticVector)
@@ -179,6 +180,8 @@ end
     return ret
 end
 
+#--------------------------------------------------
+# Norms
 @inline LinearAlgebra.norm_sqr(v::StaticVector) = mapreduce(abs2, +, v; init=zero(real(eltype(v))))
 
 @inline norm(a::StaticArray) = _norm(Size(a), a)
@@ -240,9 +243,7 @@ end
 
 @inline tr(a::StaticMatrix) = _tr(Size(a), a)
 @generated function _tr(::Size{S}, a::StaticMatrix) where {S}
-    if S[1] != S[2]
-        throw(DimensionMismatch("matrix is not square"))
-    end
+    checksquare(a)
 
     if S[1] == 0
         return :(zero(eltype(a)))
@@ -257,6 +258,10 @@ end
     end
 end
 
+
+#--------------------------------------------------
+# Outer products
+
 const _length_limit = Length(200)
 
 @inline kron(a::StaticMatrix, b::StaticMatrix) = _kron(_length_limit, Size(a), Size(b), a, b)
@@ -414,11 +419,9 @@ end
     end
 end
 
-# some micro-optimizations (TODO check these make sense for v0.6+)
-@inline LinearAlgebra.checksquare(::SM) where {SM<:StaticMatrix} = _checksquare(Size(SM))
-@inline LinearAlgebra.checksquare(::Type{SM}) where {SM<:StaticMatrix} = _checksquare(Size(SM))
 
-@pure _checksquare(::Size{S}) where {S} = (S[1] == S[2] || throw(DimensionMismatch("matrix is not square: dimensions are $S")); S[1])
+#--------------------------------------------------
+# Some shimming for special linear algebra matrix types
+@inline LinearAlgebra.Symmetric(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Symmetric{eltype(A),typeof(A)}(A, uplo))
+@inline LinearAlgebra.Hermitian(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Hermitian{eltype(A),typeof(A)}(A, uplo))
 
-@inline LinearAlgebra.Symmetric(A::StaticMatrix, uplo::Char='U') = (LinearAlgebra.checksquare(A);Symmetric{eltype(A),typeof(A)}(A, uplo))
-@inline LinearAlgebra.Hermitian(A::StaticMatrix, uplo::Char='U') = (LinearAlgebra.checksquare(A);Hermitian{eltype(A),typeof(A)}(A, uplo))

src/qr.jl

@@ -98,7 +98,7 @@ end
     Q = [Symbol("Q_$(i)_$(j)") for i = 1:m, j = 1:m]
     R = [Symbol("R_$(i)_$(j)") for i = 1:m, j = 1:n]
 
-    initQ = [:($(Q[i, j]) = $(i == j ? one : zero)(T)) for i = 1:m, j = 1:m]  # Q .= eye(A)
+    initQ = [:($(Q[i, j]) = $(i == j ? one : zero)(T)) for i = 1:m, j = 1:m]  # Q .= I
     initR = [:($(R[i, j]) = T(A[$i, $j])) for i = 1:m, j = 1:n]               # R .= A
 
     code = quote end

src/util.jl

@@ -108,3 +108,4 @@ Base.@propagate_inbounds function invperm(p::StaticVector)
      ip[p] = 1:length(p)
      similar_type(p)(ip)
 end
+

test/linalg.jl

@@ -1,5 +1,7 @@
 using StaticArrays, Test, LinearAlgebra
 
+using LinearAlgebra: checksquare
+
 @testset "Linear algebra" begin
 
     @testset "SArray as a (mathematical) vector space" begin
@@ -100,7 +102,7 @@ using StaticArrays, Test, LinearAlgebra
         @test @inferred(one(MMatrix{2,2}))::MMatrix == @MMatrix [1.0 0.0; 0.0 1.0]
         @test @inferred(one(MMatrix{2}))::MMatrix == @MMatrix [1.0 0.0; 0.0 1.0]
 
-        @test_throws ErrorException one(MMatrix{2,4})
+        @test_throws DimensionMismatch one(MMatrix{2,4})
     end
 
     @testset "cross()" begin
@@ -265,4 +267,12 @@ using StaticArrays, Test, LinearAlgebra
         @test @inferred(kron(A,transpose(b)))::SizedMatrix{10,210} == kron(P,transpose(q))
 
     end
+
+    @testset "checksquare" begin
+        m22 = SA[1 2; 3 4]
+        @test @inferred(checksquare(m22)) === 2
+        @test_inlined checksquare(m22)
+        m23 = SA[1 2 3; 4 5 6]
+        @test_inlined checksquare(m23) false
+    end
 end

test/testutil.jl

@@ -58,7 +58,7 @@ Note that LLVM IR can contain `call` instructions to intrinsics which don't
 make it into the native code, so this can be overly eager in declaring a
 a lack of complete inlining.
 """
-macro test_inlined(ex)
+macro test_inlined(ex, should_inline=true)
     ex_orig = ex
     ex = macroexpand(@__MODULE__, :(@code_llvm $ex))
     expr = quote
@@ -71,7 +71,10 @@ macro test_inlined(ex)
         # TODO: Figure out some better pattern matching; LLVM IR can contain
         # calls to intrinsics, so this will sometimes/often fail even when the
         # native code has no call instructions.
-        @test !occursin("call", code_str)
+        $(should_inline ?
+          :(@test !occursin("call", code_str)) :
+          :(@test occursin("call", code_str))
+        )
     end
     @assert expr.args[4].head == :macrocall
     expr.args[4].args[2] = __source__
@@ -84,9 +87,11 @@ should_not_be_inlined(x) = _should_not_be_inlined(x)
 
 @testset "@test_inlined" begin
     @test_inlined should_be_inlined(1)
+    @test_inlined should_not_be_inlined(1) false
     ts = @testset ErrorCounterTestSet "" begin
+        @test_inlined should_be_inlined(1) false
         @test_inlined should_not_be_inlined(1)
     end
-    @test ts.errorcount == 0 && ts.failcount == 1 && ts.passcount == 0
+    @test ts.errorcount == 0 && ts.failcount == 2 && ts.passcount == 0
 end