Drop phase of Factor?

How can I reserve the symb. phase?

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

1 Like

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/

1 Like

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).

1 Like

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)

1 Like

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.

1 Like

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

1 Like

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.