Address andyferris' comments

Author: mschauer

src/SDiagonal.jl

@@ -4,13 +4,13 @@
 import Base: getindex,setindex!,==,-,+,*,/,\,transpose,ctranspose,convert, size, abs, real, imag, conj, eye, inv
 import Base.LinAlg: ishermitian, issymmetric, isposdef, factorize, diag, trace, det, logdet, expm, logm, sqrtm
 
-@generated function scalem{T, M, N}(a::SMatrix{M,N, T}, b::SVector{N, T})
+@generated function scalem{M, N}(a::StaticMatrix{M,N}, b::StaticVector{N})
     expr = vec([:(a[$j,$i]*b[$i]) for j=1:M, i=1:N])
-    :(SMatrix{M,N,T}($(expr...)))
+    :(let val1 = ($(expr[1])); similar_type(SMatrix{M,N},typeof(val1))(val1, $(expr[2:end]...)); end)
 end
-@generated function scalem{T, M, N}(a::SVector{M,T}, b::SMatrix{M, N, T})
+@generated function scalem{M, N}(a::StaticVector{M}, b::StaticMatrix{M, N})
     expr = vec([:(b[$j,$i]*a[$j]) for j=1:M, i=1:N])
-    :(SMatrix{M,N,T}($(expr...)))
+    :(let val1 = ($(expr[1])); similar_type(SMatrix{M,N},typeof(val1))(val1, $(expr[2:end]...)); end)
 end
 
 struct SDiagonal{N,T} <: StaticMatrix{N, N, T}
@@ -22,35 +22,21 @@ diagtype{N}(::Type{SDiagonal{N}}) = SVector{N}
 diagtype(::Type{SDiagonal}) = SVector
 
 # this is to deal with convert.jl
-@inline (::Type{SD})(a::AbstractVector) where {SD <: SDiagonal} = SDiagonal(diagtype(SD)(a)) 
-@inline (::Type{SD})(a::Tuple) where {SD <: SDiagonal} = SDiagonal(diagtype(SD)(a)) 
+@inline (::Type{SD})(a::AbstractVector) where {SD <: SDiagonal} = SDiagonal(convert(diagtype(SD), a))
+@inline (::Type{SD})(a::Tuple) where {SD <: SDiagonal} = SDiagonal(convert(diagtype(SD), a))
 @inline (::Type{SDiagonal}){N,T}(a::SVector{N,T}) = SDiagonal{N,T}(a) 
 
-@generated function SDiagonal{N,T}(a::SMatrix{N,N,T})
+@generated function SDiagonal{N,T}(a::StaticMatrix{N,N,T})
     expr = [:(a[$i,$i]) for i=1:N]
     :(SDiagonal{N,T}($(expr...)))
 end
 
-
 convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal{N,T}) = D
-convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal) = SDiagonal{N,T}(convert(SVector{N,T}, D.diag))
-
-size{N}(D::SDiagonal{N}) = (N,N)
-
-function size{N}(D::SDiagonal{N},d::Int)
-    if d<1
-        throw(ArgumentError("dimension must be ≥ 1, got $d"))
-    end
-    return d<=2 ? N : 1
-end
+convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal{N}) = SDiagonal{N,T}(convert(SVector{N,T}, D.diag))
 
 Base.@propagate_inbounds function getindex{N,T}(D::SDiagonal{N,T}, i::Int, j::Int)  
     @boundscheck checkbounds(D, i, j)
-    if i == j
-        @inbounds return D.diag[i]
-    else
-        zero(T)
-    end
+    @inbounds return ifelse(i == j, D.diag[i], zero(T))
 end
 
 # avoid linear indexing?
@@ -76,11 +62,12 @@ factorize(D::SDiagonal) = D
 *{T<:Number}(D::SDiagonal, x::T) = SDiagonal(D.diag * x)
 /{T<:Number}(D::SDiagonal, x::T) = SDiagonal(D.diag / x)
 *(Da::SDiagonal, Db::SDiagonal) = SDiagonal(Da.diag .* Db.diag)
-*(D::SDiagonal, V::SVector) = D.diag .* V
-*(V::SVector, D::SDiagonal) = D.diag .* V
-*(A::SMatrix, D::SDiagonal) = scalem(A,D.diag)
-*(D::SDiagonal, A::SMatrix) = scalem(D.diag,A)
-\(D::SDiagonal, b::SVector) = D.diag .\ b
+*(D::SDiagonal, V::AbstractVector) = D.diag .* V
+*(D::SDiagonal, V::StaticVector) = D.diag .* V
+*(A::StaticMatrix, D::SDiagonal) = scalem(A,D.diag)
+*(D::SDiagonal, A::StaticMatrix) = scalem(D.diag,A)
+\(D::SDiagonal, b::AbstractVector) = D.diag .\ b
+\(D::SDiagonal, b::StaticVector) = D.diag .\ b # catch ambiguity
 
 conj(D::SDiagonal) = SDiagonal(conj(D.diag))
 transpose(D::SDiagonal) = D
@@ -101,17 +88,22 @@ expm(D::SDiagonal) = SDiagonal(exp.(D.diag))
 logm(D::SDiagonal) = SDiagonal(log.(D.diag))
 sqrtm(D::SDiagonal) = SDiagonal(sqrt.(D.diag))
 
-\(D::SDiagonal, B::SMatrix) = scalem(1 ./ D.diag, B)
-/(B::SMatrix, D::SDiagonal) = scalem(1 ./ D.diag, B)
+\(D::SDiagonal, B::StaticMatrix) = scalem(1 ./ D.diag, B)
+/(B::StaticMatrix, D::SDiagonal) = scalem(1 ./ D.diag, B)
 \(Da::SDiagonal, Db::SDiagonal) = SDiagonal(Db.diag ./ Da.diag)
 /(Da::SDiagonal, Db::SDiagonal) = SDiagonal(Da.diag ./ Db.diag )
 
-function inv{N,T}(D::SDiagonal{N,T})
-    for i = 1:N
-        if D.diag[i] == zero(T)
-            throw(SingularException(i))
-        end
+
+@generated function check_singular{N,T}(D::SDiagonal{N,T})
+    expr = Expr(:block)
+    for i=1:N
+        push!(expr.args, :(@inbounds iszero(D.diag[$i]) && throw(Base.LinAlg.SingularException($i))))
     end
+    expr
+end
+
+function inv{N,T}(D::SDiagonal{N,T})
+    check_singular(D)
     SDiagonal(inv.(D.diag))
 end
 

test/SDiagonal.jl

@@ -16,6 +16,8 @@
 
     @testset "Methods" begin
 
+        @test StaticArrays.scalem(@SMatrix([1 1 1;1 1 1; 1 1 1]), @SVector [1,2,3]) === @SArray [1 2 3; 1 2 3; 1 2 3]
+    
         m = SDiagonal(@SVector [11, 12, 13, 14])
         m2 = diagm([11, 12, 13, 14])
 
@@ -54,6 +56,8 @@
         @test_throws Exception m[1] = 1
 
         @test m*b ==  @SVector [22,-12,26,14]
+        @test (b'*m)' ==  @SVector [22,-12,26,14]
+        
         @test m\b == m2\b
         @test m*m == m2*m
 
@@ -68,5 +72,7 @@
         @test m\m == eye(SDiagonal{4,Float64})
 
 
+        
+        
     end
 end