Skip to content

Commit b8cfa09

Browse files
claudeclaude
committed
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)
1 parent 25de27c commit b8cfa09

File tree

4 files changed

+160
-4
lines changed

4 files changed

+160
-4
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -129,7 +129,7 @@ Reexport = "1.2.2"
129129
SafeTestsets = "0.1"
130130
SciMLBase = "2.70"
131131
SciMLLogging = "1.3.1"
132-
SciMLOperators = "1.7.1"
132+
SciMLOperators = "1.13"
133133
Setfield = "1.1.1"
134134
SparseArrays = "1.10"
135135
Sparspak = "0.3.9"

‎ext/LinearSolveSparseArraysExt.jl‎

Lines changed: 99 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,8 @@ using LinearSolve: LinearSolve, BLASELTYPES, pattern_changed, ArrayInterface,
55
GenericLUFactorization,
66
KLUFactorization, LUFactorization, NormalCholeskyFactorization,
77
OperatorAssumptions, LinearVerbosity,
8-
QRFactorization, RFLUFactorization, UMFPACKFactorization, solve
8+
QRFactorization, RFLUFactorization, UMFPACKFactorization, solve, has_concretization
9+
using SciMLOperators: AbstractSciMLOperator
910
using ArrayInterface: ArrayInterface
1011
using LinearAlgebra: LinearAlgebra, I, Hermitian, Symmetric, cholesky, ldiv!, lu, lu!, QR
1112
using SparseArrays: SparseArrays, AbstractSparseArray, AbstractSparseMatrixCSC,
@@ -274,6 +275,103 @@ function LinearSolve.init_cacheval(
274275
0, 0, [Int32(1)], Int32[], Float64[]))
275276
end
276277

