How can I reserve the symb. phase?

(Posted and answered here: How to drop the Symbolic phase? · Issue #6 · JuliaSparse/KLU.jl · GitHub)

Note that if you use LinearSolve.jl in its caching mode, it will cache the symbolic phase automatically.

http://linearsolve.sciml.ai/dev/tutorials/caching_interface/

hanks for your note.

Can I use KLU in the caching mode of `LinearSolve.jl`

?

if so, could you please guide me on my previous MWE?

Just do exactly what it shows in the tutorial I linked.

http://linearsolve.sciml.ai/dev/tutorials/caching_interface/

```
using LinearSolve, SparseArrays
n = 100
A = sprand(n,n,0.1)
b1 = rand(n); b2 = rand(n)
prob = LinearProblem(A, b1)
linsolve = init(prob,KLUFactorization())
sol1 = solve(linsolve) # Does symbolic factorization, then numeric factorization, then backsolve
sol1.u
linsolve = LinearSolve.set_b(sol1.cache,b2)
sol2 = solve(linsolve) # only does backsolve
sol2.u
A.nzval .= rand(length(A.nzval)) # new sparse matrix, same sparsity pattern
linsolve = LinearSolve.set_A(sol2.cache,A)
sol3 = solve(linsolve) # new numerical factorization and backsolve, no new symbolic factorization
```

Note that as described in the documentation:

http://linearsolve.sciml.ai/dev/solvers/solvers/#SuiteSparse.jl

`KLUFactorization(reuse_symbolic=true)`

is the default, so changing `A`

will reuse the symbolic factorization (thus assuming the sparsity pattern is the same), and one needs to do `KLUFactorization(;reuse_symbolic=false)`

to make it symbolic refactorize with changes to `A`

.

So the behavior you want is the default (for clear reasons).

Run it with a different example. I updated it so it’s a sparse matrix which is (most likely) non-singular (though there’s a change any random sparse matrix is singular)

Drop the `if t == 1`

thing, I’m not sure what you’re trying to do with that. It’s just, if you make a new LinearProblem, i.e. everything is brand new, then it takes a new symbolic factorization (of course, it has to), numeric factorization, and backsolve. If you only change A, then it does a new numeric factorization and backsolve (if the sparsity pattern is the same). If you only change b, then it only does the backsolve. I.e., it only does the computations it has to do.

Thanks for your feedback and I have dropped the `t == 1`

thing. So, is my structure below correct?

```
A = []
result = []
linsolve = []
for t in 1:100 # Time-loop
# calculate "b"
# calculate "A"
#if (A is changed in its pattern)
prob = LinearProblem(A, b)
linsolve = init(prob,KLUFactorization())
# end
#if (A is changed in its value only)
linsolve = LinearSolve.set_A(result.cache, A)
# end
linsolve = LinearSolve.set_b(result.cache, b)
result = solve(linsolve)
X = result.u
end # for t in 1:100
```

Yes

Thank you for the confirmation.

However, following this procedure if cashing interface, I didnt have the results as without cashing interface.

Post an MWE

@ChrisRackauckas Here it is:

```
function s1()
B = ones(10);
A = spzeros(10,10);
A[diagind(A)] .= 1;
t_start = 1;
t_end = 10;
global R_Record1 = zeros(length(B),size(t_start:t_end,1));
result = []
linsolve = []
A_old = deepcopy(A);
iter = 0;
for t in t_start:t_end
iter +=1;
B[t] = t;
if iseven(t)
A[t,t] = t;
else
A[end,t] = 10;
end
if A_old != A
if t == 1
prob = LinearProblem(A, B)
linsolve = init(prob,KLUFactorization())
result = solve(linsolve)
end
if iseven(t)
linsolve = LinearSolve.set_A(result.cache, A)
else
prob = LinearProblem(A, B)
linsolve = init(prob,KLUFactorization())
result = solve(linsolve)
end
A_old = deepcopy(A);
end
linsolve = LinearSolve.set_b(result.cache, B)
result = solve(linsolve)
R_Record1[:,iter] = result.u;
end # for t in 1:100
end # function s1()
function s2()
B = ones(10);
A = spzeros(10,10);
A[diagind(A)] .= 1;
t_start = 1;
t_end = 10;
global R_Record2 = zeros(length(B),size(t_start:t_end,1));
result = [];
linsolve = [];
iter = 0;
for t in t_start:t_end
iter +=1;
B[t] = t;
if iseven(t)
A[t,t] = t;
else
A[end,t] = 10;
end
prob = LinearProblem(A, B)
linsolve = init(prob,KLUFactorization())
#linsolve = LinearSolve.set_A(result.cache, A)
#linsolve = LinearSolve.set_b(result.cache, B)
result = solve(linsolve)
R_Record2[:,iter] = result.u;
end # for t in 1:100
end # function s2()
s1()
s2()
julia> R_Record1 == R_Record2
false
```

Let me take a look at this tonight. It’s possible there’s an error, or that I need to check the matrices more on input. The refactor capability was finicky at first.