Implement Butterfly Factorization method using RecursiveFactorization.jl #785
Conversation
|
@oscardssmith does this look good? |
| B, U, V = cache.cacheval[2], cache.cacheval[3], cache.cacheval[4] | ||
| if cache.isfresh | ||
| @assert M==N "A must be square" | ||
| U, V, F, out = RecursiveFactorization.🦋workspace(A, b, B, U, V, alg.thread) |
There was a problem hiding this comment.
Why is this not just a struct and ! operation? It would be much easier to read. I assume this is all just in-place and non-allocating.
There was a problem hiding this comment.
as in make U, V, F, out a struct?
| end | ||
| end | ||
| A, B, U, V, F = cache.cacheval | ||
| sol = V * (F \ (U * b)) |
There was a problem hiding this comment.
This is allocating and not using TriangularSolve.jl?
There was a problem hiding this comment.
I think we found that TriangularSolve.jl was slower than this method, so we left it as is
|
this pr should be ready now, let me know if there are any other edits needed |
|
Core LTS tests fail. |
|
not sure the test fails are related to this, can't quite figure out what's wrong |
alg = LinearSolve.ButterflyFactorization
A = [10.0 14.0; 14.0 20.0]
Re-solve: Error During Test at /home/chrisrackauckas/github-runners/deepsea4-21/_work/LinearSolve.jl/LinearSolve.jl/test/resolve.jl:67
Test threw exception
Expression: (solve!(linsolve)).u ≈ expected
UndefVarError: `ws` not defined
Stacktrace:
[1] #solve!#3
@ ~/github-runners/deepsea4-21/_work/LinearSolve.jl/LinearSolve.jl/ext/LinearSolveRecursiveFactorizationExt.jl:119 [inlined]
[2] solve!
@ ~/github-runners/deepsea4-21/_work/LinearSolve.jl/LinearSolve.jl/ext/LinearSolveRecursiveFactorizationExt.jl:107 [inlined]
[3] #solve!#12
@ ~/github-runners/deepsea4-21/_work/LinearSolve.jl/LinearSolve.jl/src/common.jl:424 [inlined]
[4] solve!(::LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.ButterflyFactorization{true}, RecursiveFactorization.🦋workspace{Float64}, SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Float64, Bool, LinearSolve.LinearSolveAdjoint{Missing}})
@ LinearSolve ~/github-runners/deepsea4-21/_work/LinearSolve.jl/LinearSolve.jl/src/common.jl:423
[5] macro expansion
@ ~/github-runners/deepsea4-21/_work/_tool/julia/1.10.10/x64/share/julia/stdlib/v1.10/Test/src/Test.jl:669 [inlined]
[6] top-level scope
@ ~/github-runners/deepsea4-21/_work/LinearSolve.jl/LinearSolve.jl/test/resolve.jl:67that is definitely related to this |
|
yeah you're right, sorry. hopefully that's fixed now |
| M, N = size(A) | ||
| if cache.isfresh | ||
| @assert M==N "A must be square" | ||
| ws = RecursiveFactorization.🦋workspace(A, b) |
There was a problem hiding this comment.
I think we can reuse the ws even when cache is fresh. As long as the size is the same, I think we can just update A, and b and everything works.
Re-solve: Test Failed at /home/chrisrackauckas/github-runners/demeter4-21/_work/LinearSolve.jl/LinearSolve.jl/test/resolve.jl:67
Expression: (solve!(linsolve)).u ≈ expected
Evaluated: [-15.249999999999984, 10.749999999999988] ≈ [-2.0, 1.5]
Stacktrace:
[1] macro expansion
@ ~/github-runners/demeter4-21/_work/_tool/julia/1.10.10/x64/share/julia/stdlib/v1.10/Test/src/Test.jl:672 [inlined]
[2] top-level scope
@ ~/github-runners/demeter4-21/_work/LinearSolve.jl/LinearSolve.jl/test/resolve.jl:67
[3] include(mod::Module, _path::String)
@ Base ./Base.jl:495
[4] include(x::String)
@ Main.var"##Re-solve#227" ~/.julia/packages/SafeTestsets/raUNr/src/SafeTestsets.jl:28
[5] macro expansion
@ ~/.julia/packages/SafeTestsets/raUNr/src/SafeTestsets.jl:24 [inlined]
[6] macro expansion
@ ~/github-runners/demeter4-21/_work/_tool/julia/1.10.10/x64/share/julia/stdlib/v1.10/Test/src/Test.jl:1577 [inlined]
[7] top-level scope
@ ~/.julia/packages/SafeTestsets/raUNr/src/SafeTestsets.jl:24It looks like the re-solve failed. |
| (;A, b, ws, U, V, out, tmp, n) = workspace | ||
| mul!(tmp, U', b) | ||
| TriangularSolve.ldiv!(F, tmp, thread) | ||
| mul!(b, V, tmp) |
There was a problem hiding this comment.
Mutating b will make it incorrect for the second solve. It seems like you need another temporary vector to do this right?
Co-authored-by: Oscar Smith <[email protected]>
|
wait |
|
then use a view |
3 participants
Checklist
contributor guidelines, in particular the SciML Style Guide and
COLPRAC.
Additional context
Implements a new
ButterflyFactorization()algorithm