From 8646b3dc65ddf0fc08d37da56e76e6911b9f1cc7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Mon, 19 Aug 2024 19:56:24 +0200 Subject: [PATCH] Introduce AbstractSparseFactorization and AbstractDenseFactorization Both are subtypes of AbstractFactorization. This allows to set default_alias_A=true for sparse factorizations (which just cannot overwrite A anyway). This allows to create the cache parametrized with the original matrix type and we can use `reinit!(cache;A)` or `cache.A=A` if `A` is an AbstractSparseMatrixCSC, like with Krylov solvers. Before, due to https://github.com/SciML/LinearSolve.jl/pull/525, these methods would require a SparseMatrixCSC, even if the Problem was created with a AbstractSparseMatrixCSC. --- src/LinearSolve.jl | 3 +++ src/common.jl | 3 +++ src/factorization.jl | 39 ++++++++++++++++++++------------------- test/basictests.jl | 27 +++++++++++++++++++++++++++ test/resolve.jl | 3 ++- 5 files changed, 55 insertions(+), 20 deletions(-) diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index c2fb6067c..4bb47ba6f 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -62,12 +62,15 @@ using SciMLBase: _unwrap_val abstract type SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm end abstract type AbstractFactorization <: SciMLLinearSolveAlgorithm end +abstract type AbstractSparseFactorization <: AbstractFactorization end +abstract type AbstractDenseFactorization <: AbstractFactorization end abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end # Traits needs_concrete_A(alg::AbstractFactorization) = true +needs_concrete_A(alg::AbstractSparseFactorization) = true needs_concrete_A(alg::AbstractKrylovSubspaceMethod) = false needs_concrete_A(alg::AbstractSolveFunction) = false diff --git a/src/common.jl b/src/common.jl index cf91672ad..f212fe250 100644 --- a/src/common.jl +++ b/src/common.jl @@ -132,6 +132,9 @@ default_alias_b(::Any, ::Any, ::Any) = false default_alias_A(::AbstractKrylovSubspaceMethod, ::Any, ::Any) = true default_alias_b(::AbstractKrylovSubspaceMethod, ::Any, ::Any) = true +default_alias_A(::AbstractSparseFactorization, ::Any, ::Any) = true +default_alias_b(::AbstractSparseFactorization, ::Any, ::Any) = true + DEFAULT_PRECS(A, p) = IdentityOperator(size(A)[1]), IdentityOperator(size(A)[2]) function __init_u0_from_Ab(A, b) diff --git a/src/factorization.jl b/src/factorization.jl index a32f3ede3..0d339f1c4 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -49,7 +49,7 @@ Julia's built in `lu`. Equivalent to calling `lu!(A)` - pivot: The choice of pivoting. Defaults to `LinearAlgebra.RowMaximum()`. The other choice is `LinearAlgebra.NoPivot()`. """ -Base.@kwdef struct LUFactorization{P} <: AbstractFactorization +Base.@kwdef struct LUFactorization{P} <: AbstractDenseFactorization pivot::P = LinearAlgebra.RowMaximum() reuse_symbolic::Bool = true check_pattern::Bool = true # Check factorization re-use @@ -70,7 +70,7 @@ Has low overhead and is good for small matrices. - pivot: The choice of pivoting. Defaults to `LinearAlgebra.RowMaximum()`. The other choice is `LinearAlgebra.NoPivot()`. """ -struct GenericLUFactorization{P} <: AbstractFactorization +struct GenericLUFactorization{P} <: AbstractDenseFactorization pivot::P end @@ -177,7 +177,7 @@ Julia's built in `qr`. Equivalent to calling `qr!(A)`. - On CuMatrix, it will use a CUDA-accelerated QR from CuSolver. - On BandedMatrix and BlockBandedMatrix, it will use a banded QR. """ -struct QRFactorization{P} <: AbstractFactorization +struct QRFactorization{P} <: AbstractDenseFactorization pivot::P blocksize::Int inplace::Bool @@ -260,7 +260,7 @@ Julia's built in `cholesky`. Equivalent to calling `cholesky!(A)`. - shift: the shift argument in CHOLMOD. Only used for sparse matrices. - perm: the perm argument in CHOLMOD. Only used for sparse matrices. """ -struct CholeskyFactorization{P, P2} <: AbstractFactorization +struct CholeskyFactorization{P, P2} <: AbstractDenseFactorization pivot::P tol::Int shift::Float64 @@ -319,7 +319,7 @@ end ## LDLtFactorization -struct LDLtFactorization{T} <: AbstractFactorization +struct LDLtFactorization{T} <: AbstractDenseFactorization shift::Float64 perm::T end @@ -361,7 +361,7 @@ Julia's built in `svd`. Equivalent to `svd!(A)`. which by default is OpenBLAS but will use MKL if the user does `using MKL` in their system. """ -struct SVDFactorization{A} <: AbstractFactorization +struct SVDFactorization{A} <: AbstractDenseFactorization full::Bool alg::A end @@ -410,7 +410,7 @@ Only for Symmetric matrices. - rook: whether to perform rook pivoting. Defaults to false. """ -Base.@kwdef struct BunchKaufmanFactorization <: AbstractFactorization +Base.@kwdef struct BunchKaufmanFactorization <: AbstractDenseFactorization rook::Bool = false end @@ -464,7 +464,7 @@ factorization API. Quoting from Base: - fact_alg: the factorization algorithm to use. Defaults to `LinearAlgebra.factorize`, but can be swapped to choices like `lu`, `qr` """ -struct GenericFactorization{F} <: AbstractFactorization +struct GenericFactorization{F} <: AbstractDenseFactorization fact_alg::F end @@ -781,7 +781,7 @@ patterns with “more structure”. `A` has the same sparsity pattern as the previous `A`. If this algorithm is to be used in a context where that assumption does not hold, set `reuse_symbolic=false`. """ -Base.@kwdef struct UMFPACKFactorization <: AbstractFactorization +Base.@kwdef struct UMFPACKFactorization <: AbstractSparseFactorization reuse_symbolic::Bool = true check_pattern::Bool = true # Check factorization re-use end @@ -860,7 +860,7 @@ A fast sparse LU-factorization which specializes on sparsity patterns with “le `A` has the same sparsity pattern as the previous `A`. If this algorithm is to be used in a context where that assumption does not hold, set `reuse_symbolic=false`. """ -Base.@kwdef struct KLUFactorization <: AbstractFactorization +Base.@kwdef struct KLUFactorization <: AbstractSparseFactorization reuse_symbolic::Bool = true check_pattern::Bool = true end @@ -941,7 +941,7 @@ Only supports sparse matrices. - shift: the shift argument in CHOLMOD. - perm: the perm argument in CHOLMOD """ -Base.@kwdef struct CHOLMODFactorization{T} <: AbstractFactorization +Base.@kwdef struct CHOLMODFactorization{T} <: AbstractSparseFactorization shift::Float64 = 0.0 perm::T = nothing end @@ -993,7 +993,7 @@ implementation, usually outperforming OpenBLAS and MKL for smaller matrices (<500x500), but currently optimized only for Base `Array` with `Float32` or `Float64`. Additional optimization for complex matrices is in the works. """ -struct RFLUFactorization{P, T} <: AbstractFactorization +struct RFLUFactorization{P, T} <: AbstractDenseFactorization RFLUFactorization(::Val{P}, ::Val{T}) where {P, T} = new{P, T}() end @@ -1064,7 +1064,7 @@ be applied to well-conditioned matrices. - pivot: Defaults to RowMaximum(), but can be NoPivot() """ -struct NormalCholeskyFactorization{P} <: AbstractFactorization +struct NormalCholeskyFactorization{P} <: AbstractDenseFactorization pivot::P end @@ -1152,7 +1152,7 @@ be applied to well-conditioned matrices. - rook: whether to perform rook pivoting. Defaults to false. """ -struct NormalBunchKaufmanFactorization <: AbstractFactorization +struct NormalBunchKaufmanFactorization <: AbstractDenseFactorization rook::Bool end @@ -1189,7 +1189,7 @@ end A special implementation only for solving `Diagonal` matrices fast. """ -struct DiagonalFactorization <: AbstractFactorization end +struct DiagonalFactorization <: AbstractDenseFactorization end function init_cacheval(alg::DiagonalFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) @@ -1225,7 +1225,7 @@ end The FastLapackInterface.jl version of the LU factorization. Notably, this version does not allow for choice of pivoting method. """ -struct FastLUFactorization <: AbstractFactorization end +struct FastLUFactorization <: AbstractDenseFactorization end function init_cacheval(::FastLUFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, @@ -1255,7 +1255,7 @@ end The FastLapackInterface.jl version of the QR factorization. """ -struct FastQRFactorization{P} <: AbstractFactorization +struct FastQRFactorization{P} <: AbstractDenseFactorization pivot::P blocksize::Int end @@ -1329,7 +1329,7 @@ dispatch to route around standard BLAS routines in the case e.g. of arbitrary-pr floating point numbers or ForwardDiff.Dual. This e.g. allows for Automatic Differentiation (AD) of a sparse-matrix solve. """ -Base.@kwdef struct SparspakFactorization <: AbstractFactorization +Base.@kwdef struct SparspakFactorization <: AbstractSparseFactorization reuse_symbolic::Bool = true end @@ -1388,7 +1388,8 @@ function SciMLBase.solve!(cache::LinearCache, alg::SparspakFactorization; kwargs SciMLBase.build_linear_solution(alg, y, nothing, cache) end -for alg in InteractiveUtils.subtypes(AbstractFactorization) +for alg in vcat(InteractiveUtils.subtypes(AbstractDenseFactorization), + InteractiveUtils.subtypes(AbstractSparseFactorization)) @eval function init_cacheval(alg::$alg, A::MatrixOperator, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) diff --git a/test/basictests.jl b/test/basictests.jl index 2e69eedc1..6d274a07f 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -547,9 +547,36 @@ end N = 10_000 A = spdiagm(1 => -ones(N - 1), 0 => fill(10.0, N), -1 => -ones(N - 1)) u0 = ones(size(A, 2)) + b = A * u0 B = MySparseMatrixCSC(A) pr = LinearProblem(B, b) + + # test default algorithn @time "solve MySparseMatrixCSC" u=solve(pr) @test norm(u - u0, Inf) < 1.0e-13 + + # test Krylov algorithm with reinit! + pr = LinearProblem(B, b) + solver=KrylovJL_CG() + cache=init(pr,solver,maxiters=1000,reltol=1.0e-10) + u=solve!(cache) + A1 = spdiagm(1 => -ones(N - 1), 0 => fill(100.0, N), -1 => -ones(N - 1)) + b1=A1*u0 + B1= MySparseMatrixCSC(A1) + @test norm(u - u0, Inf) < 1.0e-8 + reinit!(cache; A=B1, b=b1) + u=solve!(cache) + @test norm(u - u0, Inf) < 1.0e-8 + + # test factorization with reinit! + pr = LinearProblem(B, b) + solver=SparspakFactorization() + cache=init(pr,solver) + u=solve!(cache) + @test norm(u - u0, Inf) < 1.0e-8 + reinit!(cache; A=B1, b=b1) + u=solve!(cache) + @test norm(u - u0, Inf) < 1.0e-8 + end diff --git a/test/resolve.jl b/test/resolve.jl index bd2922f54..a567a90b1 100644 --- a/test/resolve.jl +++ b/test/resolve.jl @@ -1,6 +1,7 @@ using LinearSolve, LinearAlgebra, SparseArrays, InteractiveUtils, Test +using LinearSolve: AbstractDenseFactorization, AbstractSparseFactorization -for alg in subtypes(LinearSolve.AbstractFactorization) +for alg in vcat(InteractiveUtils.subtypes(AbstractDenseFactorization),InteractiveUtils.subtypes(AbstractSparseFactorization)) @show alg if !(alg in [ DiagonalFactorization,