SciML / LinearSolve.jl

LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface.
https://docs.sciml.ai/LinearSolve/stable/
Other
248 stars 53 forks source link

`QRFactorization` causes change in the `cacheval` after solve #299

Closed utkarsh530 closed 1 year ago

utkarsh530 commented 1 year ago

MWE:

using LinearSolve

prob = LinearProblem(rand(100,100), rand(100))

linsolve = init(prob, QRFactorization())

@show typeof(linsolve.cacheval)
# 0×0 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}
sol = solve(prob, QRFactorization())

@show typeof(sol.cache.cacheval)
# LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}

Causes a problem with OrdinaryDiffEq.jl due to the type change.

Possible reason:

there are two different factorization dispatch here: https://github.com/SciML/LinearSolve.jl/blob/f9389edccd2ca156d3f18c05481035a88a0d16f0/src/factorization.jl#L92-L106 This causes problems because init calls init_cacheval, which gives a different cache, and since solve doesn't have separate dispatch for QR, it calls do_factorization and sets a different cache.

ChrisRackauckas commented 1 year ago

So is qr_instance giving the wrong type? Is there an MWE of that?

utkarsh530 commented 1 year ago

Yes:

using LinearSolve

prob = LinearProblem(rand(100,100), rand(100))

linsolve = init(prob, QRFactorization())

@show typeof(linsolve.cacheval)

https://github.com/JuliaArrays/ArrayInterface.jl/blob/2be78bdfd7565151c698eb05f8417a50d5de37bd/src/ArrayInterface.jl#L576

qr_instance returns this type for Matrices.

The problem happens as follows:

  1. solve(linearprob) lowers to solve(init(linearprob...)), which sets linsolve.cacheval to the type returned by qr_instance.
  2. When solve(linsolve) happens it dispatches to: https://github.com/SciML/LinearSolve.jl/blob/f9389edccd2ca156d3f18c05481035a88a0d16f0/src/factorization.jl#L9 Which causes do_factorization to set the cache by qr as mentioned in the first comment.

Possible workaround: Changeqr_instance to initialize LinearAlgebra.QRCompactWY instead of LinearAlgebra.QRCompactWYQ.

rayegun commented 1 year ago

QR just got a bit of a rework in LinearAlgebra so the change is likely related to that?

ChrisRackauckas commented 1 year ago

Yes, this is fixed upstream https://github.com/JuliaArrays/ArrayInterface.jl/pull/413