From b8cfa09c6287138c97fd629621f70d4412f9ede8 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 11 Dec 2025 06:29:06 +0000 Subject: [PATCH 1/4] Add has_concretization support for AbstractSciMLOperator in sparse factorizations - Add `has_concretization` import from SciMLOperators to enable checking operator concretization support - Add `init_cacheval` methods for `AbstractSciMLOperator` types with sparse factorization algorithms: - KLUFactorization - UMFPACKFactorization - LUFactorization - CHOLMODFactorization - GenericFactorization - GenericLUFactorization - QRFactorization - NormalCholeskyFactorization - When an operator has `has_concretization`, it is converted to an AbstractMatrix and the appropriate sparse solver is used - Add tests for MatrixOperator and WOperator with sparse matrices - Bump SciMLOperators lower bound to v1.13 (WOperator was added in v1.13.0) --- Project.toml | 2 +- ext/LinearSolveSparseArraysExt.jl | 100 +++++++++++++++++++++++++++++- src/LinearSolve.jl | 2 +- test/basictests.jl | 60 +++++++++++++++++- 4 files changed, 160 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 6e652c6fa..b7f4990d8 100644 --- a/Project.toml +++ b/Project.toml @@ -129,7 +129,7 @@ Reexport = "1.2.2" SafeTestsets = "0.1" SciMLBase = "2.70" SciMLLogging = "1.3.1" -SciMLOperators = "1.7.1" +SciMLOperators = "1.13" Setfield = "1.1.1" SparseArrays = "1.10" Sparspak = "0.3.9" diff --git a/ext/LinearSolveSparseArraysExt.jl b/ext/LinearSolveSparseArraysExt.jl index 17249618c..fbcdb9dc8 100644 --- a/ext/LinearSolveSparseArraysExt.jl +++ b/ext/LinearSolveSparseArraysExt.jl @@ -5,7 +5,8 @@ using LinearSolve: LinearSolve, BLASELTYPES, pattern_changed, ArrayInterface, GenericLUFactorization, KLUFactorization, LUFactorization, NormalCholeskyFactorization, OperatorAssumptions, LinearVerbosity, - QRFactorization, RFLUFactorization, UMFPACKFactorization, solve + QRFactorization, RFLUFactorization, UMFPACKFactorization, solve, has_concretization +using SciMLOperators: AbstractSciMLOperator using ArrayInterface: ArrayInterface using LinearAlgebra: LinearAlgebra, I, Hermitian, Symmetric, cholesky, ldiv!, lu, lu!, QR using SparseArrays: SparseArrays, AbstractSparseArray, AbstractSparseMatrixCSC, @@ -274,6 +275,103 @@ function LinearSolve.init_cacheval( 0, 0, [Int32(1)], Int32[], Float64[])) end +# AbstractSciMLOperator handling for sparse factorizations +function LinearSolve.init_cacheval( + alg::KLUFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("KLUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + +function LinearSolve.init_cacheval( + alg::UMFPACKFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("UMFPACKFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + +function LinearSolve.init_cacheval( + alg::LUFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("LUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + +function LinearSolve.init_cacheval( + alg::CHOLMODFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("CHOLMODFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + +function LinearSolve.init_cacheval( + alg::GenericFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("GenericFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + +function LinearSolve.init_cacheval( + alg::GenericLUFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("GenericLUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + +function LinearSolve.init_cacheval( + alg::QRFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("QRFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + +function LinearSolve.init_cacheval( + alg::NormalCholeskyFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions) + if has_concretization(A) + return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions) + else + error("NormalCholeskyFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + end +end + function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::KLUFactorization; kwargs...) A = cache.A A = convert(AbstractMatrix, A) diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index a0e0c1345..7f45a3783 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -20,7 +20,7 @@ using SciMLBase: SciMLBase, LinearAliasSpecifier, AbstractSciMLOperator, init, solve!, reinit!, solve, ReturnCode, LinearProblem using SciMLOperators: SciMLOperators, AbstractSciMLOperator, IdentityOperator, MatrixOperator, - has_ldiv!, issquare + has_ldiv!, issquare, has_concretization using SciMLLogging: SciMLLogging, @SciMLMessage, verbosity_to_int, AbstractVerbositySpecifier, AbstractMessageLevel, AbstractVerbosityPreset, Silent, InfoLevel, WarnLevel, CustomLevel, None, Minimal, Standard, Detailed, All using Setfield: @set, @set! diff --git a/test/basictests.jl b/test/basictests.jl index dc7e0ad1f..018921ea1 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -1,5 +1,6 @@ using LinearSolve, LinearAlgebra, SparseArrays, MultiFloats, ForwardDiff -using SciMLOperators, RecursiveFactorization, Sparspak, FastLapackInterface +using SciMLOperators: SciMLOperators, MatrixOperator, FunctionOperator, WOperator +using RecursiveFactorization, Sparspak, FastLapackInterface using IterativeSolvers, KrylovKit, MKL_jll, KrylovPreconditioners using Test import CliqueTrees, Random @@ -561,6 +562,63 @@ end @test sol13.u ≈ sol23.u end + @testset "Operators with has_concretization" begin + n = 4 + Random.seed!(42) + A_sparse = sprand(n, n, 0.8) + I + b = rand(n) + + # Create a MatrixOperator wrapping the sparse matrix + A_op = MatrixOperator(A_sparse) + + prob_matrix = LinearProblem(A_sparse, b) + prob_operator = LinearProblem(A_op, b) + + # Test KLU with operator + sol_matrix = solve(prob_matrix, KLUFactorization()) + sol_operator = solve(prob_operator, KLUFactorization()) + @test sol_matrix.u ≈ sol_operator.u + + # Test UMFPACK with operator + sol_matrix = solve(prob_matrix, UMFPACKFactorization()) + sol_operator = solve(prob_operator, UMFPACKFactorization()) + @test sol_matrix.u ≈ sol_operator.u + + # Test LU with operator + sol_matrix = solve(prob_matrix, LUFactorization()) + sol_operator = solve(prob_operator, LUFactorization()) + @test sol_matrix.u ≈ sol_operator.u + + # Test WOperator with sparse Jacobian + n_w = 8 + M = sparse(I(n_w) * 1.0) + gamma = 1 / 2.0 + J = sprand(n_w, n_w, 0.5) + sparse(I(n_w) * 10.0) # Make it diagonally dominant + u = rand(n_w) + b_w = rand(n_w) + + W = WOperator{true}(M, gamma, J, u) + W_matrix = convert(AbstractMatrix, W) + + prob_woperator = LinearProblem(W, b_w) + prob_wmatrix = LinearProblem(W_matrix, b_w) + + # Test KLU with WOperator + sol_woperator = solve(prob_woperator, KLUFactorization()) + sol_wmatrix = solve(prob_wmatrix, KLUFactorization()) + @test sol_woperator.u ≈ sol_wmatrix.u + + # Test UMFPACK with WOperator + sol_woperator = solve(prob_woperator, UMFPACKFactorization()) + sol_wmatrix = solve(prob_wmatrix, UMFPACKFactorization()) + @test sol_woperator.u ≈ sol_wmatrix.u + + # Test LU with WOperator + sol_woperator = solve(prob_woperator, LUFactorization()) + sol_wmatrix = solve(prob_wmatrix, LUFactorization()) + @test sol_woperator.u ≈ sol_wmatrix.u + end + @testset "Solve Function" begin A1 = rand(n) |> Diagonal b1 = rand(n) From ad9386262c88b387bdde9b7db605c33fe248a2c7 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 11 Dec 2025 06:44:11 +0000 Subject: [PATCH 2/4] Add RequiresConcreteMatrixError struct for better error handling Replace plain error() calls with a custom exception type when operators that don't support concretization are passed to factorization algorithms. --- ext/LinearSolveSparseArraysExt.jl | 18 +++++++++--------- src/LinearSolve.jl | 26 ++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 9 deletions(-) diff --git a/ext/LinearSolveSparseArraysExt.jl b/ext/LinearSolveSparseArraysExt.jl index fbcdb9dc8..43572277d 100644 --- a/ext/LinearSolveSparseArraysExt.jl +++ b/ext/LinearSolveSparseArraysExt.jl @@ -4,7 +4,7 @@ using LinearSolve: LinearSolve, BLASELTYPES, pattern_changed, ArrayInterface, @get_cacheval, CHOLMODFactorization, GenericFactorization, GenericLUFactorization, KLUFactorization, LUFactorization, NormalCholeskyFactorization, - OperatorAssumptions, LinearVerbosity, + OperatorAssumptions, LinearVerbosity, RequiresConcreteMatrixError, QRFactorization, RFLUFactorization, UMFPACKFactorization, solve, has_concretization using SciMLOperators: AbstractSciMLOperator using ArrayInterface: ArrayInterface @@ -284,7 +284,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("KLUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("KLUFactorization")) end end @@ -296,7 +296,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("UMFPACKFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("UMFPACKFactorization")) end end @@ -308,7 +308,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("LUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("LUFactorization")) end end @@ -320,7 +320,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("CHOLMODFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("CHOLMODFactorization")) end end @@ -332,7 +332,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("GenericFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("GenericFactorization")) end end @@ -344,7 +344,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("GenericLUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("GenericLUFactorization")) end end @@ -356,7 +356,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("QRFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("QRFactorization")) end end @@ -368,7 +368,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - error("NormalCholeskyFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") + throw(RequiresConcreteMatrixError("NormalCholeskyFactorization")) end end diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 7f45a3783..63f79394a 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -230,6 +230,30 @@ for integrating custom algorithms or simple solve strategies. """ abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end +# Errors + +""" + RequiresConcreteMatrixError(alg::String) + +Error thrown when a factorization algorithm that requires a concrete matrix +receives an operator that does not support concretization (i.e., `has_concretization(A)` returns `false`). + +## Example + +```julia +throw(RequiresConcreteMatrixError("KLUFactorization")) +``` + +This error suggests using a Krylov method instead, which can work with abstract operators. +""" +struct RequiresConcreteMatrixError <: Exception + alg::String +end + +function Base.showerror(io::IO, e::RequiresConcreteMatrixError) + print(io, "$(e.alg) requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") +end + # Traits """ @@ -519,6 +543,8 @@ export MetalOffload32MixedLUFactorization export OperatorAssumptions, OperatorCondition +export RequiresConcreteMatrixError + export LinearSolveAdjoint export LinearVerbosity From b3acd983ef582d42a17ce6b70dad12c42ac0c735 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 11 Dec 2025 09:17:13 +0000 Subject: [PATCH 3/4] Fix has_concretization import in SparseArrays extension Import has_concretization directly from SciMLOperators instead of LinearSolve, since LinearSolve doesn't re-export it. --- ext/LinearSolveSparseArraysExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/LinearSolveSparseArraysExt.jl b/ext/LinearSolveSparseArraysExt.jl index 43572277d..bef19c006 100644 --- a/ext/LinearSolveSparseArraysExt.jl +++ b/ext/LinearSolveSparseArraysExt.jl @@ -5,8 +5,8 @@ using LinearSolve: LinearSolve, BLASELTYPES, pattern_changed, ArrayInterface, GenericLUFactorization, KLUFactorization, LUFactorization, NormalCholeskyFactorization, OperatorAssumptions, LinearVerbosity, RequiresConcreteMatrixError, - QRFactorization, RFLUFactorization, UMFPACKFactorization, solve, has_concretization -using SciMLOperators: AbstractSciMLOperator + QRFactorization, RFLUFactorization, UMFPACKFactorization, solve +using SciMLOperators: AbstractSciMLOperator, has_concretization using ArrayInterface: ArrayInterface using LinearAlgebra: LinearAlgebra, I, Hermitian, Symmetric, cholesky, ldiv!, lu, lu!, QR using SparseArrays: SparseArrays, AbstractSparseArray, AbstractSparseMatrixCSC, From f7bd1fa45bc8db2d00adf522e0a5d0eab93d4a59 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 11 Dec 2025 09:27:58 +0000 Subject: [PATCH 4/4] Revert error throwing to return nothing for default algorithm support The init_cacheval methods need to return nothing when concretization is not supported, so that the default algorithm selection can fall back to other methods like Krylov solvers. --- ext/LinearSolveSparseArraysExt.jl | 18 +++++++++--------- src/LinearSolve.jl | 26 -------------------------- 2 files changed, 9 insertions(+), 35 deletions(-) diff --git a/ext/LinearSolveSparseArraysExt.jl b/ext/LinearSolveSparseArraysExt.jl index bef19c006..ae3fc27f2 100644 --- a/ext/LinearSolveSparseArraysExt.jl +++ b/ext/LinearSolveSparseArraysExt.jl @@ -4,7 +4,7 @@ using LinearSolve: LinearSolve, BLASELTYPES, pattern_changed, ArrayInterface, @get_cacheval, CHOLMODFactorization, GenericFactorization, GenericLUFactorization, KLUFactorization, LUFactorization, NormalCholeskyFactorization, - OperatorAssumptions, LinearVerbosity, RequiresConcreteMatrixError, + OperatorAssumptions, LinearVerbosity, QRFactorization, RFLUFactorization, UMFPACKFactorization, solve using SciMLOperators: AbstractSciMLOperator, has_concretization using ArrayInterface: ArrayInterface @@ -284,7 +284,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("KLUFactorization")) + nothing end end @@ -296,7 +296,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("UMFPACKFactorization")) + nothing end end @@ -308,7 +308,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("LUFactorization")) + nothing end end @@ -320,7 +320,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("CHOLMODFactorization")) + nothing end end @@ -332,7 +332,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("GenericFactorization")) + nothing end end @@ -344,7 +344,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("GenericLUFactorization")) + nothing end end @@ -356,7 +356,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("QRFactorization")) + nothing end end @@ -368,7 +368,7 @@ function LinearSolve.init_cacheval( return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters, abstol, reltol, verbose, assumptions) else - throw(RequiresConcreteMatrixError("NormalCholeskyFactorization")) + nothing end end diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 63f79394a..7f45a3783 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -230,30 +230,6 @@ for integrating custom algorithms or simple solve strategies. """ abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end -# Errors - -""" - RequiresConcreteMatrixError(alg::String) - -Error thrown when a factorization algorithm that requires a concrete matrix -receives an operator that does not support concretization (i.e., `has_concretization(A)` returns `false`). - -## Example - -```julia -throw(RequiresConcreteMatrixError("KLUFactorization")) -``` - -This error suggests using a Krylov method instead, which can work with abstract operators. -""" -struct RequiresConcreteMatrixError <: Exception - alg::String -end - -function Base.showerror(io::IO, e::RequiresConcreteMatrixError) - print(io, "$(e.alg) requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.") -end - # Traits """ @@ -543,8 +519,6 @@ export MetalOffload32MixedLUFactorization export OperatorAssumptions, OperatorCondition -export RequiresConcreteMatrixError - export LinearSolveAdjoint export LinearVerbosity