278+
# AbstractSciMLOperator handling for sparse factorizations
279+
function LinearSolve.init_cacheval(
280+
alg::KLUFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
281+
maxiters::Int, abstol, reltol,
282+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
283+
if has_concretization(A)
284+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
285+
maxiters, abstol, reltol, verbose, assumptions)
286+
else
287+
error("KLUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
288+
end
289+
end
290+
291+
function LinearSolve.init_cacheval(
292+
alg::UMFPACKFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
293+
maxiters::Int, abstol, reltol,
294+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
295+
if has_concretization(A)
296+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
297+
maxiters, abstol, reltol, verbose, assumptions)
298+
else
299+
error("UMFPACKFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
300+
end
301+
end
302+
303+
function LinearSolve.init_cacheval(
304+
alg::LUFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
305+
maxiters::Int, abstol, reltol,
306+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
307+
if has_concretization(A)
308+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
309+
maxiters, abstol, reltol, verbose, assumptions)
310+
else
311+
error("LUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
312+
end
313+
end
314+
315+
function LinearSolve.init_cacheval(
316+
alg::CHOLMODFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
317+
maxiters::Int, abstol, reltol,
318+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
319+
if has_concretization(A)
320+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
321+
maxiters, abstol, reltol, verbose, assumptions)
322+
else
323+
error("CHOLMODFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
324+
end
325+
end
326+
327+
function LinearSolve.init_cacheval(
328+
alg::GenericFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
329+
maxiters::Int, abstol, reltol,
330+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
331+
if has_concretization(A)
332+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
333+
maxiters, abstol, reltol, verbose, assumptions)
334+
else
335+
error("GenericFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
336+
end
337+
end
338+
339+
function LinearSolve.init_cacheval(
340+
alg::GenericLUFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
341+
maxiters::Int, abstol, reltol,
342+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
343+
if has_concretization(A)
344+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
345+
maxiters, abstol, reltol, verbose, assumptions)
346+
else
347+
error("GenericLUFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
348+
end
349+
end
350+
351+
function LinearSolve.init_cacheval(
352+
alg::QRFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
353+
maxiters::Int, abstol, reltol,
354+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
355+
if has_concretization(A)
356+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
357+
maxiters, abstol, reltol, verbose, assumptions)
358+
else
359+
error("QRFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
360+
end
361+
end
362+
363+
function LinearSolve.init_cacheval(
364+
alg::NormalCholeskyFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
365+
maxiters::Int, abstol, reltol,
366+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
367+
if has_concretization(A)
368+
return LinearSolve.init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
369+
maxiters, abstol, reltol, verbose, assumptions)
370+
else
371+
error("NormalCholeskyFactorization requires a concrete matrix. The provided operator does not support concretization. Use a Krylov method instead.")
372+
end
373+
end
374+
277375
function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::KLUFactorization; kwargs...)
278376
A = cache.A
279377
A = convert(AbstractMatrix, A)

‎src/LinearSolve.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ using SciMLBase: SciMLBase, LinearAliasSpecifier, AbstractSciMLOperator,
2020
init, solve!, reinit!, solve, ReturnCode, LinearProblem
2121
using SciMLOperators: SciMLOperators, AbstractSciMLOperator, IdentityOperator,
2222
MatrixOperator,
23-
has_ldiv!, issquare
23+
has_ldiv!, issquare, has_concretization
2424
using SciMLLogging: SciMLLogging, @SciMLMessage, verbosity_to_int, AbstractVerbositySpecifier, AbstractMessageLevel, AbstractVerbosityPreset,
2525
Silent, InfoLevel, WarnLevel, CustomLevel, None, Minimal, Standard, Detailed, All
2626
using Setfield: @set, @set!

‎test/basictests.jl‎

Lines changed: 59 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
using LinearSolve, LinearAlgebra, SparseArrays, MultiFloats, ForwardDiff
2-
using SciMLOperators, RecursiveFactorization, Sparspak, FastLapackInterface
2+
using SciMLOperators: SciMLOperators, MatrixOperator, FunctionOperator, WOperator
3+
using RecursiveFactorization, Sparspak, FastLapackInterface
34
using IterativeSolvers, KrylovKit, MKL_jll, KrylovPreconditioners
45
using Test
56
import CliqueTrees, Random
@@ -561,6 +562,63 @@ end
561562
@test sol13.u ≈ sol23.u
562563
end
563564

565+
@testset "Operators with has_concretization" begin
566+
n = 4
567+
Random.seed!(42)
568+
A_sparse = sprand(n, n, 0.8) + I
569+
b = rand(n)
570+
571+
# Create a MatrixOperator wrapping the sparse matrix
572+
A_op = MatrixOperator(A_sparse)
573+
574+
prob_matrix = LinearProblem(A_sparse, b)
575+
prob_operator = LinearProblem(A_op, b)
576+
577+
# Test KLU with operator
578+
sol_matrix = solve(prob_matrix, KLUFactorization())
579+
sol_operator = solve(prob_operator, KLUFactorization())
580+
@test sol_matrix.u ≈ sol_operator.u
581+
582+
# Test UMFPACK with operator
583+
sol_matrix = solve(prob_matrix, UMFPACKFactorization())
584+
sol_operator = solve(prob_operator, UMFPACKFactorization())
585+
@test sol_matrix.u ≈ sol_operator.u
586+
587+
# Test LU with operator
588+
sol_matrix = solve(prob_matrix, LUFactorization())
589+
sol_operator = solve(prob_operator, LUFactorization())
590+
@test sol_matrix.u ≈ sol_operator.u
591+
592+
# Test WOperator with sparse Jacobian
593+
n_w = 8
594+
M = sparse(I(n_w) * 1.0)
595+
gamma = 1 / 2.0
596+
J = sprand(n_w, n_w, 0.5) + sparse(I(n_w) * 10.0) # Make it diagonally dominant
597+
u = rand(n_w)
598+
b_w = rand(n_w)
599+
600+
W = WOperator{true}(M, gamma, J, u)
601+
W_matrix = convert(AbstractMatrix, W)
602+
603+
prob_woperator = LinearProblem(W, b_w)
604+
prob_wmatrix = LinearProblem(W_matrix, b_w)
605+
606+
# Test KLU with WOperator
607+
sol_woperator = solve(prob_woperator, KLUFactorization())
608+
sol_wmatrix = solve(prob_wmatrix, KLUFactorization())
609+
@test sol_woperator.u ≈ sol_wmatrix.u
610+
611+
# Test UMFPACK with WOperator
612+
sol_woperator = solve(prob_woperator, UMFPACKFactorization())
613+
sol_wmatrix = solve(prob_wmatrix, UMFPACKFactorization())
614+
@test sol_woperator.u ≈ sol_wmatrix.u
615+
616+
# Test LU with WOperator
617+
sol_woperator = solve(prob_woperator, LUFactorization())
618+
sol_wmatrix = solve(prob_wmatrix, LUFactorization())
619+
@test sol_woperator.u ≈ sol_wmatrix.u
620+
end
621+
564622
@testset "Solve Function" begin
565623
A1 = rand(n) |> Diagonal
566624
b1 = rand(n)

0 commit comments

Comments
 (0)