Merge branch 'master' into pull-request/a0d168f3

# Conflicts:
#	src/det.jl
#	test/det.jl

Author: TsurHerman

REQUIRE

@@ -1 +1 @@
-julia 0.6.0-pre
+julia 0.6.0

src/StaticArrays.jl

@@ -8,7 +8,7 @@ import Base: getindex, setindex!, size, similar, vec, show,
              length, convert, promote_op, promote_rule, map, map!, reduce, reducedim, mapreducedim,
              mapreduce, broadcast, broadcast!, conj, transpose, ctranspose,
              hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill,
-             fill!, det, logdet, inv, eig, eigvals, expm, logm, sqrtm, lyap, trace, diag, vecnorm, norm, dot, diagm, diag,
+             fill!, det, logdet, inv, eig, eigvals, expm, logm, sqrtm, lyap, trace, kron, diag, vecnorm, norm, dot, diagm, diag,
              lu, svd, svdvals, svdfact, factorize, ishermitian, issymmetric, isposdef,
              sum, diff, prod, count, any, all, minimum,
              maximum, extrema, mean, copy, rand, randn, randexp, rand!, randn!,

src/det.jl

@@ -1,64 +1,26 @@
-@inline det(A::StaticMatrix) = _det(Size(A), A)
-@inline logdet(A::StaticMatrix) = _logdet(Size(A), A)
+@inline function det(A::StaticMatrix)
+    T = eltype(A)
+    S = typeof((one(T)*zero(T) + zero(T))/one(T))
+    _det(Size(A),A,S)
+end
 
-@inline _det(::Size{(1,1)}, A::StaticMatrix) = @inbounds return A[1]
+@inline logdet(A::StaticMatrix) = log(det(A))
 
-@inline function _det(::Size{(2,2)}, A::StaticMatrix)
-    @inbounds return A[1]*A[4] - A[3]*A[2]
-end
+@inline _det(::Size{(1,1)}, A::StaticMatrix,S::Type) = @inbounds return convert(S,A[1])
 
-@inline function _det(::Size{(2,2)}, A::StaticMatrix{<:Any, <:Any, <:Unsigned})
-    @inbounds return Signed(A[1]*A[4]) - Signed(A[3]*A[2])
+@inline function _det(::Size{(2,2)}, A::StaticMatrix, S::Type)
+    A = similar_type(A,S)(A)
+    @inbounds return A[1]*A[4] - A[3]*A[2]
 end
 
-@inline function _det(::Size{(3,3)}, A::StaticMatrix)
+@inline function _det(::Size{(3,3)}, A::StaticMatrix, S::Type)
+    A = similar_type(A,S)(A)
     @inbounds x0 = SVector(A[1], A[2], A[3])
     @inbounds x1 = SVector(A[4], A[5], A[6])
     @inbounds x2 = SVector(A[7], A[8], A[9])
     return vecdot(x0, cross(x1, x2))
 end
 
