LinearSolve and SparseArrays: reuse factorization for multiple right-hand side vectors?

Hi everyone,

I am implementing several Finite Element procedures in Julia and I am trying to switch all the implementation to use LinearSolve. I would like to know if it is possible to reuse a symbolic factorization of the coefficient matrix when the right-hand side is also changing, and I was hoping someone could help me out.
I have a fairly complex example in which I slice a large sparse matrix (using @view) and solve the LinearProblem for given right-hand side vectors. For simplicity, this is an equivalent MWE:

using LinearAlgebra
using LinearSolve
using SparseArrays

A = sparse([1. 0 0 2 0 0; 0 3 1 0 0 0; 0 1 1 0 0 1; 2 0 0 5 0 0; 0 0 0 0 1 0; 0 0 1 0 0 2])
id_dof = BitVector([1, 0, 1, 1, 0, 1])
A_11 = @view A[id_dof, id_dof]
b_1 = [3., 2., 6., 1.]
problem = LinearProblem(A_11, b_1)
solution = solve(problem)

If anyone is curious, we generally construct a large, sparse stiffness matrix (which is usually preconditioned), and then solve the problem over a sub-domain (the free nodes). Up to here, everything works great. After each iteration , however, we re-compute the right-hand side, update A (which does not change its sparsity pattern), and perform the next iteration / step.

Reading the LinearSolve Documentation and watching various tutorials, I could only find examples for re-setting either A or b using the cache of the previously solved step. I tried updating both in series, but obviously it doesn’t work:

A[1,1] = 12 # Update A
b_1 = [6., 3. 1.5, 0.5]
my_problem = LinearSolve.set_A(solution.cache, A_11) # It works up to here
my_problem = LinearSolve.set_b(solution.cache, A_11) # As expected, this reverts A to the previous instance

Since the sparsity pattern of A does not change, I was wondering if it would be worthwhile (in terms of computational speed) to reuse such information in the next step or not, and if that would be possible using LinearSolve.

Hi Alessandro, I also had this question and was wondering if you made progress on this? I noted that the SuiteSparse interface has a reuse_symbolic keyword that is intended for reuse of the sparsity pattern.

1 Like

Check out GitHub - PetrKryslUCSD/Sparspak.jl: Solution of large systems of linear algebraic equations : reusing the structure of the matrix and multiple right hand sides are no problem.