solving the same problem twice is broken

[karajan9]

On Julia 1.9.0 with LinearSolve 1.42.0 solving the same problem twice doesn't give the same result:

A = [1.0 2.0; 3.0 4.0]
b = [1.0, 2.0]
prob = LinearProblem(A, b)
linsolve = init(prob)
sol1 = solve(linsolve)
# julia> sol1.u
# 2-element Vector{Float64}:
#  0.0
#  0.5
sol2 = solve(linsolve)
# julia> sol2.u
# 2-element Vector{Float64}:
#  -10.999999999999995
#    8.499999999999996

This is also true when set_b is used in between.

A = [1.0 2.0; 3.0 4.0]
b = [1.0, 2.0]
prob = LinearProblem(A, b)
linsolve = init(prob)
sol1 = solve(linsolve)
# julia> sol1.u
# 2-element Vector{Float64}:
#  0.0
#  0.5
linsolve = LinearSolve.set_b(linsolve, [2.0, 3.0])
sol2 = solve(linsolve)
# julia> sol2.u
# 2-element Vector{Float64}:
# -15.999999999999993
#  12.499999999999995
# should be [-1.0, 1.5]

That means that the examples in the README are currently broken.

[DanielVandH]

Could be related to the cache being still treated as fresh after solving:

julia> A = [1.0 2.0; 3.0 4.0]
       b = [1.0, 2.0]
       prob = LinearProblem(A, b)
       linsolve = init(prob)
       solve(linsolve)
       linsolve.isfresh
true

[karajan9]

Interestingly,

sol = solve(linsolve)
sol.isfresh == false

[karajan9]

I also noticed that

A2 = [1.0 3.0; 3.0 4.0]
linsolve = LinearSolve.set_A(linsolve, A2)
solve(linsolve)

julia> A2
2×2 Matrix{Float64}:
 3.0       4.0
 0.333333  1.66667

if you don't deepcopy it.

Also, basically all tests fail as soon as you replace A in

# top of `basictests.jl`
n = 8
A = Matrix(I, n, n)
b = ones(n)
A1 = rand(n, n);
b1 = rand(n);
x1 = zero(b);
A2 = A / 2;
b2 = rand(n);
x2 = zero(b);

prob1 = LinearProblem(A1, b1; u0 = x1)
prob2 = LinearProblem(A2, b2; u0 = x2)

with a random matrix, for example. It baffles me a bit that apparently a lot of tests were only done with an identity matrix.

[DanielVandH]

Interesting... Indeed those tests should definitely include some random testing.

That first point is because set_A just does @set! cache.A = A, in contrast to the initial use of init which has an alias_A keyword argument to allow deepcopy to be used by default - the set_A and set_b functions need alias_A = default_alias_A(cache.alg, cache.A, cache.b) and alias_b = ... keyword arguments, just like init has. Could be a good issue to make separate to this one.

[DanielVandH]

These issues could be the same thing/related to what we're seeing here: #199, #181 (though IterativeSolvers.jl is apparently problematic, so maybe not this one), #179,

[karajan9]

Oh dear. The earliest version this code runs on is 0.1.1 (Nov 2021), it has the same issue. I don't want to know how many wrong answers may have be silently produced during that time.

[ChrisRackauckas]

The linsolve cache is immutable, so note that the documented interface is:

using LinearSolve
A = [1.0 2.0; 3.0 4.0]
b = [1.0, 2.0]
prob = LinearProblem(A, b)
linsolve = init(prob, LUFactorization(), alias_A = false, alias_b = false)
sol1 = solve(linsolve)
sol2 = solve(sol1.cache)
sol2.cache.u

which is fine and is the one shown in https://docs.sciml.ai/LinearSolve/stable/tutorials/caching_interface/, and is why #199 is fine.

This issue is that with an immutable cache it cannot flip the isfresh. Without mutability we cannot prevent this either without the linear algebra routines also being immutable (which would then not optimize in the way that is wanted, ruining the purpose for the immutability). We may need to just get rid of the immutability since people are seeming to treat LinearSolve caches like the mutable ones.

