|
1 | 1 | # LU decomposition
|
2 | 2 | function lu(A::StaticMatrix, pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{true})
|
3 |
| - L,U,p = _lu(A, pivot) |
| 3 | + L,U,p = _lu(Size(A), A, pivot) |
4 | 4 | (L,U,p)
|
5 | 5 | end
|
6 | 6 |
|
7 | 7 | # For the square version, return explicit lower and upper triangular matrices.
|
8 | 8 | # We would do this for the rectangular case too, but Base doesn't support that.
|
9 | 9 | function lu(A::StaticMatrix{N,N}, pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{true}) where {N}
|
10 |
| - L,U,p = _lu(A, pivot) |
| 10 | + L,U,p = _lu(Size(A), A, pivot) |
11 | 11 | (LowerTriangular(L), UpperTriangular(U), p)
|
12 | 12 | end
|
13 | 13 |
|
14 | 14 |
|
15 |
| -@inline function _lu(A::StaticMatrix, pivot) |
| 15 | +@inline function _lu(::Size{S}, A::StaticMatrix, pivot) where {S} |
16 | 16 | # For now, just call through to Base.
|
17 | 17 | # TODO: statically sized LU without allocations!
|
18 | 18 | f = lufact(Matrix(A), pivot)
|
19 | 19 | T = eltype(A)
|
| 20 | + # Trick to get the output eltype - can't rely on the result of f[:L] as |
| 21 | + # it's not type inferrable. |
20 | 22 | T2 = typeof((one(T)*zero(T) + zero(T))/one(T))
|
21 | 23 | L = similar_type(A, T2, Size(Size(A)[1], diagsize(A)))(f[:L])
|
22 | 24 | U = similar_type(A, T2, Size(diagsize(A), Size(A)[2]))(f[:U])
|
23 |
| - p = similar_type(A, Size(Size(A)[1]))(f[:p]) |
| 25 | + p = similar_type(A, Int, Size(Size(A)[1]))(f[:p]) |
24 | 26 | (L,U,p)
|
25 | 27 | end
|
26 | 28 |
|
|
0 commit comments