`LinearSolve.solve!` throws error when using `KLUFactorization` in `LinearSolve.init`

[mariekeosw]

Hi there,
I was working with a linear system defined by an ExtendableSparseMatrix, which I wanted to solve using LinearSolve.jl. However, when I tried to solve the system using the cache, I received the following error. The code and error message are documented below.

Note that using linsolve = init(prob) or linsolve = init(prob, UMFPACKFactorization()) runs successfully.

function minimal_example(N; usecache=true)

	b = zeros(N)
	A = ExtendableSparseMatrix(N,N)

	for i in 1:N
		A[i,i] = 1
	end
	flush!(A)
	
	prob = LinearProblem(A, b)
	
	if usecache
		linsolve = LinearSolve.init(prob, KLUFactorization())
		LinearSolve.solve!(linsolve)
	else
		LinearSolve.solve(prob)
	end
	
	return "success return"
end;

Resulting error when using usecache=true:

ERROR: DimensionMismatch: 
Stacktrace:
 [1] klu!(K::LinearSolveSparseArraysExt.KLU.KLUFactorization{…}, nzval::Vector{…}; check::Bool, allowsingular::Bool)
   @ LinearSolveSparseArraysExt.KLU C:\Users\marie\.julia\packages\LinearSolve\wrTJy\src\KLU\klu.jl:640
 [2] klu!
   @ C:\Users\marie\.julia\packages\LinearSolve\wrTJy\src\KLU\klu.jl:638 [inlined]
 [3] solve!(cache::LinearSolve.LinearCache{…}, alg::KLUFactorization; kwargs::@Kwargs{})
   @ LinearSolveSparseArraysExt C:\Users\marie\.julia\packages\LinearSolve\wrTJy\ext\LinearSolveSparseArraysExt.jl:247
 [4] solve!
   @ C:\Users\marie\.julia\packages\LinearSolve\wrTJy\ext\LinearSolveSparseArraysExt.jl:235 [inlined]
 [5] #solve!#11
   @ C:\Users\marie\.julia\packages\LinearSolve\wrTJy\src\common.jl:299 [inlined]
 [6] solve!(::LinearSolve.LinearCache{…})
   @ LinearSolve C:\Users\marie\.julia\packages\LinearSolve\wrTJy\src\common.jl:298
 [7] minimal_example(N::Int64; usecache::Bool)
   @ Main c:\Users\marie\Desktop\WIAS\04_Minimal-working-examples\Minimal-working-example_Compatibility-ExtendableSparse-LinearSolve.jl:20
 [8] top-level scope
   @ c:\Users\marie\Desktop\WIAS\04_Minimal-working-examples\Minimal-working-example_Compatibility-ExtendableSparse-LinearSolve.jl:30
Some type information was truncated. Use `show(err)` to see complete types

[j-fu]

Ok I put this onto my list for after my vacancies...

[jpthiele]

I had a quick look and get the same error.
However, when I use the SparseArrays.SparseMatrixCSC that is definitely fully populated after flush! everything works fine (somehow).
So a workaround while we figure this out would be:

prob = LinearProblem(A.cscmatrix,b)

[j-fu]

This is an issue with LinearSolve, here is another MWE which uses the (not completely exported) AbstractSparseMatrixCSC interface (which is followed by ExtendableSparse):

using LinearAlgebra, SparseArrays, LinearSolve

using LinearSolve
using Pardiso
using Sparspak


struct MySparseMatrixCSC{Tv, Ti} <: SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}
    csc::SparseMatrixCSC{Tv, Ti}
end

Base.size(m::MySparseMatrixCSC) = size(m.csc)
SparseArrays.getcolptr(m::MySparseMatrixCSC) = SparseArrays.getcolptr(m.csc)
SparseArrays.rowvals(m::MySparseMatrixCSC) = SparseArrays.rowvals(m.csc)
SparseArrays.nonzeros(m::MySparseMatrixCSC) = SparseArrays.nonzeros(m.csc)

function main(; N = 10000, solver = UMFPACKFactorization())
    A = spdiagm(1 => -ones(N - 1), 0 => fill(10.0, N), -1 => -ones(N - 1))
    u0 = ones(size(A, 2))
    b = A * u0
    MyA = MySparseMatrixCSC(A)

    prob = LinearProblem(MyA, b)

    linsolve = LinearSolve.init(prob, solver)
    u = LinearSolve.solve!(linsolve)

    @show n = norm(u - u0, Inf)
    if n < 1.0e-12
        println("$(typeof(solver)) ok")
    end
    return nothing
end


main()
main(solver = MKLPardisoFactorize())
main(solver = SparspakFactorization())
main(solver = KLUFactorization())

So this should become an issue with LinearSolve, or, better, a PR which adds this
as a test. I already added some MySparseMatrixCSC stuff there, this could be amended:
https://github.com/SciML/LinearSolve.jl/blob/00d4764bef7251997282afea4c3ccf23475e5db3/test/basictests.jl#L580

I can do this next wee after return from my vacancy trip.

[j-fu]

See SciML/LinearSolve.jl#737