From a0d168f3cdd0aebcdf2f0b1283cf1d764d133ef0 Mon Sep 17 00:00:00 2001 From: Tsur Herman Date: Sat, 29 Jul 2017 16:39:26 +0300 Subject: [PATCH 1/8] replace inv and det for 4x4 matrices, --- src/det.jl | 19 +++++++++++++++++-- src/inv.jl | 47 ++++++++++++++++++++++------------------------- test/det.jl | 4 ++-- test/inv.jl | 4 ++-- 4 files changed, 43 insertions(+), 31 deletions(-) diff --git a/src/det.jl b/src/det.jl index eaee9ad9..46efa889 100644 --- a/src/det.jl +++ b/src/det.jl @@ -25,6 +25,22 @@ end 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)] @@ -44,6 +60,5 @@ for (symb, f) in [(:_det, :det), (:_logdet, :logdet)] return $($f)(lufact(AA)) end end - end) + end) end - diff --git a/src/inv.jl b/src/inv.jl index 42861ea4..9b225fcb 100644 --- a/src/inv.jl +++ b/src/inv.jl @@ -30,32 +30,29 @@ end @inbounds return newtype((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3])) end -Base.@pure function splitrange(r::SUnitRange) - mid = (first(r) + last(r)) รท 2 - (SUnitRange(first(r), mid), SUnitRange(mid+1, last(r))) -end +@inline function _inv(::Size{(4,4)}, A) + idet = 1/det(A) + B = @SMatrix [ + (A[2,3]*A[3,4]*A[4,2] - A[2,4]*A[3,3]*A[4,2] + A[2,4]*A[3,2]*A[4,3] - A[2,2]*A[3,4]*A[4,3] - A[2,3]*A[3,2]*A[4,4] + A[2,2]*A[3,3]*A[4,4]) * idet + (A[2,4]*A[3,3]*A[4,1] - A[2,3]*A[3,4]*A[4,1] - A[2,4]*A[3,1]*A[4,3] + A[2,1]*A[3,4]*A[4,3] + A[2,3]*A[3,1]*A[4,4] - A[2,1]*A[3,3]*A[4,4]) * idet + (A[2,2]*A[3,4]*A[4,1] - A[2,4]*A[3,2]*A[4,1] + A[2,4]*A[3,1]*A[4,2] - A[2,1]*A[3,4]*A[4,2] - A[2,2]*A[3,1]*A[4,4] + A[2,1]*A[3,2]*A[4,4]) * idet + (A[2,3]*A[3,2]*A[4,1] - A[2,2]*A[3,3]*A[4,1] - A[2,3]*A[3,1]*A[4,2] + A[2,1]*A[3,3]*A[4,2] + A[2,2]*A[3,1]*A[4,3] - A[2,1]*A[3,2]*A[4,3]) * idet + + (A[1,4]*A[3,3]*A[4,2] - A[1,3]*A[3,4]*A[4,2] - A[1,4]*A[3,2]*A[4,3] + A[1,2]*A[3,4]*A[4,3] + A[1,3]*A[3,2]*A[4,4] - A[1,2]*A[3,3]*A[4,4]) * idet + (A[1,3]*A[3,4]*A[4,1] - A[1,4]*A[3,3]*A[4,1] + A[1,4]*A[3,1]*A[4,3] - A[1,1]*A[3,4]*A[4,3] - A[1,3]*A[3,1]*A[4,4] + A[1,1]*A[3,3]*A[4,4]) * idet + (A[1,4]*A[3,2]*A[4,1] - A[1,2]*A[3,4]*A[4,1] - A[1,4]*A[3,1]*A[4,2] + A[1,1]*A[3,4]*A[4,2] + A[1,2]*A[3,1]*A[4,4] - A[1,1]*A[3,2]*A[4,4]) * idet + (A[1,2]*A[3,3]*A[4,1] - A[1,3]*A[3,2]*A[4,1] + A[1,3]*A[3,1]*A[4,2] - A[1,1]*A[3,3]*A[4,2] - A[1,2]*A[3,1]*A[4,3] + A[1,1]*A[3,2]*A[4,3]) * idet + + (A[1,3]*A[2,4]*A[4,2] - A[1,4]*A[2,3]*A[4,2] + A[1,4]*A[2,2]*A[4,3] - A[1,2]*A[2,4]*A[4,3] - A[1,3]*A[2,2]*A[4,4] + A[1,2]*A[2,3]*A[4,4]) * idet + (A[1,4]*A[2,3]*A[4,1] - A[1,3]*A[2,4]*A[4,1] - A[1,4]*A[2,1]*A[4,3] + A[1,1]*A[2,4]*A[4,3] + A[1,3]*A[2,1]*A[4,4] - A[1,1]*A[2,3]*A[4,4]) * idet + (A[1,2]*A[2,4]*A[4,1] - A[1,4]*A[2,2]*A[4,1] + A[1,4]*A[2,1]*A[4,2] - A[1,1]*A[2,4]*A[4,2] - A[1,2]*A[2,1]*A[4,4] + A[1,1]*A[2,2]*A[4,4]) * idet + (A[1,3]*A[2,2]*A[4,1] - A[1,2]*A[2,3]*A[4,1] - A[1,3]*A[2,1]*A[4,2] + A[1,1]*A[2,3]*A[4,2] + A[1,2]*A[2,1]*A[4,3] - A[1,1]*A[2,2]*A[4,3]) * idet -@noinline function _inv(::Size{(4,4)}, A) - # Partition matrix into four 2x2 blocks. For 4x4 matrices this seems to be - # more stable than directly using the adjugate expansion. - # See http://www.freevec.org/function/inverse_matrix_4x4_using_partitioning - # - # TODO: This decomposition works in higher dimensions, but numerical - # stability doesn't seem good for badly conditioned matrices. Can be - # fixed? - (i1,i2) = splitrange(SUnitRange(1,4)) - @inbounds P = A[i1,i1] - @inbounds Q = A[i1,i2] - @inbounds R = A[i2,i1] - @inbounds S = A[i2,i2] - invP = inv(P) - invP_Q = invP*Q - S2 = inv(S - R*invP_Q) - R2 = -S2*(R*invP) - Q2 = -invP_Q*S2 - P2 = invP - invP_Q*R2 - [[P2 Q2]; - [R2 S2]] + (A[1,4]*A[2,3]*A[3,2] - A[1,3]*A[2,4]*A[3,2] - A[1,4]*A[2,2]*A[3,3] + A[1,2]*A[2,4]*A[3,3] + A[1,3]*A[2,2]*A[3,4] - A[1,2]*A[2,3]*A[3,4]) * idet + (A[1,3]*A[2,4]*A[3,1] - A[1,4]*A[2,3]*A[3,1] + A[1,4]*A[2,1]*A[3,3] - A[1,1]*A[2,4]*A[3,3] - A[1,3]*A[2,1]*A[3,4] + A[1,1]*A[2,3]*A[3,4]) * idet + (A[1,4]*A[2,2]*A[3,1] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] + A[1,1]*A[2,4]*A[3,2] + A[1,2]*A[2,1]*A[3,4] - A[1,1]*A[2,2]*A[3,4]) * idet + (A[1,2]*A[2,3]*A[3,1] - A[1,3]*A[2,2]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] + A[1,1]*A[2,2]*A[3,3]) * idet] + return SMatrix{4,4,eltype(B),16}(B) end @inline function _inv(::Size, A) diff --git a/test/det.jl b/test/det.jl index 2485b233..f9ba57c6 100644 --- a/test/det.jl +++ b/test/det.jl @@ -3,13 +3,13 @@ @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 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 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) === 400.0 + @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 diff --git a/test/inv.jl b/test/inv.jl index 25f00f17..bcca0f0d 100644 --- a/test/inv.jl +++ b/test/inv.jl @@ -66,8 +66,8 @@ end 3.42945636322714076e-01 3.29329508494837497e-01 2.25635541033857107e-01 6.03912636987153917e-02 4.77828437366344727e-01 4.86406974015710591e-01 1.95415684569693188e-01 6.74775080892497797e-02] sm = SMatrix{4,4}(m) - @test_broken norm(Matrix(inv(sm)*sm - eye(4))) < 10*norm(inv(m)*m - eye(4)) - @test_broken norm(Matrix(sm*inv(sm) - eye(4))) < 10*norm(m*inv(m) - eye(4)) + @test norm(Matrix(inv(sm)*sm - eye(4))) < 10*norm(inv(m)*m - eye(4)) + @test norm(Matrix(sm*inv(sm) - eye(4))) < 10*norm(m*inv(m) - eye(4)) end From 02b4890bde9c03374251ae6e3d76294e93717dc2 Mon Sep 17 00:00:00 2001 From: Tsur Date: Sun, 30 Jul 2017 15:48:11 +0300 Subject: [PATCH 2/8] Relax test for poorly conditioned 4x4 matrices --- test/inv.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/inv.jl b/test/inv.jl index bcca0f0d..14a02fc0 100644 --- a/test/inv.jl +++ b/test/inv.jl @@ -55,8 +55,8 @@ end 1.63327582123091175e-01 1.70552365412100865e-01 4.55720450934200327e-01 1.23299968650232419e-01 ] sm = SMatrix{4,4}(m) - @test norm(Matrix(sm*inv(sm) - eye(4))) < 10*norm(m*inv(m) - eye(4)) - @test norm(Matrix(inv(sm)*sm - eye(4))) < 10*norm(inv(m)*m - eye(4)) + @test norm(Matrix(sm*inv(sm) - eye(4))) < 12*norm(m*inv(m) - eye(4)) + @test norm(Matrix(inv(sm)*sm - eye(4))) < 12*norm(inv(m)*m - eye(4)) # Poorly conditioned matrix, generated by fuzz testing with # almost_singular_matrix(N, 3, 1e-7) @@ -66,8 +66,8 @@ end 3.42945636322714076e-01 3.29329508494837497e-01 2.25635541033857107e-01 6.03912636987153917e-02 4.77828437366344727e-01 4.86406974015710591e-01 1.95415684569693188e-01 6.74775080892497797e-02] sm = SMatrix{4,4}(m) - @test norm(Matrix(inv(sm)*sm - eye(4))) < 10*norm(inv(m)*m - eye(4)) - @test norm(Matrix(sm*inv(sm) - eye(4))) < 10*norm(m*inv(m) - eye(4)) + @test norm(Matrix(inv(sm)*sm - eye(4))) < 12*norm(inv(m)*m - eye(4)) + @test norm(Matrix(sm*inv(sm) - eye(4))) < 12*norm(m*inv(m) - eye(4)) end From 0cd325ac197dde9c521dff613281b8ddafcf36cf Mon Sep 17 00:00:00 2001 From: Tsur Herman Date: Sat, 5 Aug 2017 09:55:47 +0300 Subject: [PATCH 3/8] Make inv similar to det and clear noise --- src/det.jl | 17 +++++++++++++++++ src/inv.jl | 40 +++++++++++++--------------------------- 2 files changed, 30 insertions(+), 27 deletions(-) diff --git a/src/det.jl b/src/det.jl index 9b0caa58..8bf3b68a 100644 --- a/src/det.jl +++ b/src/det.jl @@ -21,6 +21,23 @@ end return vecdot(x0, cross(x1, x2)) end +@inline function _det(::Size{(4,4)}, A::StaticMatrix,S::Type) + A = similar_type(A,S)(A) + @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 function _det(::Size, A::StaticMatrix,::Type) return det(Matrix(A)) end diff --git a/src/inv.jl b/src/inv.jl index a28f11b4..2f4ecf3d 100644 --- a/src/inv.jl +++ b/src/inv.jl @@ -1,19 +1,18 @@ -@inline inv(A::StaticMatrix) = _inv(Size(A), A) - -@inline _inv(::Size{(1,1)}, A) = similar_type(typeof(A), typeof(inv(one(eltype(A)))))(inv(A[1])) - -@inline function _inv(::Size{(2,2)}, A) +@inline function inv(A::StaticMatrix) T = eltype(A) S = typeof((one(T)*zero(T) + zero(T))/one(T)) - newtype = similar_type(A, S) + _inv(Size(A),A,S) +end + +@inline _inv(::Size{(1,1)}, A, S::Type) = inv(A[1]) +@inline function _inv(::Size{(2,2)}, A, S::Type) + newtype = similar_type(A, S) d = det(A) @inbounds return newtype((A[4]/d, -(A[2]/d), -(A[3]/d), A[1]/d)) end -@inline function _inv(::Size{(3,3)}, A) - T = eltype(A) - S = typeof((one(T)*zero(T) + zero(T))/one(T)) +@inline function _inv(::Size{(3,3)}, A,S::Type) newtype = similar_type(A, S) @inbounds x0 = SVector{3,S}(A[1], A[2], A[3]) @@ -30,7 +29,8 @@ end @inbounds return newtype((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3])) end -@inline function _inv(::Size{(4,4)}, A) +@inline function _inv(::Size{(4,4)}, A, S::Type) + A = similar_type(A,S)(A) idet = 1/det(A) B = @SMatrix [ (A[2,3]*A[3,4]*A[4,2] - A[2,4]*A[3,3]*A[4,2] + A[2,4]*A[3,2]*A[4,3] - A[2,2]*A[3,4]*A[4,3] - A[2,3]*A[3,2]*A[4,4] + A[2,2]*A[3,3]*A[4,4]) * idet @@ -52,23 +52,9 @@ end (A[1,3]*A[2,4]*A[3,1] - A[1,4]*A[2,3]*A[3,1] + A[1,4]*A[2,1]*A[3,3] - A[1,1]*A[2,4]*A[3,3] - A[1,3]*A[2,1]*A[3,4] + A[1,1]*A[2,3]*A[3,4]) * idet (A[1,4]*A[2,2]*A[3,1] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] + A[1,1]*A[2,4]*A[3,2] + A[1,2]*A[2,1]*A[3,4] - A[1,1]*A[2,2]*A[3,4]) * idet (A[1,2]*A[2,3]*A[3,1] - A[1,3]*A[2,2]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] + A[1,1]*A[2,2]*A[3,3]) * idet] - return SMatrix{4,4,eltype(B),16}(B) + return A = similar_type(A,S)(B) end -@inline function _inv(::Size, A) - T = eltype(A) - S = typeof((one(T)*zero(T) + zero(T))/one(T)) - AA = convert(Array{S}, A) # lufact() doesn't work with StaticArrays at the moment... and the branches below must be type-stable - if istriu(A) - Ai_ut = inv(UpperTriangular(AA)) - # TODO double check these routines leave the parent in a clean (upper triangular) state - return Size(A)(parent(Ai_ut)) # Return a `SizedArray` - elseif istril(A) - Ai_lt = inv(LowerTriangular(AA)) - # TODO double check these routines leave the parent in a clean (lower triangular) state - return Size(A)(parent(Ai_lt)) # Return a `SizedArray` - else - Ai_lu = inv(lufact(AA)) - return Size(A)(Ai_lu) # Return a `SizedArray` - end +@inline function _inv(::Size, A ,::Type) + inv(Matrix(A)) end From 3fbb8547b8e64604fbc005e13132b4a275ba2477 Mon Sep 17 00:00:00 2001 From: Tsur Herman Date: Sat, 5 Aug 2017 09:59:43 +0300 Subject: [PATCH 4/8] speed improvement for 2x2 inv --- src/inv.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/inv.jl b/src/inv.jl index 2f4ecf3d..b625c836 100644 --- a/src/inv.jl +++ b/src/inv.jl @@ -7,9 +7,10 @@ end @inline _inv(::Size{(1,1)}, A, S::Type) = inv(A[1]) @inline function _inv(::Size{(2,2)}, A, S::Type) - newtype = similar_type(A, S) - d = det(A) - @inbounds return newtype((A[4]/d, -(A[2]/d), -(A[3]/d), A[1]/d)) + newtype = similar_type(A,S) + A = newtype(A) + idet = 1/det(A) + @inbounds return newtype((A[4]*idet, -(A[2]*idet), -(A[3]*idet), A[1]*idet)) end @inline function _inv(::Size{(3,3)}, A,S::Type) From cc34e5e76b07717fb544b127e96dd0c8342cdae1 Mon Sep 17 00:00:00 2001 From: Tsur Herman Date: Sat, 5 Aug 2017 10:22:31 +0300 Subject: [PATCH 5/8] similar_type convert on all inv methods including >4 and 1 --- src/inv.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/inv.jl b/src/inv.jl index b625c836..18e8c11d 100644 --- a/src/inv.jl +++ b/src/inv.jl @@ -4,7 +4,7 @@ _inv(Size(A),A,S) end -@inline _inv(::Size{(1,1)}, A, S::Type) = inv(A[1]) +@inline _inv(::Size{(1,1)}, A, S::Type) = similar_type(A,S)(inv(A[1])) @inline function _inv(::Size{(2,2)}, A, S::Type) newtype = similar_type(A,S) @@ -56,6 +56,6 @@ end return A = similar_type(A,S)(B) end -@inline function _inv(::Size, A ,::Type) - inv(Matrix(A)) +@inline function _inv(::Size, A , S::Type) + similar_type(A,S)(inv(Matrix(A))) end From e578b9b929446ea4ed42f4d58e4fe74de2e6ad25 Mon Sep 17 00:00:00 2001 From: Tsur Herman Date: Mon, 7 Aug 2017 09:40:37 +0300 Subject: [PATCH 6/8] Avoid variable name changing type, moved promotion to dispatching function --- src/det.jl | 16 +++++++--------- src/inv.jl | 29 ++++++++++++++--------------- 2 files changed, 21 insertions(+), 24 deletions(-) diff --git a/src/det.jl b/src/det.jl index 8bf3b68a..c05bac3a 100644 --- a/src/det.jl +++ b/src/det.jl @@ -1,28 +1,26 @@ @inline function det(A::StaticMatrix) T = eltype(A) S = typeof((one(T)*zero(T) + zero(T))/one(T)) - _det(Size(A),A,S) + A_S = similar_type(A,S)(A) + _det(Size(A_S),A_S) end @inline logdet(A::StaticMatrix) = log(det(A)) -@inline _det(::Size{(1,1)}, A::StaticMatrix,S::Type) = @inbounds return convert(S,A[1]) +@inline _det(::Size{(1,1)}, A::StaticMatrix) = @inbounds return A[1] -@inline function _det(::Size{(2,2)}, A::StaticMatrix, S::Type) - A = similar_type(A,S)(A) +@inline function _det(::Size{(2,2)}, A::StaticMatrix) @inbounds return A[1]*A[4] - A[3]*A[2] end -@inline function _det(::Size{(3,3)}, A::StaticMatrix, S::Type) - A = similar_type(A,S)(A) +@inline function _det(::Size{(3,3)}, A::StaticMatrix) @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{(4,4)}, A::StaticMatrix,S::Type) - A = similar_type(A,S)(A) +@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] + @@ -38,6 +36,6 @@ end A[5] * A[2] * A[11] * A[16] + A[1] * A[6] * A[11] * A[16]) end -@inline function _det(::Size, A::StaticMatrix,::Type) +@inline function _det(::Size, A::StaticMatrix) return det(Matrix(A)) end diff --git a/src/inv.jl b/src/inv.jl index 18e8c11d..ca42ae87 100644 --- a/src/inv.jl +++ b/src/inv.jl @@ -1,24 +1,24 @@ @inline function inv(A::StaticMatrix) T = eltype(A) S = typeof((one(T)*zero(T) + zero(T))/one(T)) - _inv(Size(A),A,S) + A_S = similar_type(A,S)(A) + _inv(Size(A_S),A_S) end -@inline _inv(::Size{(1,1)}, A, S::Type) = similar_type(A,S)(inv(A[1])) +@inline _inv(::Size{(1,1)}, A) = similar_type(A)(inv(A[1])) -@inline function _inv(::Size{(2,2)}, A, S::Type) - newtype = similar_type(A,S) - A = newtype(A) +@inline function _inv(::Size{(2,2)}, A) + newtype = similar_type(A) idet = 1/det(A) @inbounds return newtype((A[4]*idet, -(A[2]*idet), -(A[3]*idet), A[1]*idet)) end -@inline function _inv(::Size{(3,3)}, A,S::Type) - newtype = similar_type(A, S) +@inline function _inv(::Size{(3,3)}, A) + newtype = similar_type(A) - @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]) + @inbounds x0 = SVector{3}(A[1], A[2], A[3]) + @inbounds x1 = SVector{3}(A[4], A[5], A[6]) + @inbounds x2 = SVector{3}(A[7], A[8], A[9]) y0 = cross(x1,x2) d = vecdot(x0, y0) @@ -30,8 +30,7 @@ end @inbounds return newtype((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3])) end -@inline function _inv(::Size{(4,4)}, A, S::Type) - A = similar_type(A,S)(A) +@inline function _inv(::Size{(4,4)}, A) idet = 1/det(A) B = @SMatrix [ (A[2,3]*A[3,4]*A[4,2] - A[2,4]*A[3,3]*A[4,2] + A[2,4]*A[3,2]*A[4,3] - A[2,2]*A[3,4]*A[4,3] - A[2,3]*A[3,2]*A[4,4] + A[2,2]*A[3,3]*A[4,4]) * idet @@ -53,9 +52,9 @@ end (A[1,3]*A[2,4]*A[3,1] - A[1,4]*A[2,3]*A[3,1] + A[1,4]*A[2,1]*A[3,3] - A[1,1]*A[2,4]*A[3,3] - A[1,3]*A[2,1]*A[3,4] + A[1,1]*A[2,3]*A[3,4]) * idet (A[1,4]*A[2,2]*A[3,1] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] + A[1,1]*A[2,4]*A[3,2] + A[1,2]*A[2,1]*A[3,4] - A[1,1]*A[2,2]*A[3,4]) * idet (A[1,2]*A[2,3]*A[3,1] - A[1,3]*A[2,2]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] + A[1,1]*A[2,2]*A[3,3]) * idet] - return A = similar_type(A,S)(B) + return similar_type(A)(B) end -@inline function _inv(::Size, A , S::Type) - similar_type(A,S)(inv(Matrix(A))) +@inline function _inv(::Size, A) + similar_type(A)(inv(Matrix(A))) end From 045c1cfc2e6d347a63fc91f02de4b077e23cf8ba Mon Sep 17 00:00:00 2001 From: Tsur Herman Date: Tue, 8 Aug 2017 18:25:59 +0300 Subject: [PATCH 7/8] explicitly call convert instead of constructor --- src/det.jl | 2 +- src/inv.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/det.jl b/src/det.jl index c05bac3a..ccadb81d 100644 --- a/src/det.jl +++ b/src/det.jl @@ -1,7 +1,7 @@ @inline function det(A::StaticMatrix) T = eltype(A) S = typeof((one(T)*zero(T) + zero(T))/one(T)) - A_S = similar_type(A,S)(A) + A_S = convert(similar_type(A,S),A) _det(Size(A_S),A_S) end diff --git a/src/inv.jl b/src/inv.jl index ca42ae87..b74ad02a 100644 --- a/src/inv.jl +++ b/src/inv.jl @@ -1,7 +1,7 @@ @inline function inv(A::StaticMatrix) T = eltype(A) S = typeof((one(T)*zero(T) + zero(T))/one(T)) - A_S = similar_type(A,S)(A) + A_S = convert(similar_type(A,S),A) _inv(Size(A_S),A_S) end From 801b43cf1155500dc8a6c745e238e8d3c22d828e Mon Sep 17 00:00:00 2001 From: Chris Foster Date: Sat, 19 Aug 2017 06:03:57 +1000 Subject: [PATCH 8/8] Add a 4x4 permutation matrix as an inv() test --- test/inv.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/inv.jl b/test/inv.jl index 59b552d8..104a3f36 100644 --- a/test/inv.jl +++ b/test/inv.jl @@ -47,6 +47,14 @@ end sm = SMatrix{4,4}(m) @test isapprox(inv(sm)::StaticMatrix, inv(m), rtol=2e-15) + # A permutation matrix which can cause problems for inversion methods + # without pivoting, eg 2x2 block decomposition (see #250) + mperm = [1 0 0 0; + 0 0 1 0; + 0 1 0 0; + 0 0 0 1] + @test inv(SMatrix{4,4}(mperm))::StaticMatrix โ‰ˆ inv(mperm) rtol=2e-16 + # Poorly conditioned matrix; almost_singular_matrix(4, 3, 1e-7) m = [ 2.83056817904263402e-01 1.26822318692296848e-01 2.59665505365002547e-01 1.24524798964590747e-01