-@inline function _det(::Size{(3,3)}, A::StaticMatrix{<:Any, <:Any, <:Unsigned})
-    @inbounds x0 = SVector(Signed(A[1]), Signed(A[2]), Signed(A[3]))
-    @inbounds x1 = SVector(Signed(A[4]), Signed(A[5]), Signed(A[6]))
-    @inbounds x2 = SVector(Signed(A[7]), Signed(A[8]), Signed(A[9]))
-    return vecdot(x0, cross(x1, x2))
-end
-
-@inline function _det(::Size{(4,4)}, A::StaticMatrix)
-    @inbounds return (
-        A[13] * A[10]  * A[7]  * A[4]  - A[9] * A[14] * A[7]  * A[4]   -
-        A[13] * A[6]   * A[11] * A[4]  + A[5] * A[14] * A[11] * A[4]   +
-        A[9]  * A[6]   * A[15] * A[4]  - A[5] * A[10] * A[15] * A[4]   -
-        A[13] * A[10]  * A[3]  * A[8]  + A[9] * A[14] * A[3]  * A[8]   +
-        A[13] * A[2]   * A[11] * A[8]  - A[1] * A[14] * A[11] * A[8]   -
-        A[9]  * A[2]   * A[15] * A[8]  + A[1] * A[10] * A[15] * A[8]   +
-        A[13] * A[6]   * A[3]  * A[12] - A[5] * A[14] * A[3]  * A[12]  -
-        A[13] * A[2]   * A[7]  * A[12] + A[1] * A[14] * A[7]  * A[12]  +
-        A[5]  * A[2]   * A[15] * A[12] - A[1] * A[6]  * A[15] * A[12]  -
-        A[9]  * A[6]   * A[3]  * A[16] + A[5] * A[10] * A[3]  * A[16]  +
-        A[9]  * A[2]   * A[7]  * A[16] - A[1] * A[10] * A[7]  * A[16]  -
-        A[5]  * A[2]   * A[11] * A[16] + A[1] * A[6]  * A[11] * A[16])
-end
-
-@inline _logdet(S::Union{Size{(1,1)},Size{(2,2)},Size{(3,3)}}, A::StaticMatrix) = log(_det(S, A))
-
-for (symb, f) in [(:_det, :det), (:_logdet, :logdet)]
-    eval(quote
-        @generated function $symb{S}(::Size{S}, A::StaticMatrix)
-            if S[1] != S[2]
-                throw(DimensionMismatch("matrix is not square"))
-            end
-            return quote # Implementation from Base
-                @_inline_meta
-                T = eltype(A)
-                T2 = typeof((one(T)*zero(T) + zero(T))/one(T))
-                if istriu(A) || istril(A)
-                    return convert(T2, $($f)(UpperTriangular(A))) # Is this a Julia bug that a convert is not type stable??
-                end
-                AA = convert(Array{T2}, A)
-                return $($f)(lufact(AA))
-            end
-        end
-    end)
+@inline function _det(::Size, A::StaticMatrix,::Type)
+    return det(Matrix(A))
 end

src/inv.jl

@@ -16,9 +16,9 @@ end
     S = typeof((one(T)*zero(T) + zero(T))/one(T))
     newtype = similar_type(A, S)
 
-    @inbounds x0 = SVector(A[1], A[2], A[3])
-    @inbounds x1 = SVector(A[4], A[5], A[6])
-    @inbounds x2 = SVector(A[7], A[8], A[9])
+    @inbounds x0 = SVector{3,S}(A[1], A[2], A[3])
+    @inbounds x1 = SVector{3,S}(A[4], A[5], A[6])
+    @inbounds x2 = SVector{3,S}(A[7], A[8], A[9])
 
     y0 = cross(x1,x2)
     d  = vecdot(x0, y0)

src/linalg.jl

@@ -172,11 +172,18 @@ end
     end
 end
 
-@inline diagm(v::StaticVector) = _diagm(Size(v), v)
-@generated function _diagm(::Size{S}, v::StaticVector) where {S}
-    Snew = (S[1], S[1])
+@inline diagm(v::StaticVector, k::Type{Val{D}}=Val{0}) where {D} = _diagm(Size(v), v, k)
+@generated function _diagm(::Size{S}, v::StaticVector, ::Type{Val{D}}) where {S,D}
+    S1 = S[1]
+    Snew1 = S1+abs(D)
+    Snew = (Snew1, Snew1)
+    Lnew = Snew1 * Snew1
     T = eltype(v)
-    exprs = [i == j ? :(v[$i]) : zero(T) for i = 1:S[1], j = 1:S[1]]
+    ind = diagind(Snew1, Snew1, D)
+    exprs = fill(:(zero($T)), Lnew)
+    for n = 1:S[1]
+        exprs[ind[n]] = :(v[$n])
+    end
     return quote
         $(Expr(:meta, :inline))
         @inbounds return similar_type($v, Size($Snew))(tuple($(exprs...)))
@@ -310,6 +317,21 @@ end
     end
 end
 
