|
| 1 | +@inline det(A::StaticMatrix) = _det(Size(A), A) |
1 | 2 |
|
2 |
| -@generated function det{T}(A::StaticMatrix{T}) |
3 |
| - if size(A) == (1,1) |
4 |
| - return quote |
5 |
| - $(Expr(:meta, :inline)) |
6 |
| - @inbounds return A[1] |
7 |
| - end |
8 |
| - elseif size(A) == (2,2) |
9 |
| - return quote |
10 |
| - $(Expr(:meta, :inline)) |
11 |
| - @inbounds return A[1]*A[4] - A[3]*A[2] |
12 |
| - end |
13 |
| - elseif size(A) == (3,3) |
14 |
| - return quote |
15 |
| - $(Expr(:meta, :inline)) |
16 |
| - #@inbounds a = A[5]*A[9] - A[8]*A[6] |
17 |
| - #@inbounds b = A[8]*A[3] - A[2]*A[9] |
18 |
| - #@inbounds c = A[2]*A[6] - A[5]*A[3] |
19 |
| - #@inbounds return A[1]*a + A[4]*b + A[7]*c |
| 3 | +""" |
| 4 | + det(Size(m,m), mat) |
20 | 5 |
|
21 |
| - @inbounds x0 = SVector(A[1], A[2], A[3]) |
22 |
| - @inbounds x1 = SVector(A[4], A[5], A[6]) |
23 |
| - @inbounds x2 = SVector(A[7], A[8], A[9]) |
24 |
| - return vecdot(x0, cross(x1, x2)) |
25 |
| - end |
26 |
| - else |
27 |
| - S = typeof((one(T)*zero(T) + zero(T))/one(T)) |
28 |
| - return quote # Implementation from Base |
29 |
| - if istriu(A) || istril(A) |
30 |
| - return convert($S, det(UpperTriangular(A)))::$S # Is this a Julia bug that a convert is not type stable?? |
31 |
| - end |
32 |
| - AA = convert(Array{$S}, A) |
33 |
| - return det(lufact(AA)) |
| 6 | +Calculate the matrix determinate using an algorithm specialized on the size of |
| 7 | +the `m`×`m` matrix `mat`, which is much faster for small matrices. |
| 8 | +""" |
| 9 | +@inline _det(::Size{(1,1)}, A::AbstractMatrix) = @inbounds return A[1] |
| 10 | + |
| 11 | +@inline function _det(::Size{(2,2)}, A::AbstractMatrix) |
| 12 | + @inbounds return A[1]*A[4] - A[3]*A[2] |
| 13 | +end |
| 14 | + |
| 15 | +@inline function _det(::Size{(3,3)}, A::AbstractMatrix) |
| 16 | + @inbounds x0 = SVector(A[1], A[2], A[3]) |
| 17 | + @inbounds x1 = SVector(A[4], A[5], A[6]) |
| 18 | + @inbounds x2 = SVector(A[7], A[8], A[9]) |
| 19 | + return vecdot(x0, cross(x1, x2)) |
| 20 | +end |
| 21 | + |
| 22 | +@generated function _det{S,T}(::Size{S}, A::AbstractMatrix{T}) |
| 23 | + if S[1] != S[2] |
| 24 | + throw(DimensionMismatch("matrix is not square")) |
| 25 | + end |
| 26 | + T2 = typeof((one(T)*zero(T) + zero(T))/one(T)) |
| 27 | + return quote # Implementation from Base |
| 28 | + if istriu(A) || istril(A) |
| 29 | + return convert($T2, det(UpperTriangular(A)))::$T2 # Is this a Julia bug that a convert is not type stable?? |
34 | 30 | end
|
| 31 | + AA = convert(Array{$T2}, A) |
| 32 | + return det(lufact(AA)) |
35 | 33 | end
|
36 | 34 | end
|
0 commit comments