LinearSolve is changing their solutions

Specific Domains Numerics
1

I am solving some Linear Systems using LinearSolve.jl and its caching system. I was fine doing fine with this example

using LinearSolve
n = 4
A = rand(n, n)
b = [rand(n) for i ∈ 1:3]

## everything is 'true'
let
    linsolve_1 = LinearProblem(A, b[1])
    sol_1 = solve(linsolve_1)
    x1 = A\b[1]
    @info all(sol_1.u .≈ x1)


    linsolve_2 = LinearSolve.set_b(sol_1.cache, b[2])
    sol_2 = solve(linsolve_2)
    x2 = A\b[2]
    @info all(sol_2.u .≈ x2)


    linsolve_3 = LinearSolve.set_b(sol_1.cache, b[3])
    sol_3 = solve(linsolve_3)
    x3 = A\b[3]
    @info all(sol_3.u .≈ x3)
end

# returns
#[ Info: true
#[ Info: true
#[ Info: true

Until, I just crop my verification of each solution, and paste them in the end of the code, Suddendly, I have some false.

let
    linsolve_1 = LinearProblem(A, b[1])
    sol_1 = solve(linsolve_1)
    x1 = A\b[1]
    println(sol_1.u)

    linsolve_2 = LinearSolve.set_b(sol_1.cache, b[2])
    sol_2 = solve(linsolve_2)
    x2 = A\b[2]
    println(sol_2.u)

    linsolve_3 = LinearSolve.set_b(sol_1.cache, b[3])
    sol_3 = solve(linsolve_3)
    x3 = A\b[3]
    println(sol_3.u)

    println("---------")
    println(sol_1.u)
    println(sol_2.u)
    println(sol_3.u)

    @info all(sol_1.u .≈ x1)
    @info all(sol_2.u .≈ x2)
    @info all(sol_3.u .≈ x3)
end

# returns
#[-1.0672073084102143, -2.5268971546641725, 3.084186997230096, 1.1307573345753459]
#[-1.0438221505570733, -0.08229940246406439, 1.891254274527264, -0.13214959434703005]
#[-1.6954927011416538, -4.2716079880714295, 4.98481354047552, 0.9501886647386011]
#---------
#[-1.6954927011416538, -4.2716079880714295, 4.98481354047552, 0.9501886647386011]
#[-1.6954927011416538, -4.2716079880714295, 4.98481354047552, 0.9501886647386011]
#[-1.6954927011416538, -4.2716079880714295, 4.98481354047552, 0.9501886647386011]
#[ Info: false
#[ Info: false
#[ Info: true

As we can see by the println, something had changed in sol_1 and sol_2 to have the value of sol_3.

Is this the expected behavior? Because I would consider a bug, or at least, some ! symbol should appear somewhere.

Note: I tested with Julia v1.8.5 and 1.9.0-beta4

2

You’re reusing the same cache. That’s the purpose of the interface for reusing the cache! Don’t reuse the cache if you don’t want the cache reused :sweat_smile:.

3

For me it was not clear that I was using the cache. Let me point out what I though:

Here I setup the problem

linsolve_1 = LinearProblem(A, b[1])
sol_1 = solve(linsolve_1)

Next, naively I just saw sol_1.cache being explicitly called, than I was expecting that I was copying it into a new object called linsolve_2.

linsolve_2 = LinearSolve.set_b(sol_1.cache, b[2])
sol_2 = solve(linsolve_2)

Instead, if I would see a !, like LinearSolve.set_b!, I would see an indication that I am changing sol_1.cache, and in ultimately, that my results would be tied to sol_1

4

You can do set_u to swap out the u vector.