+const _length_limit = Length(200)
+
+@inline kron(a::StaticMatrix, b::StaticMatrix) = _kron(_length_limit, Size(a), Size(b), a, b)
+@generated function _kron(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
+    outsize = SA .* SB
+    if prod(outsize) > length_limit
+        return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
+    end
+    rows = [:(hcat($([:(a[$(sub2ind(SA,i,j))]*b) for j=1:SA[2]]...))) for i=1:SA[1]]
+    return quote
+        @_inline_meta
+        @inbounds return vcat($(rows...))
+    end
+end
+
 @inline Size(::Union{RowVector{T, SA}, Type{RowVector{T, SA}}}) where {T, SA <: StaticArray} = Size(1, Size(SA)[1])
 @inline Size(::Union{RowVector{T, CA}, Type{RowVector{T, CA}}} where CA <: ConjVector{<:Any, SA}) where {T, SA <: StaticArray} = Size(1, Size(SA)[1])
 @inline Size(::Union{Symmetric{T,SA}, Type{Symmetric{T,SA}}}) where {T,SA<:StaticArray} = Size(SA)
@@ -323,4 +345,3 @@ end
 
 @inline Base.LinAlg.Symmetric(A::StaticMatrix, uplo::Char='U') = (Base.LinAlg.checksquare(A);Symmetric{eltype(A),typeof(A)}(A, uplo))
 @inline Base.LinAlg.Hermitian(A::StaticMatrix, uplo::Char='U') = (Base.LinAlg.checksquare(A);Hermitian{eltype(A),typeof(A)}(A, uplo))
-

src/util.jl

@@ -59,3 +59,34 @@ end
 
     return nothing
 end
+
+
+# Trivial view used to drop static dimensions to override dispatch
+struct TrivialView{A,T,N} <: AbstractArray{T,N}
+    a::A
+end
+
+size(a::TrivialView) = size(a.a)
+getindex(a::TrivialView, inds...) = getindex(a.a, inds...)
+setindex!(a::TrivialView, inds...) = setindex!(a.a, inds...)
+Base.IndexStyle(::Type{<:TrivialView{A}}) where {A} = IndexStyle(A)
+
+TrivialView(a::AbstractArray{T,N}) where {T,N} = TrivialView{typeof(a),T,N}(a)
+
+
+# Remove the static dimensions from an array
+
+"""
+    drop_sdims(a)
+
+Return an `AbstractArray` with the same elements as `a`, but with static
+dimensions removed (ie, not a `StaticArray`).
+
+This is useful if you want to override dispatch to call the `Base` version of
+operations such as `kron` instead of the implementation in `StaticArrays`.
+Normally you shouldn't need to do this, but it can be more efficient for
+certain algorithms where the number of elements of the output is a lot larger
+than the input.
+"""
+@inline drop_sdims(a::StaticArray) = TrivialView(a)
+@inline drop_sdims(a) = a

test/SizedArray.jl

@@ -68,6 +68,8 @@
         @test convert(Vector, SizedArray{Tuple{4}, Int, 1}(collect(3:6))) == collect(3:6)
         @test Matrix(SMatrix{2,2}((1,2,3,4))) == [1 3; 2 4]
         @test convert(Matrix, SMatrix{2,2}((1,2,3,4))) == [1 3; 2 4]
+        # Conversion after reshaping
+        @test_broken Array(SizedMatrix{2,2}([1,2,3,4])) == [1 3; 2 4]
     end
 
     @testset "promotion" begin

test/det.jl

@@ -1,21 +1,24 @@
 @testset "Determinant" begin
-    @test det(@SMatrix [1]) === 1
-    @test logdet(@SMatrix [1]) === 0.0
-    @test det(@SMatrix [0 1; 1 0]) === -1
+    @test det(@SMatrix [1]) == 1
+    @test logdet(@SMatrix [1]) == 0.0
+    @test det(@SMatrix [0 1; 1 0]) == -1
     @test logdet(@SMatrix Complex{Float64}[0 1; 1 0]) == log(det(@SMatrix Complex{Float64}[0 1; 1 0]))
 
-    @test det(@SMatrix [0 1 0; 1 0 0; 0 0 1]) === -1
+    @test det(@SMatrix [0 1 0; 1 0 0; 0 0 1]) == -1
     m = randn(Float64, 4,4)
     @test det(SMatrix{4,4}(m)) ≈ det(m)
     #triu/tril
-    @test det(@SMatrix [1 2; 0 3]) === 3
+    @test det(@SMatrix [1 2; 0 3]) == 3
     @test det(@SMatrix [1 2 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) == 400.0
     @test logdet(@SMatrix [1 2 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) ≈ log(400.0)
     @test @inferred(det(ones(SMatrix{10,10,Complex{Float64}}))) == 0
 
-    # Unsigned specializations
-    @test det(@SMatrix [0x00 0x01; 0x01 0x00])::Int8 == -1
-    @test det(@SMatrix [0x00 0x01 0x00; 0x01 0x00 0x00; 0x00 0x00 0x01])::Int8 == -1
+    # Unsigned specializations , compare to Base
+    M = @SMatrix [1 2 3 4; 200 5 6 7; 0 0 8 9; 0 0 0 10]
+    for sz in (2,3,4), typ in (UInt8,UInt16,UInt32,UInt64)
+        Mtag = SMatrix{sz,sz,typ}(M[1:sz,1:sz])
+        @test det(Mtag) == det(Array(Mtag))
+    end
 
     @test_throws DimensionMismatch det(@SMatrix [0; 1])
 end

test/inv.jl

@@ -70,6 +70,14 @@ end
     @test norm(Matrix(sm*inv(sm) - eye(4))) < 12*norm(m*inv(m) - eye(4))
 end
 
+@testset "Matrix inverse 5x5" begin
+    m = randn(Float64, 5,5) + eye(5)
+    @test inv(SMatrix{5,5}(m))::StaticMatrix ≈ inv(m)
+    m = triu(randn(Float64, 5,5) + eye(5))
+    @test inv(SMatrix{5,5}(m))::StaticMatrix ≈ inv(m)
+    m = tril(randn(Float64, 5,5) + eye(5))
+    @test inv(SMatrix{5,5}(m))::StaticMatrix ≈ inv(m)
+end
 
 #-------------------------------------------------------------------------------
 # More comprehensive but qualitiative testing for inv() accuracy

test/linalg.jl

@@ -40,6 +40,8 @@ using StaticArrays, Base.Test
 
     @testset "diagm()" begin
         @test @inferred(diagm(SVector(1,2))) === @SMatrix [1 0; 0 2]
+        @test @inferred(diagm(SVector(1,2,3), Val{2}))::SMatrix == diagm([1,2,3], 2)
+        @test @inferred(diagm(SVector(1,2,3), Val{-2}))::SMatrix == diagm([1,2,3], -2)
     end
 
     @testset "diag()" begin
@@ -147,4 +149,19 @@ using StaticArrays, Base.Test
         @test vecnorm(SVector{0, Float64}(()), 1) === 0.
         @test trace(SMatrix{0,0,Float64}(())) === 0.
     end
+
+    @testset "kron" begin
+        @test @inferred(kron(@SMatrix([1 2; 3 4]), @SMatrix([0 1 0; 1 0 1]))) ==
+            SMatrix{4,6,Int}([0 1 0 0 2 0;
+                              1 0 1 2 0 2;
+                              0 3 0 0 4 0;
+                              3 0 3 4 0 4])
+        @test @inferred(kron(@SMatrix([1 2; 3 4]), @SMatrix([2.0]))) === @SMatrix [2.0 4.0; 6.0 8.0]
+
+        # Output should be heap allocated into a SizedArray when it gets large
+        # enough.
+        M1 = collect(1:20)
+        M2 = collect(20:-1:1).'
+        @test @inferred(kron(SMatrix{20,1}(M1),SMatrix{1,20}(M2)))::SizedMatrix{20,20} == kron(M1,M2)
+    end
 end