Unrolled cholesky factorization #243
Conversation
|
I'm also not sure if the |
Would using our own struct be suitable, in your opinion? (You can use
I think that might be helpful. Benchmarking over a range of matrix sizes is the way to answer this question. (There's some pretty unclear code in the |
|
Awesome PR, btw. I haven't looked at the code in detail yet, but thanks :) |
| _chol(Size(A), A.data) | ||
| end | ||
| @inline _chol(A::StaticMatrix{<:Any,<:Any,T}, ::Type{UpperTriangular}) where {T} = _chol(Size(A), A, UpperTriangular) | ||
| @inline _chol(A::StaticMatrix{<:Any,<:Any,T}, ::Type{LowerTriangular}) where {T} = _chol(Size(A), A, LowerTriangular) |
There was a problem hiding this comment.
Seems like these two methods could simply be:
@inline _chol(A::StaticMatrix, ::Type{UL}) where {UL} = _chol(Size(A), A, UL)| newtype = similar_type(A,T) | ||
| @generated function _chol(::Size{s}, A::StaticMatrix{<:Any,<:Any,T}, ::Type{UpperTriangular}) where {s, T} | ||
| if s[1] != s[2] | ||
| error("matrix must be square") |
There was a problem hiding this comment.
This check is not necessary; the ishermitian check does this before we enter this function.
| error("matrix must be square") | ||
| end | ||
| n = s[1] | ||
| TX = promote_type(typeof(chol(one(T), UpperTriangular)), Float32) |
There was a problem hiding this comment.
TX = typeof(chol(one(T)))| newtype = similar_type(A,T) | ||
| @generated function _chol(::Size{s}, A::StaticMatrix{<:Any, <:Any, T}, ::Type{LowerTriangular}) where {s, T} | ||
| if s[1] != s[2] | ||
| error("matrix must be square") |
| end | ||
|
|
||
| n = s[1] | ||
| TX = promote_type(typeof(chol(one(T), LowerTriangular)), Float32) |
| T = promote_type(typeof(sqrt(one(eltype(A)))), Float32) | ||
| newtype = similar_type(A,T) | ||
| ## Numbers | ||
| _chol(x::Number, uplo) = Base.LinAlg._chol!(x, uplo) |
There was a problem hiding this comment.
Accidentally left this in from Base?
Yeah, I think this is probably the best way to go. I'll get the done and make the changes suggested by @fredrikekre. Thanks for the feedback by the way. |
|
Where is this at? |
|
Sorry about the delay, I've been away for a bit. I'll hopefully finish this up by the end of the week. |
|
There's a few conflicts but it appears there's still a lot of potential in this PR if it was modernized. We should have some benchmarks to decide when to stop unrolling. |
Unify getindex and view to use the same approach. This helps keep component types better than the Base slicing fallback.
This PR extends the
cholto higher dimensions, adds acholfactmethod and a\based on the triangular solvers in #222. Also,cholwill now return anUpperTriangularfixing #194.The
Choleskystruct in Base changed in v0.7 so I've targeted this later version with this PR. I'm not sure how to get this working on v0.6 and v0.7 simultaneously.I've kept the old tests for the time being (just to show that everything is consistent) but these should be removed when everything is ready to be merged.