LinearSolve interface

[vpuri3]

TODO

• interoperability solve(prob::ODEProblem,Rodas5(linsolve=alg::SciMLLinearSolveAlgorithm))
• wrap and test a KSP libraries - Krylov.jl, IterativeSolvers.jl, KrylovKit.jl. handle in-place, and out-of-place algs separately? save KSP solver cache in LinearCache.cacheval
• add present solution, iteration count, other functionality to LinearCache

ref #3

[vpuri3]

@ChrisRackauckas how to handle this?

from test/runtests.jl

    prob = ODEProblemLibrary.prob_ode_linear
    sol  = solve(prob, Rodas5(linsolve=SVDFactorization()); saveat=0.1)

  MethodError: no method matching copy(::SciMLBase.UDerivativeWrapper{ODEFunction
{false, DiffEqProblemLibrary.ODEProblemLibrary.var"#1#2", UniformScaling{Bool}, D
iffEqProblemLibrary.ODEProblemLibrary.var"#3#4", Nothing, Nothing, Nothing, Nothi
ng, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBa
se.DEFAULT_OBSERVED), Nothing}, Float64, Float64})                               
  Closest candidates are:                                                        
    copy(::DataStructures.MutableLinkedList{T}) where T at /Users/vp/.julia/packa
ges/DataStructures/nBjdy/src/mutable_list.jl:110                                 
    copy(::BandedMatrices.BandedLU{T, S}) where {T, S} at /Users/vp/.julia/packag
es/BandedMatrices/2lafE/src/banded/BandedLU.jl:48                                
    copy(::FillArrays.Zeros) at /Users/vp/.julia/packages/FillArrays/VLeUk/src/Fi
llArrays.jl:280                                                                  
    ...                                                                          
  Stacktrace:                                                                        [1] init(::LinearProblem{DataType, false, SciMLBase.UDerivativeWrapper{ODEFun
ction{false, DiffEqProblemLibrary.ODEProblemLibrary.var"#1#2", UniformScaling{Boo
l}, DiffEqProblemLibrary.ODEProblemLibrary.var"#3#4", Nothing, Nothing, Nothing, 
Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(Sc
iMLBase.DEFAULT_OBSERVED), Nothing}, Float64, Float64}, Float64, SciMLBase.NullPa
rameters, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::SVDFa
ctorization{LinearAlgebra.DivideAndConquer}; alias_A::Bool, alias_b::Bool, kwargs
::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})                 
      @ LinearSolve ~/.julia/dev/LinearSolve/src/common.jl:62

do we just define additional methods for these wrappers? if so, where to define them?

[ChrisRackauckas]

The ODE solvers would need to be updated to use LinearSolve.jl before you can pass a LinearSolve.jl solver to it.

[vpuri3]

how to go about this? or is it a later step and i should write more wrappers first?

[ChrisRackauckas]

There's a lot of random changes here. Please make the PR towards a concise point, and split to multiple PRs if necessary.

[vpuri3]

yeah i made changes at a few different places. will break it into pieces.

[ChrisRackauckas]

You can always retrigger CI by closing and reopening the PR BTW.

[vpuri3]

the tests are failing because cholesky! is modifying the A matrix which is being used in tests as @test $A * $x \approx $b. fixing this

[ChrisRackauckas]

these arguments should be handling the aliasing

[vpuri3]

line 60-61 of the same file - common.jl

    A = alias_A ? A : deepcopy(A)
    b = alias_b ? b : deepcopy(b)

[vpuri3]

here's the issue!

julia> using LinearAlgebra; n=8; A=Matrix(I,n,n)/2; fa=cholesky!(A);

julia> fa.L * fa.U
8×8 Matrix{Float64}:
 0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.5  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.5  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5

julia> A
8×8 Matrix{Float64}:
 0.707107  0.0       0.0       0.0       …  0.0       0.0       0.0
 0.0       0.707107  0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.707107  0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.707107     0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0       …  0.707107  0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.707107  0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.707107

julia> A * ones(n) # should be 0.5*ones(n)
8-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865476
 0.7071067811865476
 0.7071067811865476
 0.7071067811865476
 0.7071067811865476
 0.7071067811865476
 0.7071067811865476

A is being modified to something that it is not correct. I've created an issue here: JuliaLang/julia#43212. Will comment out the cholesky! test for the meantime.

[ChrisRackauckas]

No, the entire purpose of the mutating cholesky decomposition is that... it's mutating. Are you mutating the copied one?

[vpuri3]

ok so when we do fa = cholesky!(A), the space in A is where the factors are stored. So we can't just do A * x and the tests are written incorrectly.

[ChrisRackauckas]

I am digging in and noticed that the reason is that the API was completely bungled. Making the cache and algs callable are not only not part of the API, but it's also a really bad idea and the reason for the failures. This is all really far off so I'm taking this PR over to put it back on the right track.

[vpuri3]

I was attempting to mimic the api in diffeqbase/linear_nonlinear where you could do DefaultLinSolve()(x, A, b), where DefaultLinSolve is an algorithm, and the calls to set_A, set_b, and set_u would just be within alg()(x,A,b).

[ChrisRackauckas]

The design is not supposed to mimic that. It's supposed to be matching the common interface design.

https://scimlbase.sciml.ai/dev/#Common-Interface-High-Level

The whole purpose is to move linear solvers away from a hacky side design to something that follows the same flow as nonlinear solvers, diffeq solvers, optimization solvers, etc. so that the whole universe of numerical problems follows the same exact interface.

[vpuri3]

i agree with the point on api, and apologies for messing with the API without discussing.

based on last commit, the way to do things is:

cache = set_A(cache, A)
cache = set_b(cache, b)
cache = set_u(cache, u)
solve(cache)

would it make sense to allow for this?

function solve(prob, cache)
    cache = set_A(cache, prob.A)
    cache = set_b(cache, prob.b)
    cache = set_u(cache, prob.u0)
    return solve(cache)
end

prob = LinearProblem(A,b;u0=u)
solve(prob,cache)

[ChrisRackauckas]

I don't understand the use case. What is it?

[codecov]

Codecov Report

Merging #21 (c273237) into main (eb61248) will decrease coverage by 11.61%.
The diff coverage is 71.82%.

@@             Coverage Diff             @@
##             main      #21       +/-   ##
===========================================
- Coverage   84.74%   73.13%   -11.62%     
===========================================
  Files           3        5        +2     
  Lines          59      201      +142     
===========================================
+ Hits           50      147       +97     
- Misses          9       54       +45     

Impacted Files	Coverage Δ
src/default.jl	0.00% <0.00%> (ø)
src/LinearSolve.jl	40.00% <40.00%> (ø)
src/common.jl	79.41% <75.00%> (+5.33%)	⬆️
src/wrappers.jl	84.25% <84.25%> (ø)
src/factorization.jl	93.10% <100.00%> (+1.79%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update eb61248...c273237. Read the comment docs.