Cache issues with sequential calls to solve

[yosinlpet]

Hello,

I'm trying to solve a constrained linear problem $Ax = b$ under the set of constraints $C^Tx = d$ as a step in a fluid simulation. Using Lagrange optimization, one finds the least-square approximated solution as the solution of the linear system:

$[2A^TA \quad C; C^T \quad 0] \cdot [x; \quad \lambda] = [2A^Tb;\quad d]$

This system need be solved for different RHS so using caching abilities of LinearSolve.jl seems adequate. My code is therefore something like that:

function myfunc(time)
     # some code

     A=somedensematrix(time, args)
     b1, b2, b3=somerhs(time, args)
     C, d1, d2, d3=someconstraints(nconstr)
     
     Aᵀ = 2A'                                                                    
     Ã = vcat(hcat(Aᵀ*A, C), hcat(C', zeros(nconstr, nconstr)))                              
                                                                                     
     display(Ã)                                                                                  
                                                                                     
     prob = LinearProblem(Ã, vcat(Aᵀ*b1, d1))                          
     linsolve = init(prob)                                                       
     sol1 = solve(prob)            
                                                                                                                          
     linsolve = LinearSolve.set_b(sol1.cache, vcat(Aᵀ*b2, d2))                   
     sol2 = solve(linsolve)                                                      
                                                                                  
     linsolve = LinearSolve.set_b(sol1.cache, vcat(Aᵀ*b3, d3))                   
     sol3 = solve(linsolve)
end

t = 0
while t < T
    myfunc(t)
    t += dt
end

which seems to work for a little bit and then, I run into a

ERROR: LoadError: LinearAlgebra.SingularException(1)

but what tickles me is that display(Ã) prompts a matrix filled with zeros, as though there were some memory leak(?) because there is no way somedensematrix() could return a zero matrix.

I do not get the issue when I create 3 separate LinearProblems for each RHS. Am I missing something?

I'm running julia 1.8.1 on MacOS Monterey 12.5.1

[ChrisRackauckas]

You did not give a runnable code and I could not reproduce the error. When I replaced the undefined functions with functions that give dense matrices, it worked just fine. Here's what I ran:

using LinearSolve

function myfunc(time)
    # some code

    A=rand(2,2)
    b1, b2, b3= rand(2), rand(2), rand(2)
    C, d1, d2, d3= rand(2,2), rand(2), rand(2), rand(2)
    
    Aᵀ = 2A'                                                                    
    Ã = vcat(hcat(Aᵀ*A, C), hcat(C', zeros(nconstr, nconstr)))                              
                                                                                    
    display(Ã)                                                                                  
                                                                                    
    prob = LinearProblem(Ã, vcat(Aᵀ*b1, d1))                          
    linsolve = init(prob)                                                       
    sol1 = solve(prob)            
                                                                                                                         
    linsolve = LinearSolve.set_b(sol1.cache, vcat(Aᵀ*b2, d2))                   
    sol2 = solve(linsolve)                                                      
                                                                                 
    linsolve = LinearSolve.set_b(sol1.cache, vcat(Aᵀ*b3, d3))                   
    sol3 = solve(linsolve)
end

t = 0
nconstr = 2

for i in 1:1000
    myfunc(0.0)
end

#=
4×4 Matrix{Float64}:
 2.22186    1.70077    0.0891254  0.829554
 1.70077    1.91723    0.953059   0.0255396
 0.0891254  0.953059   0.0        0.0
 0.829554   0.0255396  0.0        0.0

4×4 Matrix{Float64}:
 2.16978   1.68399   0.271874  0.725556
 1.68399   1.51925   0.395511  0.123642
 0.271874  0.395511  0.0       0.0
 0.725556  0.123642  0.0       0.0

4×4 Matrix{Float64}:
 0.0201618  0.188836  0.212852  0.992754
 0.188836   2.30555   0.812849  0.843758
 0.212852   0.812849  0.0       0.0
 0.992754   0.843758  0.0       0.0

4×4 Matrix{Float64}:
 2.0505      1.5875     0.00971727  0.348122
 1.5875      1.47593    0.0600411   0.645411
 0.00971727  0.0600411  0.0         0.0
 0.348122    0.645411   0.0         0.0

4×4 Matrix{Float64}:
 1.95893   1.89333   0.90924   0.272689
 1.89333   1.90412   0.777559  0.860229
 0.90924   0.777559  0.0       0.0
 0.272689  0.860229  0.0       0.0

4×4 Matrix{Float64}:
 0.344479  0.674005  0.715219  0.558816
 0.674005  1.33954   0.859303  0.644183
 0.715219  0.859303  0.0       0.0
 0.558816  0.644183  0.0       0.0

4×4 Matrix{Float64}:
 1.228     1.15063  0.951111  0.618644
 1.15063   1.07818  0.56325   0.86681
 0.951111  0.56325  0.0       0.0
 0.618644  0.86681  0.0       0.0

4×4 Matrix{Float64}:
 1.65317    0.763111  0.389893  0.0952648
 0.763111   0.355675  0.967553  0.971726
 0.389893   0.967553  0.0       0.0
 0.0952648  0.971726  0.0       0.0

4×4 Matrix{Float64}:
 1.58617   0.889494  0.787629  0.699808
 0.889494  0.623557  0.26398   0.482757
 0.787629  0.26398   0.0       0.0
 0.699808  0.482757  0.0       0.0

4×4 Matrix{Float64}:
 1.82908   0.939493   0.820719  0.314143
 0.939493  0.4827     0.533588  0.0440287
 0.820719  0.533588   0.0       0.0
 0.314143  0.0440287  0.0       0.0

4×4 Matrix{Float64}:
 1.46614   0.687362  0.123804  0.747498
 0.687362  0.492805  0.334168  0.681448
 0.123804  0.334168  0.0       0.0
 0.747498  0.681448  0.0       0.0

4×4 Matrix{Float64}:
 0.310385  0.914029  0.384903  0.974926
 0.914029  2.80172   0.348352  0.718685
 0.384903  0.348352  0.0       0.0
 0.974926  0.718685  0.0       0.0

4×4 Matrix{Float64}:
 1.70027   1.21622   0.280238  0.442055
 1.21622   1.04509   0.464349  0.176051
 0.280238  0.464349  0.0       0.0
 0.442055  0.176051  0.0       0.0
=#

Please share a runnable code that can be copy-pasted and display the error.