[ChrisRackauckas]

The v2 will be fixed and tested against this. Thanks for the report.

[karajan9]

Thanks for getting this fixed!
May I suggest updating the README? While the docs were correct, I got the wrong version from the README. I tried to make a PR but two things came up:

1. If you set A with cache.A = A2 then A2 will be modified on solve!. This is covered in the tests by doing cache.A = deepcopy(A2) but not mentioned in the docs. Which should be the suggested version?

2. The line sol2 = solve!(linsolve, IterativeSolversJL_GMRES()) errors while it works fine with the solve(prob... version.

[karajan9]

I will also add that while the test have been fixed now, they have been using the wrong README-version the whole time, which wasn't caught because the test in basictests used (and are using) only identity matrices.
Would it make sense to change that as well, or expand the tests? I can open another issue if that is preferred.

[ChrisRackauckas]

This isn't released yet. The README will get an update with the v2 upgrade notes.

If you set A with cache.A = A2 then A2 will be modified on solve!. This is covered in the tests by doing cache.A = deepcopy(A2) but not mentioned in the docs. Which should be the suggested version?

It depends on what you're trying to do. Do you want to modify existing memory? Then don't copy. That's required for a lot of contexts.

The line sol2 = solve!(linsolve, IterativeSolversJL_GMRES()) errors while it works fine with the solve(prob... version.

Example?

Would it make sense to change that as well, or expand the tests? I can open another issue if that is preferred.

Yes, it would make sense to do that.

[karajan9]

This isn't released yet. The README will get an update with the v2 upgrade notes.

In case there is a misunderstanding, I was talking about the v1 README which led me to write linsolve = LinearSolve.set_b(linsolve, [2.0, 3.0]) in my starting post. I wanted to update it to the v2 syntax.

Example?

I hope I'm using the new version correctly here.

] add LinearSolve#main IterativeSolvers

using IterativeSolvers
using LinearSolve
A = [1.0 2.0; 3.0 4.0]
b = [1.0, 2.0]
prob = LinearProblem(A, b)
linsolve = init(prob)

sol = solve!(linsolve, IterativeSolversJL_GMRES())
ERROR: TypeError: in setfield!, expected LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, got a value of type IterativeSolvers.GMRESIterable{SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Vector{Float64}, Vector{Float64}, Vector{Float64}, IterativeSolvers.ArnoldiDecomp{Float64, Matrix{Float64}}, IterativeSolvers.Residual{Float64, Float64}, Float64, ModifiedGramSchmidt}
Stacktrace:
 [1] setproperty!(cache::LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, name::Symbol, x::IterativeSolvers.GMRESIterable{SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Vector{Float64}, Vector{Float64}, Vector{Float64}, IterativeSolvers.ArnoldiDecomp{Float64, Matrix{Float64}}, IterativeSolvers.Residual{Float64, Float64}, Float64, ModifiedGramSchmidt})
   @ LinearSolve ~/.julia/packages/LinearSolve/5rdIU/src/common.jl:87
 [2] solve!(cache::LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, alg::IterativeSolversJL{typeof(IterativeSolvers.gmres_iterable!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ LinearSolve ~/.julia/packages/LinearSolve/5rdIU/src/iterative_wrappers.jl:272
 [3] solve!(cache::LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, SciMLOperators.IdentityOperator, SciMLOperators.IdentityOperator, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, alg::IterativeSolversJL{typeof(IterativeSolvers.gmres_iterable!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
   @ LinearSolve ~/.julia/packages/LinearSolve/5rdIU/src/iterative_wrappers.jl:267
 [4] top-level scope
   @ REPL[11]:1

sol = solve(prob, IterativeSolversJL_GMRES())
u: 2-element Vector{Float64}:
 2.498001805406602e-16
 0.4999999999999999

Yes, it would make sense to do that.

Alright 👍