Merge pull request #847 from ChrisRackauckas/claude/open-linearsol-pr-019f2nZ8STSAH5PaVG3PVNep

Setup use of has_concretization and test with SciMLOperators/WOperator

Author: ChrisRackauckas

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"

ext/LinearSolveSparseArraysExt.jl

@@ -6,6 +6,7 @@ using LinearSolve: LinearSolve, BLASELTYPES, pattern_changed, ArrayInterface,
                    KLUFactorization, LUFactorization, NormalCholeskyFactorization,
                    OperatorAssumptions, LinearVerbosity,
                    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,
@@ -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
+        nothing
+    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
+        nothing
+    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
+        nothing
+    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
+        nothing
+    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
+        nothing
+    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
+        nothing
+    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
+        nothing
+    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
+        nothing
+    end
+end
+
 function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::KLUFactorization; kwargs...)
     A = cache.A
     A = convert(AbstractMatrix, A)

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!

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)