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..ae3fc27f2 100644 --- a/ext/LinearSolveSparseArraysExt.jl +++ b/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) 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)