LinearSolve missing methods for do_factorization on SparspakFactorization

[ErikQQY]

CI in NonlinearSolve indicates the new version of LinearSolve is breaking NonlinearSolve, BoundaryValueDiffEq, etc, see https://github.com/SciML/NonlinearSolve.jl/actions/runs/13265094215/job/37031436009

Possibly related to #573 and caused by #570

[ErikQQY]

OK, it turns out LinearSolve is failing to solve sparse linear problems with version 3.1.0

using LinearSolve, SparseArrays, LinearAlgebra
n = 8
A = Matrix(I, n, n)/1
b = rand(n);
x = zero(b);
prob = LinearProblem(sparse(A), b; u0 = x)
sol = solve(prob)

Stacktrace:

ERROR: MethodError: no method matching do_factorization(::SparspakFactorization, ::SparseMatrixCSC{Float64, Int64}, ::Vector{Float64}, ::Vector{Float64})
The function `do_factorization` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  do_factorization(::BunchKaufmanFactorization, ::Any, ::Any, ::Any)
   @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/factorization.jl:439
  do_factorization(::QRFactorization, ::Any, ::Any, ::Any)
   @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/factorization.jl:214
  do_factorization(::GenericLUFactorization, ::Any, ::Any, ::Any)
   @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/factorization.jl:147
  ...

Stacktrace:
  [1] init_cacheval(alg::SparspakFactorization, A::SparseMatrixCSC{…}, b::Vector{…}, u::Vector{…}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::Bool, assumptions::OperatorAssumptions{…})
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/factorization.jl:29
  [2] macro expansion
    @ ~/.julia/packages/LinearSolve/SFMdw/src/default.jl:308 [inlined]
  [3] init_cacheval(alg::LinearSolve.DefaultLinearSolver, A::SparseMatrixCSC{…}, b::Vector{…}, u::Vector{…}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::Bool, assump::OperatorAssumptions{…})
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/default.jl:291
  [4] init(::LinearProblem{…}, ::LinearSolve.DefaultLinearSolver; alias::LinearAliasSpecifier, abstol::Float64, reltol::Float64, maxiters::Int64, verbose::Bool, Pl::Nothing, Pr::Nothing, assumptions::OperatorAssumptions{…}, sensealg::LinearSolveAdjoint{…}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:239
  [5] init(::LinearProblem{Vector{Float64}, true, SparseMatrixCSC{Float64, Int64}, Vector{Float64}, SciMLBase.NullParameters, @Kwargs{}}, ::LinearSolve.DefaultLinearSolver)
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:140
  [6] solve(::LinearProblem{Vector{…}, true, SparseMatrixCSC{…}, Vector{…}, SciMLBase.NullParameters, @Kwargs{}}, ::LinearSolve.DefaultLinearSolver; kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:303
  [7] solve(::LinearProblem{Vector{Float64}, true, SparseMatrixCSC{Float64, Int64}, Vector{Float64}, SciMLBase.NullParameters, @Kwargs{}}, ::LinearSolve.DefaultLinearSolver)
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:301
  [8] solve(::LinearProblem{Vector{…}, true, SparseMatrixCSC{…}, Vector{…}, SciMLBase.NullParameters, @Kwargs{}}, ::Nothing; assump::OperatorAssumptions{Bool}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:298
  [9] solve(::LinearProblem{Vector{Float64}, true, SparseMatrixCSC{Float64, Int64}, Vector{Float64}, SciMLBase.NullParameters, @Kwargs{}}, ::Nothing)
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:296
 [10] solve(::LinearProblem{Vector{Float64}, true, SparseMatrixCSC{Float64, Int64}, Vector{Float64}, SciMLBase.NullParameters, @Kwargs{}}; kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:293
 [11] solve(::LinearProblem{Vector{Float64}, true, SparseMatrixCSC{Float64, Int64}, Vector{Float64}, SciMLBase.NullParameters, @Kwargs{}})
    @ LinearSolve ~/.julia/packages/LinearSolve/SFMdw/src/common.jl:292
 [12] top-level scope
    @ REPL[19]:1
Some type information was truncated. Use `show(err)` to see complete types.

Tested on an empty environment with

(jl_Aapfoi) pkg> st
Status `/private/var/folders/zy/dgghnxqx0n15lmdlpjrhk6400000gn/T/jl_Aapfoi/Project.toml`
  [7ed4a6bd] LinearSolve v3.1.0

and versioninfo()

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M3
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m3)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)

[j-fu]

If you add "using Sparspak", the example should work. However, in this case the default should be to use UMFPACK/Suitesparse via SparseArrays as the default for Float64 sparse matrices. So indeed there is something broken.

[ErikQQY]

Yeah, since SparseArrays and Sparspak are moved to extension to avoid giant precompiling time, we would have to first using them, I will fix the downstream errors then.

[ErikQQY]

Seems b9668a4 will solve the problem on [email protected], this issue can be safely closed then

[Ickaser]

Ran into this today after updating package updates--had to do some searching to find out I just needed using Sparspak. Might be good to document this somewhere?

[ChrisRackauckas]

Can I get an MWE? It's supposed to be throwing a better error, but for some reason it's not being hit in some people's examples and I haven't been able to find out why.

[Ickaser]

This runs into the error for me:

using SparseArrays
using LinearSolve

function setupsolve_pde(n)
    dx = 1/n
    dx2 = dx^-2
    vals = Vector{Float64}(undef, 0)
    cols = Vector{Int}(undef, 0)
    rows = Vector{Int}(undef, 0)
    for i in 1:n
        if i != 1
            push!(vals, dx2)
            push!(cols, i-1)
            push!(rows, i)
        end
        push!(vals, -2dx2)
        push!(cols, i)
        push!(rows, i)
        if i != n
            push!(vals, dx2)
            push!(cols, i+1)
            push!(rows, i)
        end
    end
    mat = sparse(rows, cols, vals, n, n)
    rhs = zeros(n)
    rhs[begin] = rhs[end] = -2
    prob = LinearProblem(mat, rhs)
    sol = solve(prob, SparspakFactorization()).u
end

setupsolve_pde(10)

LinearSolve 3.4.0, SparseArrays, 1.11.0; using Sparspak with Sparspak 0.3.9 fixes.
Versioninfo:

Julia Version 1.11.1
Commit 8f5b7ca12a (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, tigerlake)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_PYTHONCALL_EXE = C:\Users\iwheeler\Miniconda3\python.exe

[Ickaser]

Looking at #580 , my example doesn't run into the informative error because I explicitly specified SparspakFactorization() rather than expecting a default method. So providing a better error might require some of the same gymnastics as SciML/OrdinaryDiffEq.jl#2554 , although in this case we can explicitly check for the extension rather than check if the method exists.