Skip to content

Commit 0cd325a

Browse files
committed
Make inv similar to det and clear noise
1 parent 6301885 commit 0cd325a

File tree

2 files changed

+30
-27
lines changed

2 files changed

+30
-27
lines changed

src/det.jl

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,23 @@ end
2121
return vecdot(x0, cross(x1, x2))
2222
end
2323

24+
@inline function _det(::Size{(4,4)}, A::StaticMatrix,S::Type)
25+
A = similar_type(A,S)(A)
26+
@inbounds return (
27+
A[13] * A[10] * A[7] * A[4] - A[9] * A[14] * A[7] * A[4] -
28+
A[13] * A[6] * A[11] * A[4] + A[5] * A[14] * A[11] * A[4] +
29+
A[9] * A[6] * A[15] * A[4] - A[5] * A[10] * A[15] * A[4] -
30+
A[13] * A[10] * A[3] * A[8] + A[9] * A[14] * A[3] * A[8] +
31+
A[13] * A[2] * A[11] * A[8] - A[1] * A[14] * A[11] * A[8] -
32+
A[9] * A[2] * A[15] * A[8] + A[1] * A[10] * A[15] * A[8] +
33+
A[13] * A[6] * A[3] * A[12] - A[5] * A[14] * A[3] * A[12] -
34+
A[13] * A[2] * A[7] * A[12] + A[1] * A[14] * A[7] * A[12] +
35+
A[5] * A[2] * A[15] * A[12] - A[1] * A[6] * A[15] * A[12] -
36+
A[9] * A[6] * A[3] * A[16] + A[5] * A[10] * A[3] * A[16] +
37+
A[9] * A[2] * A[7] * A[16] - A[1] * A[10] * A[7] * A[16] -
38+
A[5] * A[2] * A[11] * A[16] + A[1] * A[6] * A[11] * A[16])
39+
end
40+
2441
@inline function _det(::Size, A::StaticMatrix,::Type)
2542
return det(Matrix(A))
2643
end

src/inv.jl

Lines changed: 13 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,18 @@
1-
@inline inv(A::StaticMatrix) = _inv(Size(A), A)
2-
3-
@inline _inv(::Size{(1,1)}, A) = similar_type(typeof(A), typeof(inv(one(eltype(A)))))(inv(A[1]))
4-
5-
@inline function _inv(::Size{(2,2)}, A)
1+
@inline function inv(A::StaticMatrix)
62
T = eltype(A)
73
S = typeof((one(T)*zero(T) + zero(T))/one(T))
8-
newtype = similar_type(A, S)
4+
_inv(Size(A),A,S)
5+
end
6+
7+
@inline _inv(::Size{(1,1)}, A, S::Type) = inv(A[1])
98

9+
@inline function _inv(::Size{(2,2)}, A, S::Type)
10+
newtype = similar_type(A, S)
1011
d = det(A)
1112
@inbounds return newtype((A[4]/d, -(A[2]/d), -(A[3]/d), A[1]/d))
1213
end
1314

14-
@inline function _inv(::Size{(3,3)}, A)
15-
T = eltype(A)
16-
S = typeof((one(T)*zero(T) + zero(T))/one(T))
15+
@inline function _inv(::Size{(3,3)}, A,S::Type)
1716
newtype = similar_type(A, S)
1817

1918
@inbounds x0 = SVector{3,S}(A[1], A[2], A[3])
@@ -30,7 +29,8 @@ end
3029
@inbounds return newtype((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3]))
3130
end
3231

33-
@inline function _inv(::Size{(4,4)}, A)
32+
@inline function _inv(::Size{(4,4)}, A, S::Type)
33+
A = similar_type(A,S)(A)
3434
idet = 1/det(A)
3535
B = @SMatrix [
3636
(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
5252
(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
5353
(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
5454
(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]
55-
return SMatrix{4,4,eltype(B),16}(B)
55+
return A = similar_type(A,S)(B)
5656
end
5757

58-
@inline function _inv(::Size, A)
59-
T = eltype(A)
60-
S = typeof((one(T)*zero(T) + zero(T))/one(T))
61-
AA = convert(Array{S}, A) # lufact() doesn't work with StaticArrays at the moment... and the branches below must be type-stable
62-
if istriu(A)
63-
Ai_ut = inv(UpperTriangular(AA))
64-
# TODO double check these routines leave the parent in a clean (upper triangular) state
65-
return Size(A)(parent(Ai_ut)) # Return a `SizedArray`
66-
elseif istril(A)
67-
Ai_lt = inv(LowerTriangular(AA))
68-
# TODO double check these routines leave the parent in a clean (lower triangular) state
69-
return Size(A)(parent(Ai_lt)) # Return a `SizedArray`
70-
else
71-
Ai_lu = inv(lufact(AA))
72-
return Size(A)(Ai_lu) # Return a `SizedArray`
73-
end
58+
@inline function _inv(::Size, A ,::Type)
59+
inv(Matrix(A))
7460
end

0 commit comments

Comments
 (0)