From bd03cf18390c263af68677dec06b1188513a837f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 20 Oct 2024 22:30:09 +0200 Subject: [PATCH 01/24] Min version of LinearSolve in test environments to 2.36 --- docs/Project.toml | 1 + test/Project.toml | 1 + 2 files changed, 2 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index 5106403..6179c00 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -14,3 +14,4 @@ Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" [compat] Documenter = "1.0" IterativeSolvers = "0.9" +LinearSolve = "2.36.0" diff --git a/test/Project.toml b/test/Project.toml index 76e3ae8..b780793 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -25,3 +25,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] ExtendableGrids = "1.9" IterativeSolvers = "0.9" +LinearSolve = "2.36.0" From addcd8d249e499ad67b11e9824e182aa3c2d80ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 20 Oct 2024 22:30:45 +0200 Subject: [PATCH 02/24] rework test_linearsolve.jl --- test/test_linearsolve.jl | 70 +++++++++++++++++++++++----------------- test/test_pardiso.jl | 3 -- 2 files changed, 41 insertions(+), 32 deletions(-) diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index db161cd..14df0a7 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -7,6 +7,7 @@ using LinearAlgebra using LinearSolve using ForwardDiff using MultiFloats +import Pardiso f64(x::ForwardDiff.Dual{T}) where {T} = Float64(ForwardDiff.value(x)) f64(x::Number) = Float64(x) @@ -17,59 +18,70 @@ function test_ls1(T, k, l, m; linsolver = SparspakFactorization()) b = A*ones(k * l * m) x0 = A \ b p = LinearProblem(T.(A), T.(b)) - x1 = solve(p, linsolver) + x1 = solve(p, linsolver, abstol=1.0e-12) + @show norm(x0-x1) x0 ≈ x1 end -@test test_ls1(Float64, 10, 10, 10, linsolver = KLUFactorization()) -@test test_ls1(Float64, 25, 40, 1, linsolver = KLUFactorization()) -@test test_ls1(Float64, 100, 1, 1, linsolver = KLUFactorization()) - -for T in [Float32, Float64, Float64x1, Float64x2, Dual64] - println("$T:") - @test test_ls1(T, 10, 10, 10, linsolver = SparspakFactorization()) - @test test_ls1(T, 25, 40, 1, linsolver = SparspakFactorization()) - @test test_ls1(T, 100, 1, 1, linsolver = SparspakFactorization()) -end - function test_ls2(T, k, l, m; linsolver = SparspakFactorization()) A = fdrand(T, k, l, m; rand = () -> 1, matrixtype = ExtendableSparseMatrix) b = T.(rand(k * l * m)) p = LinearProblem(A, b) x0 = solve(p, linsolver) cache = x0.cache - x0 = A \ b - for i = 4:(k * l * m - 3) - A[i, i + 3] -= 1.0e-4 - A[i - 3, i] -= 1.0e-4 + x0=copy(x0) + nonzeros(A).-=1.0e-4 + for i = 1:k*l*m + A[i, i] += 1.0e-4 end - LinearSolve.set_A(cache, A) - x1 = solve(p, linsolver) - x1 = A \ b + reinit!(cache; A, reuse_precs=true) + x1 = solve!(cache, linsolver) all(x0 .< x1) end -@test test_ls1(Float64, 10, 10, 10, linsolver = KLUFactorization()) -@test test_ls1(Float64, 25, 40, 1, linsolver = KLUFactorization()) -@test test_ls1(Float64, 100, 1, 1, linsolver = KLUFactorization()) + for T in [Float32, Float64, Float64x1, Float64x2, Dual64] println("$T:") @test test_ls1(T, 10, 10, 10, linsolver = SparspakFactorization()) @test test_ls1(T, 25, 40, 1, linsolver = SparspakFactorization()) @test test_ls1(T, 100, 1, 1, linsolver = SparspakFactorization()) -end -@test test_ls2(Float64, 10, 10, 10, linsolver = KLUFactorization()) -@test test_ls2(Float64, 25, 40, 1, linsolver = KLUFactorization()) -@test test_ls2(Float64, 100, 1, 1, linsolver = KLUFactorization()) - -for T in [Float32, Float64, Float64x1, Float64x2, Dual64] - println("$T:") @test test_ls2(T, 10, 10, 10, linsolver = SparspakFactorization()) @test test_ls2(T, 25, 40, 1, linsolver = SparspakFactorization()) @test test_ls2(T, 100, 1, 1, linsolver = SparspakFactorization()) + end + +for factorization in [UMFPACKFactorization(), + KLUFactorization(reuse_symbolic=false), + MKLPardisoFactorize()] + println("$factorization:") + @test test_ls1(Float64, 10, 10, 10, linsolver = factorization) + @test test_ls1(Float64, 25, 40, 1, linsolver = factorization) + @test test_ls1(Float64, 100, 1, 1, linsolver = factorization) + + @test test_ls2(Float64, 10, 10, 10, linsolver = factorization) + @test test_ls2(Float64, 25, 40, 1, linsolver = factorization) + @test test_ls2(Float64, 100, 1, 1, linsolver = factorization) +end + + + +for iteration in [KrylovJL_GMRES(precs= (A,p) -> (JacobiPreconditioner(A),I))] + println("$iteration:") + @test test_ls1(Float64, 10, 10, 10, linsolver = iteration) + @test test_ls1(Float64, 25, 40, 1, linsolver = iteration) + @test test_ls1(Float64, 100, 1, 1, linsolver = iteration) + + @test test_ls2(Float64, 10, 10, 10, linsolver = iteration) + @test test_ls2(Float64, 25, 40, 1, linsolver = iteration) + @test test_ls2(Float64, 100, 1, 1, linsolver = iteration) +end + + + + end diff --git a/test/test_pardiso.jl b/test/test_pardiso.jl index 24f6081..1bc07b7 100644 --- a/test/test_pardiso.jl +++ b/test/test_pardiso.jl @@ -14,7 +14,4 @@ include("test_lu.jl") @test test_lu2(Float64, 10, 10, 10, lufac = PardisoLU()) @test test_lu2(Float64, 25, 40, 1, lufac = PardisoLU()) @test test_lu2(Float64, 1000, 1, 1, lufac = PardisoLU()) -# @test test_lu2(10,10,10,lufac=PardisoLU(mtype=2)) -# @test test_lu2(25,40,1,lufac=PardisoLU(mtype=2)) -# @test test_lu2(1000,1,1,lufac=PardisoLU(mtype=2)) end From 2c034026fd2c8c78ea11f157e961c02d02bbdc26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 2 Nov 2024 21:34:12 +0100 Subject: [PATCH 03/24] start adaptation to new `precs` API of LinearSolve --- Project.toml | 2 +- ext/ExtendableSparseAMGCLWrapExt.jl | 4 +++ ext/ExtendableSparseIncompleteLUExt.jl | 6 ++++ src/ExtendableSparse.jl | 5 ++- src/factorizations/blockpreconditioner.jl | 16 ++++++++- src/precs.jl | 7 ++++ test/Project.toml | 1 + test/test_block.jl | 3 +- test/test_linearsolve.jl | 42 +++++++++++++++++++++-- 9 files changed, 80 insertions(+), 6 deletions(-) create mode 100644 src/precs.jl diff --git a/Project.toml b/Project.toml index 2add3ca..020312f 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ ExtendableSparseIncompleteLUExt = "IncompleteLU" ExtendableSparsePardisoExt = "Pardiso" [compat] -AMGCLWrap = "0.4,1" +AMGCLWrap = "2" AlgebraicMultigrid = "0.4,0.5,0.6" DocStringExtensions = "0.8, 0.9" ILUZero = "0.2" diff --git a/ext/ExtendableSparseAMGCLWrapExt.jl b/ext/ExtendableSparseAMGCLWrapExt.jl index 3389ce8..87b020b 100644 --- a/ext/ExtendableSparseAMGCLWrapExt.jl +++ b/ext/ExtendableSparseAMGCLWrapExt.jl @@ -10,6 +10,8 @@ mutable struct AMGCL_AMGPreconditioner <: AbstractPreconditioner factorization::AMGCLWrap.AMGPrecon kwargs function ExtendableSparse.AMGCL_AMGPreconditioner(; kwargs...) + Base.depwarn("AMGCL_AMGPreconditioner() is deprecated. Use LinearSolve with `precs=AMGCLWrap.AMGPreconBuilder()` instead.", + :AMGCL_AMGPreconditioner) precon = new() precon.kwargs = kwargs precon @@ -35,6 +37,8 @@ mutable struct AMGCL_RLXPreconditioner <: AbstractPreconditioner factorization::AMGCLWrap.RLXPrecon kwargs function ExtendableSparse.AMGCL_RLXPreconditioner(; kwargs...) + Base.depwarn("AMGCL_RLXPreconditioner() is deprecated. Use LinearSolve with `precs=AMGCLWrap.RLXPreconBuilder()` instead.", + :AMGCL_RLXPreconditioner) precon = new() precon.kwargs = kwargs precon diff --git a/ext/ExtendableSparseIncompleteLUExt.jl b/ext/ExtendableSparseIncompleteLUExt.jl index dc6ff25..08b70ac 100644 --- a/ext/ExtendableSparseIncompleteLUExt.jl +++ b/ext/ExtendableSparseIncompleteLUExt.jl @@ -1,9 +1,14 @@ module ExtendableSparseIncompleteLUExt using ExtendableSparse using IncompleteLU +using LinearAlgebra: I +using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC, getcolptr, rowvals, nonzeros + import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update! +import ExtendableSparse: IncompleteLUPrecs + mutable struct ILUTPreconditioner <: AbstractPreconditioner A::ExtendableSparseMatrix factorization::IncompleteLU.ILUFactorization @@ -25,5 +30,6 @@ function update!(precon::ILUTPreconditioner) @inbounds flush!(A) precon.factorization = IncompleteLU.ilu(A.cscmatrix; τ = precon.droptol) end + end diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 9beee11..373a299 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -3,7 +3,7 @@ module ExtendableSparse using DocStringExtensions: DocStringExtensions, SIGNATURES, TYPEDEF,TYPEDFIELDS using ILUZero: ILUZero, ldiv!, nnz using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, Symmetric, Tridiagonal, - cholesky, cholesky!, convert, lu!, mul!, norm, transpose + cholesky, cholesky!, convert, lu!, mul!, norm, transpose, I using SparseArrays: SparseArrays, AbstractSparseMatrix, SparseMatrixCSC, dropzeros!, findnz, nzrange, sparse, spzeros using Sparspak: Sparspak, sparspaklu, sparspaklu! @@ -52,6 +52,9 @@ include("factorizations/factorizations.jl") include("factorizations/simple_iteration.jl") export simple, simple! +include("precs.jl") +export SparspakPrecs, UMFPACKPrecs, EquationBlockPrecs + include("matrix/sprand.jl") export sprand!, sprand_sdd!, fdrand, fdrand!, fdrand_coo, solverbenchmark diff --git a/src/factorizations/blockpreconditioner.jl b/src/factorizations/blockpreconditioner.jl index af13986..8017966 100644 --- a/src/factorizations/blockpreconditioner.jl +++ b/src/factorizations/blockpreconditioner.jl @@ -30,7 +30,7 @@ function BlockPreconditioner end Factorizations on matrix partitions within a block preconditioner may or may not work with array views. E.g. the umfpack factorization cannot work with views, while ILUZeroPreconditioner can. -Implementing a method for `allow_views` returning `false` resp. `true` allows to dispatch to the proper case. + Implementing a method for `allow_views` returning `false` resp. `true` allows to dispatch to the proper case. """ allow_views(::Any)=false @@ -98,3 +98,17 @@ function LinearAlgebra.ldiv!(u,p::BlockPreconditioner,v) end Base.eltype(p::BlockPreconditioner)=eltype(p.facts[1]) + + + +Base.@kwdef struct EquationBlockPrecs + precs=UMFPACKPrecs() + partitioning= [A -> 1:size(A,1)] +end + +function (blockprecs::EquationBlockPrecs)(A,p) + (;precs, partitioning)=blockprecs + factorization= A->precs(A,p)[1] + bp=BlockPreconditioner(A;partitioning=partitioning(A), factorization) + (bp,I) +end diff --git a/src/precs.jl b/src/precs.jl new file mode 100644 index 0000000..3ff6dd6 --- /dev/null +++ b/src/precs.jl @@ -0,0 +1,7 @@ +struct UMFPACKPrecs end +(::UMFPACKPrecs)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) + +struct SparspakPrecs end +(::SparspakPrecs)(A::AbstractSparseMatrixCSC,p)=(sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) + + diff --git a/test/Project.toml b/test/Project.toml index b780793..e6a6a68 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -15,6 +15,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5" +Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" diff --git a/test/test_block.jl b/test/test_block.jl index 1a055a1..22cabc1 100644 --- a/test/test_block.jl +++ b/test/test_block.jl @@ -5,7 +5,7 @@ using ILUZero, AlgebraicMultigrid using IterativeSolvers using LinearAlgebra using Sparspak - +using AMGCLWrap ExtendableSparse.allow_views(::typeof(ilu0))=true @@ -48,4 +48,5 @@ function main(;n=100) end main(n=100) + end diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index 14df0a7..4d3688e 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -7,6 +7,8 @@ using LinearAlgebra using LinearSolve using ForwardDiff using MultiFloats + +using AMGCLWrap import Pardiso f64(x::ForwardDiff.Dual{T}) where {T} = Float64(ForwardDiff.value(x)) @@ -19,7 +21,6 @@ function test_ls1(T, k, l, m; linsolver = SparspakFactorization()) x0 = A \ b p = LinearProblem(T.(A), T.(b)) x1 = solve(p, linsolver, abstol=1.0e-12) - @show norm(x0-x1) x0 ≈ x1 end @@ -70,7 +71,11 @@ end -for iteration in [KrylovJL_GMRES(precs= (A,p) -> (JacobiPreconditioner(A),I))] +for iteration in [ + KrylovJL_GMRES(precs=AMGCLWrap.AMGPreconBuilder()), + KrylovJL_GMRES(precs=AMGCLWrap.RLXPreconBuilder()) + + ] println("$iteration:") @test test_ls1(Float64, 10, 10, 10, linsolver = iteration) @test test_ls1(Float64, 25, 40, 1, linsolver = iteration) @@ -83,5 +88,38 @@ end +function mainprecs(;n=100) + A=fdrand(n,n) + partitioning=A->[1:2:size(A,1), 2:2:size(A,1)] + sol0=ones(n^2) + b=A*ones(n^2); + + precs=EquationBlockPrecs(;precs=UMFPACKPrecs(), partitioning) + @info precs + p=LinearProblem(A,b) + sol=solve(p, KrylovJL_CG(;precs)) + @test sol≈sol0 + + precs=EquationBlockPrecs(;precs=SparspakPrecs(), partitioning) + @info precs + p=LinearProblem(A,b) + sol=solve(p, KrylovJL_CG(;precs)) + @test sol≈sol0 + + precs=EquationBlockPrecs(;precs=AMGCLWrap.AMGPreconBuilder(), partitioning) + @info precs + p=LinearProblem(A,b) + sol=solve(p, KrylovJL_CG(;precs)) + @test sol≈sol0 + + precs=EquationBlockPrecs(;precs=AMGCLWrap.RLXPreconBuilder(), partitioning) + @info precs + p=LinearProblem(A,b) + sol=solve(p, KrylovJL_CG(;precs)) + @test sol≈sol0 + +end +mainprecs() + end From b62d551028f22099ced15cdcd86d913e6e10644d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 00:30:35 +0100 Subject: [PATCH 04/24] Add compat.jl defining @public --- src/compat.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 src/compat.jl diff --git a/src/compat.jl b/src/compat.jl new file mode 100644 index 0000000..570f712 --- /dev/null +++ b/src/compat.jl @@ -0,0 +1,18 @@ +# Copied from MIT licensed +# https://github.com/JuliaDiff/DifferentiationInterface.jl/blob/main/DifferentiationInterface/src/compat.jl +# +macro public(ex) + if VERSION >= v"1.11.0-DEV.469" + args = if ex isa Symbol + (ex,) + elseif Base.isexpr(ex, :tuple) + ex.args + else + error("something informative") + end + esc(Expr(:public, args...)) + else + nothing + end +end + From b9b15b98e0d3893b2e9743822b6835ffd566c4da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 00:31:12 +0100 Subject: [PATCH 05/24] comprehensive precs API --- ext/ExtendableSparseAMGCLWrapExt.jl | 2 + ext/ExtendableSparseAlgebraicMultigridExt.jl | 16 ++- ext/ExtendableSparseIncompleteLUExt.jl | 6 +- src/ExtendableSparse.jl | 13 +- src/factorizations/blockpreconditioner.jl | 9 +- src/preconbuilders.jl | 42 +++++++ src/precs.jl | 7 -- test/test_linearsolve.jl | 124 +++++++++---------- 8 files changed, 137 insertions(+), 82 deletions(-) create mode 100644 src/preconbuilders.jl delete mode 100644 src/precs.jl diff --git a/ext/ExtendableSparseAMGCLWrapExt.jl b/ext/ExtendableSparseAMGCLWrapExt.jl index 87b020b..f078a0e 100644 --- a/ext/ExtendableSparseAMGCLWrapExt.jl +++ b/ext/ExtendableSparseAMGCLWrapExt.jl @@ -1,3 +1,5 @@ +# The whole extension is deprecated +# TODO remove in v2.0 module ExtendableSparseAMGCLWrapExt using ExtendableSparse using AMGCLWrap diff --git a/ext/ExtendableSparseAlgebraicMultigridExt.jl b/ext/ExtendableSparseAlgebraicMultigridExt.jl index a21b935..4588d84 100644 --- a/ext/ExtendableSparseAlgebraicMultigridExt.jl +++ b/ext/ExtendableSparseAlgebraicMultigridExt.jl @@ -1,6 +1,20 @@ module ExtendableSparseAlgebraicMultigridExt using ExtendableSparse -using AlgebraicMultigrid +using AlgebraicMultigrid: AlgebraicMultigrid, ruge_stuben, smoothed_aggregation, aspreconditioner +using SparseArrays: SparseMatrixCSC, AbstractSparseMatrixCSC +using LinearAlgebra: I + + +import ExtendableSparse: SmoothedAggregationAMGBuilder +import ExtendableSparse: RugeStubenAMGBuilder + +(b::SmoothedAggregationAMGBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I) +(b::RugeStubenAMGBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(ruge_stuben(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I) + + +#### +# Deprecated from here on +# TODO remove in v2.0 import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update! diff --git a/ext/ExtendableSparseIncompleteLUExt.jl b/ext/ExtendableSparseIncompleteLUExt.jl index 08b70ac..f777af9 100644 --- a/ext/ExtendableSparseIncompleteLUExt.jl +++ b/ext/ExtendableSparseIncompleteLUExt.jl @@ -4,11 +4,15 @@ using IncompleteLU using LinearAlgebra: I using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC, getcolptr, rowvals, nonzeros +import ExtendableSparse: ILUTBuilder + +(b::ILUTBuilder)(A::AbstractSparseMatrixCSC,p)=(IncompleteLU.ilu(SparseMatrixCSC(A); τ = b.droptol),I) + import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update! -import ExtendableSparse: IncompleteLUPrecs +# Deprecated from here mutable struct ILUTPreconditioner <: AbstractPreconditioner A::ExtendableSparseMatrix factorization::IncompleteLU.ILUFactorization diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 373a299..dda586e 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -21,6 +21,7 @@ if USE_GPL_LIBS end +include("compat.jl") # @public include("matrix/sparsematrixcsc.jl") include("matrix/abstractsparsematrixextension.jl") @@ -52,8 +53,11 @@ include("factorizations/factorizations.jl") include("factorizations/simple_iteration.jl") export simple, simple! -include("precs.jl") -export SparspakPrecs, UMFPACKPrecs, EquationBlockPrecs +include("preconbuilders.jl") +export SparspakPreconBuilder, UMFPACKPreconBuilder, EquationBlockPreconBuilder, JacobiPreconBuilder + +@public ILUZeroBuilder, SmoothedAggregationAMGBuilder, RugeStubenAMGBuilder + include("matrix/sprand.jl") export sprand!, sprand_sdd!, fdrand, fdrand!, fdrand_coo, solverbenchmark @@ -61,7 +65,9 @@ export sprand!, sprand_sdd!, fdrand, fdrand!, fdrand_coo, solverbenchmark export rawupdateindex!, updateindex! - +##### +# All of the following is deprecated in favor of the precs bases +# API export JacobiPreconditioner, ILU0Preconditioner, @@ -211,4 +217,5 @@ Pardiso.set_iparm!(plu.ps,5,13.0) function MKLPardisoLU end export MKLPardisoLU + end # module diff --git a/src/factorizations/blockpreconditioner.jl b/src/factorizations/blockpreconditioner.jl index 8017966..1d3d036 100644 --- a/src/factorizations/blockpreconditioner.jl +++ b/src/factorizations/blockpreconditioner.jl @@ -101,14 +101,19 @@ Base.eltype(p::BlockPreconditioner)=eltype(p.facts[1]) -Base.@kwdef struct EquationBlockPrecs +Base.@kwdef struct EquationBlockPreconBuilder precs=UMFPACKPrecs() partitioning= [A -> 1:size(A,1)] end -function (blockprecs::EquationBlockPrecs)(A,p) +function (blockprecs::EquationBlockPreconBuilder)(A,p) (;precs, partitioning)=blockprecs factorization= A->precs(A,p)[1] bp=BlockPreconditioner(A;partitioning=partitioning(A), factorization) (bp,I) end + +""" + Allow array for precs => different precoms +EquationBlockPreconBuilder +""" diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl new file mode 100644 index 0000000..033fb5c --- /dev/null +++ b/src/preconbuilders.jl @@ -0,0 +1,42 @@ + +struct UMFPACKPreconBuilder end +(::UMFPACKPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) + +struct SparspakPreconBuilder end +(::SparspakPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) + +struct JacobiPreconBuilder end +(::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(JacobiPreconditioner(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) + + +struct ILUZeroBuilder end +(::ILUZeroBuilder)(A,p)=(ilu0(A),I) + + + +Base.@kwdef struct ILUTBuilder + droptol::Float64=0.1 +end +(::ILUTBuilder)(A,p)= error("import IncompleteLU.jl in order to use ILUTBuilder") + +struct SmoothedAggregationAMGBuilder{Tk} + blocksize::Int + kwargs::Tk +end + +function SmoothedAggregationAMGBuilder(;blocksize=1, kwargs...) + SmoothedAggregationAMGBuilder(blocksize,kwargs) +end + +(::SmoothedAggregationAMGBuilder)(A,p)= error("import AlgebraicMultigrid in order to use SmoothedAggregationAMGBuilder") + +struct RugeStubenAMGBuilder{Tk} + blocksize::Int + kwargs::Tk +end + +function RugeStubenAMGBuilder(;blocksize=1, kwargs...) + SmoothedAggregationAMGBuilder(blocksize,kwargs) +end + +(::RugeStubenAMGBuilder)(A,p)= error("import AlgebraicMultigrid in order to use RugeStubenAMGBuilder") diff --git a/src/precs.jl b/src/precs.jl deleted file mode 100644 index 3ff6dd6..0000000 --- a/src/precs.jl +++ /dev/null @@ -1,7 +0,0 @@ -struct UMFPACKPrecs end -(::UMFPACKPrecs)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) - -struct SparspakPrecs end -(::SparspakPrecs)(A::AbstractSparseMatrixCSC,p)=(sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) - - diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index 4d3688e..f133716 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -9,6 +9,7 @@ using ForwardDiff using MultiFloats using AMGCLWrap +using ILUZero, IncompleteLU, AlgebraicMultigrid import Pardiso f64(x::ForwardDiff.Dual{T}) where {T} = Float64(ForwardDiff.value(x)) @@ -42,84 +43,71 @@ function test_ls2(T, k, l, m; linsolver = SparspakFactorization()) end - -for T in [Float32, Float64, Float64x1, Float64x2, Dual64] - println("$T:") - @test test_ls1(T, 10, 10, 10, linsolver = SparspakFactorization()) - @test test_ls1(T, 25, 40, 1, linsolver = SparspakFactorization()) - @test test_ls1(T, 100, 1, 1, linsolver = SparspakFactorization()) - - @test test_ls2(T, 10, 10, 10, linsolver = SparspakFactorization()) - @test test_ls2(T, 25, 40, 1, linsolver = SparspakFactorization()) - @test test_ls2(T, 100, 1, 1, linsolver = SparspakFactorization()) +@testset "Sparspak" begin + for T in [Float32, Float64, Float64x1, Float64x2, Dual64] + @test test_ls1(T, 10, 10, 10, linsolver = SparspakFactorization()) + @test test_ls1(T, 25, 40, 1, linsolver = SparspakFactorization()) + @test test_ls1(T, 100, 1, 1, linsolver = SparspakFactorization()) + + @test test_ls2(T, 10, 10, 10, linsolver = SparspakFactorization()) + @test test_ls2(T, 25, 40, 1, linsolver = SparspakFactorization()) + @test test_ls2(T, 100, 1, 1, linsolver = SparspakFactorization()) + end end - - -for factorization in [UMFPACKFactorization(), - KLUFactorization(reuse_symbolic=false), - MKLPardisoFactorize()] - println("$factorization:") - @test test_ls1(Float64, 10, 10, 10, linsolver = factorization) - @test test_ls1(Float64, 25, 40, 1, linsolver = factorization) - @test test_ls1(Float64, 100, 1, 1, linsolver = factorization) - - @test test_ls2(Float64, 10, 10, 10, linsolver = factorization) - @test test_ls2(Float64, 25, 40, 1, linsolver = factorization) - @test test_ls2(Float64, 100, 1, 1, linsolver = factorization) +@testset "Factorizations" begin + + for factorization in [UMFPACKFactorization(), + KLUFactorization(reuse_symbolic=false), + MKLPardisoFactorize()] + + @test test_ls1(Float64, 10, 10, 10, linsolver = factorization) + @test test_ls1(Float64, 25, 40, 1, linsolver = factorization) + @test test_ls1(Float64, 100, 1, 1, linsolver = factorization) + + @test test_ls2(Float64, 10, 10, 10, linsolver = factorization) + @test test_ls2(Float64, 25, 40, 1, linsolver = factorization) + @test test_ls2(Float64, 100, 1, 1, linsolver = factorization) + end end - - -for iteration in [ - KrylovJL_GMRES(precs=AMGCLWrap.AMGPreconBuilder()), - KrylovJL_GMRES(precs=AMGCLWrap.RLXPreconBuilder()) - - ] - println("$iteration:") - @test test_ls1(Float64, 10, 10, 10, linsolver = iteration) - @test test_ls1(Float64, 25, 40, 1, linsolver = iteration) - @test test_ls1(Float64, 100, 1, 1, linsolver = iteration) - - @test test_ls2(Float64, 10, 10, 10, linsolver = iteration) - @test test_ls2(Float64, 25, 40, 1, linsolver = iteration) - @test test_ls2(Float64, 100, 1, 1, linsolver = iteration) +allprecs=[ + AMGCLWrap.AMGPreconBuilder(), + AMGCLWrap.AMGPreconBuilder(), + AMGCLWrap.RLXPreconBuilder(), + ExtendableSparse.ILUZeroBuilder(), + ExtendableSparse.ILUTBuilder(), + ExtendableSparse.SmoothedAggregationAMGBuilder(), + ExtendableSparse.RugeStubenAMGBuilder() +] + +@testset "iterations" begin + for precs in allprecs + iteration=KrylovJL_GMRES(precs;) + + @test test_ls1(Float64, 10, 10, 10, linsolver = iteration) + @test test_ls1(Float64, 25, 40, 1, linsolver = iteration) + @test test_ls1(Float64, 100, 1, 1, linsolver = iteration) + + @test test_ls2(Float64, 10, 10, 10, linsolver = iteration) + @test test_ls2(Float64, 25, 40, 1, linsolver = iteration) + @test test_ls2(Float64, 100, 1, 1, linsolver = iteration) + end end - - -function mainprecs(;n=100) +@testset "equationblock" begin + n=100 A=fdrand(n,n) partitioning=A->[1:2:size(A,1), 2:2:size(A,1)] sol0=ones(n^2) b=A*ones(n^2); - - precs=EquationBlockPrecs(;precs=UMFPACKPrecs(), partitioning) - @info precs - p=LinearProblem(A,b) - sol=solve(p, KrylovJL_CG(;precs)) - @test sol≈sol0 - - precs=EquationBlockPrecs(;precs=SparspakPrecs(), partitioning) - @info precs - p=LinearProblem(A,b) - sol=solve(p, KrylovJL_CG(;precs)) - @test sol≈sol0 - - precs=EquationBlockPrecs(;precs=AMGCLWrap.AMGPreconBuilder(), partitioning) - @info precs - p=LinearProblem(A,b) - sol=solve(p, KrylovJL_CG(;precs)) - @test sol≈sol0 - - precs=EquationBlockPrecs(;precs=AMGCLWrap.RLXPreconBuilder(), partitioning) - @info precs - p=LinearProblem(A,b) - sol=solve(p, KrylovJL_CG(;precs)) - @test sol≈sol0 - + + for precs in allprecs + iteration=KrylovJL_CG(precs=EquationBlockPreconBuilder(;precs, partitioning)) + p=LinearProblem(A,b) + sol=solve(p, KrylovJL_CG(;precs)) + @test sol≈sol0 + end end -mainprecs() - end From eddbb04630eb43adfd41d1a290db50a0ad6b03fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 00:37:47 +0100 Subject: [PATCH 06/24] update changelog --- CHANGELOG.md | 11 +++++++++++ Project.toml | 2 +- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fc41761..43fab43 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # Changelog +## [2.0.0] - Planned + +### Breaking +- remove solver + precon API which is not based on precs or directly overloading `\`. + Fully rely on LinearSolve (besides `\`) +- Move AMGBuilder, ILUZeroBuilder to the correspondig packages (depending on the PRs) +- remove "old" SparseMatrixLNK (need to benchmark before) + +## [1.6.0] - WIP +- Support precs API of LinearSolve.jl + ## [1.5.3] - 2024-10-07 - Moved repo to WIAS-PDELib organization diff --git a/Project.toml b/Project.toml index ae35070..20345a0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableSparse" uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" authors = ["Juergen Fuhrmann "] -version = "1.5.3" +version = "1.6.0-dev" [deps] AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288" From be7a688af2d9b53c70d22e91f7e93d393f669e70 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 01:00:20 +0100 Subject: [PATCH 07/24] fix nagging (why???) on stale explicit import of I from LinearAlgebra --- src/ExtendableSparse.jl | 2 +- src/preconbuilders.jl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index dda586e..944178f 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -3,7 +3,7 @@ module ExtendableSparse using DocStringExtensions: DocStringExtensions, SIGNATURES, TYPEDEF,TYPEDFIELDS using ILUZero: ILUZero, ldiv!, nnz using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, Symmetric, Tridiagonal, - cholesky, cholesky!, convert, lu!, mul!, norm, transpose, I + cholesky, cholesky!, convert, lu!, mul!, norm, transpose using SparseArrays: SparseArrays, AbstractSparseMatrix, SparseMatrixCSC, dropzeros!, findnz, nzrange, sparse, spzeros using Sparspak: Sparspak, sparspaklu, sparspaklu! diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 033fb5c..e22e898 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -1,16 +1,16 @@ struct UMFPACKPreconBuilder end -(::UMFPACKPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) +(::UMFPACKPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) struct SparspakPreconBuilder end -(::SparspakPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) +(::SparspakPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) struct JacobiPreconBuilder end -(::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(JacobiPreconditioner(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),I) +(::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(JacobiPreconditioner(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) struct ILUZeroBuilder end -(::ILUZeroBuilder)(A,p)=(ilu0(A),I) +(::ILUZeroBuilder)(A,p)=(ilu0(A),LinearAlgebra.I) From 116b31b876537bd2e00ebc73864561b48a11ac75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 01:03:49 +0100 Subject: [PATCH 08/24] test on lts min, update actions --- .github/workflows/ci.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c99162a..fb88e6a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,9 +15,9 @@ jobs: fail-fast: false matrix: version: - - '1.9' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. + - 'lts' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. -# - 'nightly' + - 'nightly' os: - ubuntu-latest - windows-latest @@ -34,7 +34,7 @@ jobs: arch: x64 steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} @@ -59,7 +59,7 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: '1' - run: | From 52d2f6581df2c48c22dd6f331901f1b69c6f38e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 01:05:16 +0100 Subject: [PATCH 09/24] remove stale factorizations/linearsolve.jl --- src/factorizations/linearsolve.jl | 17 ----------------- 1 file changed, 17 deletions(-) delete mode 100644 src/factorizations/linearsolve.jl diff --git a/src/factorizations/linearsolve.jl b/src/factorizations/linearsolve.jl deleted file mode 100644 index 81bfe57..0000000 --- a/src/factorizations/linearsolve.jl +++ /dev/null @@ -1,17 +0,0 @@ -""" - $(SIGNATURES) - -Create LinearProblem from ExtendableSparseMatrix. -""" -function LinearSolve.LinearProblem(A::ExtendableSparseMatrix, b::AbstractArray) - LinearSolve.LinearProblem(SparseMatrixCSC(A), b) -end - -""" - $(SIGNATURES) - -Update LinearSolve cache. -""" -function LinearSolve.set_A(cache::LinearSolve.LinearCache, A::ExtendableSparseMatrix) - LinearSolve.set_A(cache, SparseMatrixCSC(A)) -end From 557bc69f7e4f78cfabc5d8784e2aa90c5effb96a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 01:19:34 +0100 Subject: [PATCH 10/24] tweek test tolerance for windows --- test/test_linearsolve.jl | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index f133716..3b8a0dd 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -55,11 +55,17 @@ end end end + +factorizations=[UMFPACKFactorization(), + KLUFactorization(reuse_symbolic=false)] + +if !Sys.isapple() + push!(factorizations,MKLPardisoFactorize()) +end + @testset "Factorizations" begin - for factorization in [UMFPACKFactorization(), - KLUFactorization(reuse_symbolic=false), - MKLPardisoFactorize()] + for factorization in factorizations @test test_ls1(Float64, 10, 10, 10, linsolver = factorization) @test test_ls1(Float64, 25, 40, 1, linsolver = factorization) @@ -105,8 +111,8 @@ end for precs in allprecs iteration=KrylovJL_CG(precs=EquationBlockPreconBuilder(;precs, partitioning)) p=LinearProblem(A,b) - sol=solve(p, KrylovJL_CG(;precs)) - @test sol≈sol0 + sol=solve(p, KrylovJL_CG(;precs), abstol=1.0e-12) + @test isapprox(sol, sol0, atol=1e-6) end end From 72fd4a5a63a5bb06cf1f47145fce679d431492a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 22:01:12 +0100 Subject: [PATCH 11/24] fix public names All callable objects describing preconditioners are now PreconBuilder. --- ext/ExtendableSparseAlgebraicMultigridExt.jl | 8 +-- ext/ExtendableSparseIncompleteLUExt.jl | 4 +- src/ExtendableSparse.jl | 4 +- src/factorizations/blockpreconditioner.jl | 20 ++++-- src/preconbuilders.jl | 71 ++++++++++++++++---- test/test_linearsolve.jl | 10 +-- 6 files changed, 87 insertions(+), 30 deletions(-) diff --git a/ext/ExtendableSparseAlgebraicMultigridExt.jl b/ext/ExtendableSparseAlgebraicMultigridExt.jl index 4588d84..9d69470 100644 --- a/ext/ExtendableSparseAlgebraicMultigridExt.jl +++ b/ext/ExtendableSparseAlgebraicMultigridExt.jl @@ -5,11 +5,11 @@ using SparseArrays: SparseMatrixCSC, AbstractSparseMatrixCSC using LinearAlgebra: I -import ExtendableSparse: SmoothedAggregationAMGBuilder -import ExtendableSparse: RugeStubenAMGBuilder +import ExtendableSparse: SmoothedAggregationPreconBuilder +import ExtendableSparse: RugeStubenPreconBuilder -(b::SmoothedAggregationAMGBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I) -(b::RugeStubenAMGBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(ruge_stuben(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I) +(b::SmoothedAggregationPreconBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I) +(b::RugeStubenPreconBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(ruge_stuben(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I) #### diff --git a/ext/ExtendableSparseIncompleteLUExt.jl b/ext/ExtendableSparseIncompleteLUExt.jl index f777af9..6074574 100644 --- a/ext/ExtendableSparseIncompleteLUExt.jl +++ b/ext/ExtendableSparseIncompleteLUExt.jl @@ -4,9 +4,9 @@ using IncompleteLU using LinearAlgebra: I using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC, getcolptr, rowvals, nonzeros -import ExtendableSparse: ILUTBuilder +import ExtendableSparse: ILUTPreconBuilder -(b::ILUTBuilder)(A::AbstractSparseMatrixCSC,p)=(IncompleteLU.ilu(SparseMatrixCSC(A); τ = b.droptol),I) +(b::ILUTPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(IncompleteLU.ilu(SparseMatrixCSC(A); τ = b.droptol),I) import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update! diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 944178f..cefa54e 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -54,9 +54,9 @@ include("factorizations/simple_iteration.jl") export simple, simple! include("preconbuilders.jl") -export SparspakPreconBuilder, UMFPACKPreconBuilder, EquationBlockPreconBuilder, JacobiPreconBuilder +export SparspakPreconBuilder, UMFPACKPreconBuilder, BlockPreconBuilder, JacobiPreconBuilder -@public ILUZeroBuilder, SmoothedAggregationAMGBuilder, RugeStubenAMGBuilder +@public ILUZeroPreconBuilder, ILUTPreconBuilder, SmoothedAggregationPreconBuilder, RugeStubenPreconBuilder include("matrix/sprand.jl") diff --git a/src/factorizations/blockpreconditioner.jl b/src/factorizations/blockpreconditioner.jl index 1d3d036..a92fad1 100644 --- a/src/factorizations/blockpreconditioner.jl +++ b/src/factorizations/blockpreconditioner.jl @@ -100,13 +100,24 @@ end Base.eltype(p::BlockPreconditioner)=eltype(p.facts[1]) +""" + BlockPreconBuilder(;precs=UMFPACKPreconBuilder(), + partitioning = A -> [1:size(A,1)] + +Return callable object constructing a left block Jacobi preconditioner +from partition of unknowns. -Base.@kwdef struct EquationBlockPreconBuilder - precs=UMFPACKPrecs() - partitioning= [A -> 1:size(A,1)] +- `partitioning(A)` shall return a vector of AbstractVectors describing the indices of the partitions of the matrix. + For a matrix of size `n x n`, e.g. partitioning could be `[ 1:n÷2, (n÷2+1):n]` or [ 1:2:n, 2:2:n]. + +- `precs(A,p)` shall return a left precondioner for a matrix block. +""" +Base.@kwdef struct BlockPreconBuilder + precs=UMFPACKPreconBuilder() + partitioning= A -> [1:size(A,1)] end -function (blockprecs::EquationBlockPreconBuilder)(A,p) +function (blockprecs::BlockPreconBuilder)(A,p) (;precs, partitioning)=blockprecs factorization= A->precs(A,p)[1] bp=BlockPreconditioner(A;partitioning=partitioning(A), factorization) @@ -115,5 +126,4 @@ end """ Allow array for precs => different precoms -EquationBlockPreconBuilder """ diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index e22e898..3311ac4 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -1,42 +1,89 @@ +""" + UMFPACKPreconBuilder() +Return callable object constructing a formal left preconditioner from an LU factorization using UMFPACK to be passed +as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. +""" struct UMFPACKPreconBuilder end +@static if USE_GPL_LIBS (::UMFPACKPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) +end + +""" + SparspakPreconBuilder() +Return callable object constructing a formal left preconditioner from an LU factorization using Sparspak to be passed +as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. +""" struct SparspakPreconBuilder end (::SparspakPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) +""" + JacobiPreconBuilder() + +Return callable object constructing a left Jacobi preconditioner +to be passed as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. +""" struct JacobiPreconBuilder end (::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(JacobiPreconditioner(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) -struct ILUZeroBuilder end -(::ILUZeroBuilder)(A,p)=(ilu0(A),LinearAlgebra.I) +""" + ILUZeroPreconBuilder() +Return callable object constructing a left zero fill-in ILU preconditioner +using [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl) +""" +struct ILUZeroPreconBuilder end +(::ILUZeroPreconBuilder)(A,p)=(ilu0(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) +""" + ILUTPreconBuilder(; droptol=0.1) -Base.@kwdef struct ILUTBuilder +Return callable object constructing a left ILUT preconditioner +using [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl) +""" +Base.@kwdef struct ILUTPreconBuilder droptol::Float64=0.1 end -(::ILUTBuilder)(A,p)= error("import IncompleteLU.jl in order to use ILUTBuilder") +(::ILUTPreconBuilder)(A,p)= error("import IncompleteLU.jl in order to use ILUTBuilder") + + -struct SmoothedAggregationAMGBuilder{Tk} +""" + SmoothedAggregationPreconBuilder(;blocksize=1, kwargs...) + +Return callable object constructing a left smoothed aggregation algebraic multigrid preconditioner +using [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl). + +Needs `import AlgebraicMultigrid` to trigger the corresponding extension. +""" +struct SmoothedAggregationPreconBuilder{Tk} blocksize::Int kwargs::Tk end -function SmoothedAggregationAMGBuilder(;blocksize=1, kwargs...) - SmoothedAggregationAMGBuilder(blocksize,kwargs) +function SmoothedAggregationPreconBuilder(;blocksize=1, kwargs...) + SmoothedAggregationPreconBuilder(blocksize,kwargs) end -(::SmoothedAggregationAMGBuilder)(A,p)= error("import AlgebraicMultigrid in order to use SmoothedAggregationAMGBuilder") +(::SmoothedAggregationPreconBuilder)(A,p)= error("import AlgebraicMultigrid in order to use SmoothedAggregationPreconBuilder") + +""" + RugeStubenPreconBuilder(;blocksize=1, kwargs...) + +Return callable object constructing a left algebraic multigrid preconditioner after Ruge & Stüben +using [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl). -struct RugeStubenAMGBuilder{Tk} +Needs `import AlgebraicMultigrid` to trigger the corresponding extension. +""" +struct RugeStubenPreconBuilder{Tk} blocksize::Int kwargs::Tk end -function RugeStubenAMGBuilder(;blocksize=1, kwargs...) - SmoothedAggregationAMGBuilder(blocksize,kwargs) +function RugeStubenPreconBuilder(;blocksize=1, kwargs...) + SmoothedAggregationPreconBuilder(blocksize,kwargs) end -(::RugeStubenAMGBuilder)(A,p)= error("import AlgebraicMultigrid in order to use RugeStubenAMGBuilder") +(::RugeStubenPreconBuilder)(A,p)= error("import AlgebraicMultigrid in order to use RugeStubenAMGBuilder") diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index 3b8a0dd..bf3618d 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -81,10 +81,10 @@ allprecs=[ AMGCLWrap.AMGPreconBuilder(), AMGCLWrap.AMGPreconBuilder(), AMGCLWrap.RLXPreconBuilder(), - ExtendableSparse.ILUZeroBuilder(), - ExtendableSparse.ILUTBuilder(), - ExtendableSparse.SmoothedAggregationAMGBuilder(), - ExtendableSparse.RugeStubenAMGBuilder() + ExtendableSparse.ILUZeroPreconBuilder(), + ExtendableSparse.ILUTPreconBuilder(), + ExtendableSparse.SmoothedAggregationPreconBuilder(), + ExtendableSparse.RugeStubenPreconBuilder() ] @testset "iterations" begin @@ -109,7 +109,7 @@ end b=A*ones(n^2); for precs in allprecs - iteration=KrylovJL_CG(precs=EquationBlockPreconBuilder(;precs, partitioning)) + iteration=KrylovJL_CG(precs=BlockPreconBuilder(;precs, partitioning)) p=LinearProblem(A,b) sol=solve(p, KrylovJL_CG(;precs), abstol=1.0e-12) @test isapprox(sol, sol0, atol=1e-6) From bd1624d8080714630739ee47aab96923e8e66819 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 22:03:02 +0100 Subject: [PATCH 12/24] Documentation for new LinearSolve/precs API --- CHANGELOG.md | 2 +- docs/make.jl | 5 ++- docs/src/iter.md | 4 ++- docs/src/linearsolve.md | 74 ++++++++++++++++++++++++++++++++--------- 4 files changed, 64 insertions(+), 21 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 43fab43..16a6b8a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,7 @@ ### Breaking - remove solver + precon API which is not based on precs or directly overloading `\`. Fully rely on LinearSolve (besides `\`) -- Move AMGBuilder, ILUZeroBuilder to the correspondig packages (depending on the PRs) +- Move AMGBuilder, ILUZeroBuilder etc. to the correspondig packages (depending on the PRs) - remove "old" SparseMatrixLNK (need to benchmark before) ## [1.6.0] - WIP diff --git a/docs/make.jl b/docs/make.jl index 2eb1e8d..c9d39b1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,4 +1,3 @@ -push!(LOAD_PATH, "../src/") using Documenter, ExtendableSparse, AlgebraicMultigrid, IncompleteLU, Sparspak, LinearAlgebra function mkdocs() @@ -12,10 +11,10 @@ function mkdocs() pages = [ "Home" => "index.md", "example.md", - "linearsolve.md", "extsparse.md", - "iter.md", + "linearsolve.md", "internal.md", + "iter.md", "changes.md", ]) end diff --git a/docs/src/iter.md b/docs/src/iter.md index 22b68d0..c3b97dd 100644 --- a/docs/src/iter.md +++ b/docs/src/iter.md @@ -1,6 +1,8 @@ # Factorizations & Preconditioners -This functionality probably will be reduced in favor of LinearSolve.jl. +This functionality is deprecated as of version 1.6 an will be removed +with version 2.0. Use the [integration with LinearSolve.jl](/linearsolve/#Integration-with-LinearSolve.jl) + ## Factorizations diff --git a/docs/src/linearsolve.md b/docs/src/linearsolve.md index 22d87f6..f1d4539 100644 --- a/docs/src/linearsolve.md +++ b/docs/src/linearsolve.md @@ -3,12 +3,9 @@ Starting with version 0.9.6, ExtendableSparse is compatible with [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl). Since version 0.9.7, this is facilitated via the -AbstractSparseMatrixCSC interface. - -```@autodocs -Modules = [ExtendableSparse] -Pages = ["linearsolve.jl"] -``` +AbstractSparseMatrixCSC interface. Since version 1.6, +preconditioner constructors can be passed to iterative solvers via the [`precs` +keyword argument](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#prec). We can create a test problem and solve it with the `\` operator. @@ -21,7 +18,7 @@ y = A \ b sum(y) ``` -The same problem can be solved by the tools available via `LinearSolve.jl`: +The same problem can be solved via `LinearSolve.jl`: ```@example using ExtendableSparse # hide @@ -33,31 +30,76 @@ y = solve(LinearProblem(A, b), SparspakFactorization()).u sum(y) ``` -Also, the iterative method interface works with ExtendableSparse. + +## Preconditioning + +LinearSolve allows to pass preconditioner constructors via the `precs` keyword +to the iterative solver specification. ```@example using ExtendableSparse # hide using LinearSolve # hide -using SparseArrays # hide -using ILUZero # hide +using ExtendableSparse: ILUZeroPreconBuilder A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) x = ones(1000) b = A * x -y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(); - Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u +y = LinearSolve.solve(LinearProblem(A, b), + KrylovJL_CG(; precs=ILUZeroPreconBuilder())).u sum(y) ``` -However, ExtendableSparse provides a number of wrappers around preconditioners -from various Julia packages. +ExtendableSparse provides constructors for preconditioners wich can be used as `precs`. +These generally return a tuple `(Pl,I)` of a left preconditioner and a trivial right preconditioner. + +ExtendableSparse has a number of package extensions which construct preconditioners +from some other packages. In the future, these packages may provide this functionality on their own. + +```@docs +ExtendableSparse.ILUZeroPreconBuilder +ExtendableSparse.ILUTPreconBuilder +ExtendableSparse.SmoothedAggregationPreconBuilder +ExtendableSparse.RugeStubenPreconBuilder +``` + +In addition, ExtendableSparse implements some preconditioners: + +```@docs +ExtendableSparse.JacobiPreconBuilder +``` + +LU factorizations of matrices from previous iteration steps may be good +preconditioners for Krylov solvers called during a nonlinear solve via +Newton's method. For this purpose, ExtendableSparse provides preconditioner constructors +which wrap sparse LU factorizations. + +```@docs +ExtendableSparse.UMFPACKPreconBuilder +ExtendableSparse.SparspakPreconBuilder +``` + +Block preconditioner constructors are provided as well + +```@docs; canonical=false +ExtendableSparse.BlockPreconBuilder +``` + +## Deprecated API +Passing a preconditioner via the `Pl` or `Pr` keyword arguments +will be deprecated in LinearSolve. ExtendableSparse used to +export a number of wrappers for preconditioners from other packages +for this purpose. This approach is deprecated as of v1.6 and will be removed +with v2.0. + ```@example using ExtendableSparse # hide using LinearSolve # hide -using ILUZero # hide +using SparseArrays # hide +using ILUZero A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) x = ones(1000) b = A * x y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(); - Pl = ILUZeroPreconditioner(A)).u + Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u sum(y) ``` + From c400eca2425d870d870d90439724897eeea827d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 22:17:50 +0100 Subject: [PATCH 13/24] fix documentation references --- docs/src/extsparse.md | 4 ++++ docs/src/internal.md | 2 +- src/ExtendableSparse.jl | 23 +++++++++++++++++++++++ 3 files changed, 28 insertions(+), 1 deletion(-) diff --git a/docs/src/extsparse.md b/docs/src/extsparse.md index c4d2cd6..e756042 100644 --- a/docs/src/extsparse.md +++ b/docs/src/extsparse.md @@ -2,6 +2,10 @@ ## Matrix creation and update API +```@docs +ExtendableSparseMatrix +``` + ```@autodocs Modules = [ExtendableSparse] Pages = ["extendable.jl"] diff --git a/docs/src/internal.md b/docs/src/internal.md index 71a10f3..477058c 100644 --- a/docs/src/internal.md +++ b/docs/src/internal.md @@ -15,7 +15,7 @@ Pages = ["sparsematrixcsc.jl"] ``` ## New API Under development - aimed at multithreading -```@autodocs +```@autodocs; canonical = false Modules = [ExtendableSparse] Pages = ["abstractsparsematrixextension.jl", "abstractextendablesparsematrixcsc.jl", diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index cefa54e..70d3a99 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -32,10 +32,33 @@ include("matrix/extendable.jl") include("matrix/genericmtextendablesparsematrixcsc.jl") include("matrix/genericextendablesparsematrixcsc.jl") +""" + ExtendableSparseMatrix + +Aliased to [`ExtendableSparseMatrixCSC`](@ref) +""" const ExtendableSparseMatrix=ExtendableSparseMatrixCSC + + +""" + MTExtendableSparseMatrixCSC + +Multithreaded extendable sparse matrix (Experimental). + +Aliased to [`GenericMTExtendableSparseMatricCSC`](@ref) with [`SparseMatrixDILNKC`](@ref) +scalar matrix parameter. +""" const MTExtendableSparseMatrixCSC{Tv,Ti}=GenericMTExtendableSparseMatrixCSC{SparseMatrixDILNKC{Tv,Ti},Tv,Ti} MTExtendableSparseMatrixCSC(m,n,args...)=MTExtendableSparseMatrixCSC{Float64,Int64}(m,n,args...) +""" + STExtendableSparseMatrixCSC + +Single threaded extendable sparse matrix (Experimental). + +Aliased to [`GenericExtendableSparseMatricCSC`](@ref) with [`SparseMatrixDILNKC`](@ref) +scalar matrix parameter. +""" const STExtendableSparseMatrixCSC{Tv,Ti}=GenericExtendableSparseMatrixCSC{SparseMatrixDILNKC{Tv,Ti},Tv,Ti} STExtendableSparseMatrixCSC(m,n,args...)=STExtendableSparseMatrixCSC{Float64,Int64}(m,n,args...) From 19b6a4808bc2dc99d4287b04c7a0584fbaaf22d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 22:32:01 +0100 Subject: [PATCH 14/24] fix missing I --- src/factorizations/blockpreconditioner.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorizations/blockpreconditioner.jl b/src/factorizations/blockpreconditioner.jl index a92fad1..d488cfa 100644 --- a/src/factorizations/blockpreconditioner.jl +++ b/src/factorizations/blockpreconditioner.jl @@ -121,7 +121,7 @@ function (blockprecs::BlockPreconBuilder)(A,p) (;precs, partitioning)=blockprecs factorization= A->precs(A,p)[1] bp=BlockPreconditioner(A;partitioning=partitioning(A), factorization) - (bp,I) + (bp,LinearAlgebra.I) end """ From 8a8cd60ed0845a20cabcf154836452ecd855bca6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 23:23:36 +0100 Subject: [PATCH 15/24] Add blocksize to ILUZeroPrecon uses some type piracy which can be easily circumvented --- src/preconbuilders.jl | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 3311ac4..ffa8d26 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -29,13 +29,33 @@ struct JacobiPreconBuilder end """ - ILUZeroPreconBuilder() + ILUZeroPreconBuilder(;blocksize=1) Return callable object constructing a left zero fill-in ILU preconditioner using [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl) """ -struct ILUZeroPreconBuilder end -(::ILUZeroPreconBuilder)(A,p)=(ilu0(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) +Base.@kwdef struct ILUZeroPreconBuilder + blocksize::Int = 1 +end + +function (b::ILUZeroPreconBuilder)(A0,p) + A=SparseMatrixCSC(size(A0)..., getcolptr(A0), rowvals(A0),nonzeros(A0)) + if b.blocksize==1 + (ILUZero.ilu0(A),LinearAlgebra.I) + else + (ILUZero.ilu0(pointblock(A,b.blocksize),SVector{b.blocksize,eltype(A)}),LinearAlgebra.I) + end +end + +# Harrr!!! ☠ +function LinearAlgebra.ldiv!(Y::Vector{Tv}, + A::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}}, + B::Vector{Tv}) where {N,NN,Tv,Ti} + BY=reinterpret(SVector{N,Tv},Y) + BB=reinterpret(SVector{N,Tv},B) + ldiv!(BY,A,BB) + Y +end """ ILUTPreconBuilder(; droptol=0.1) From db14464fd2ebab98a6763fd36470c2ce119d54a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 23:27:09 +0100 Subject: [PATCH 16/24] Deprecation warnings Also, Ensure that they occur only once. --- ext/ExtendableSparseAlgebraicMultigridExt.jl | 14 ++++++++++++++ ext/ExtendableSparseIncompleteLUExt.jl | 6 ++++++ src/factorizations/iluzero.jl | 12 ++++++++++++ src/preconbuilders.jl | 1 + test/test_linearsolve.jl | 1 + 5 files changed, 34 insertions(+) diff --git a/ext/ExtendableSparseAlgebraicMultigridExt.jl b/ext/ExtendableSparseAlgebraicMultigridExt.jl index 9d69470..7abdbc2 100644 --- a/ext/ExtendableSparseAlgebraicMultigridExt.jl +++ b/ext/ExtendableSparseAlgebraicMultigridExt.jl @@ -19,12 +19,19 @@ import ExtendableSparse: RugeStubenPreconBuilder import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update! ###################################################################################### +rswarned=false + mutable struct RS_AMGPreconditioner <: AbstractPreconditioner A::ExtendableSparseMatrix factorization::AlgebraicMultigrid.Preconditioner kwargs blocksize function ExtendableSparse.RS_AMGPreconditioner(blocksize=1; kwargs...) + global rswarned + if !rswarned + @warn "RS_AMGPreconditioner is deprecated. Use LinearSolve with `precs=RugeStubenPreconBuilder()` instead" + rswarned=true + end precon = new() precon.kwargs = kwargs precon.blocksize=blocksize @@ -45,13 +52,20 @@ end allow_views(::RS_AMGPreconditioner)=true allow_views(::Type{RS_AMGPreconditioner})=true + ###################################################################################### +sawarned=false mutable struct SA_AMGPreconditioner <: AbstractPreconditioner A::ExtendableSparseMatrix factorization::AlgebraicMultigrid.Preconditioner kwargs blocksize function ExtendableSparse.SA_AMGPreconditioner(blocksize=1; kwargs...) + global sawarned + if !sawarned + @warn "SA_AMGPreconditioner is deprecated. Use LinearSolve with `precs=SmoothedAggregationPreconBuilder()` instead" + sawarned=true + end precon = new() precon.kwargs = kwargs precon.blocksize=blocksize diff --git a/ext/ExtendableSparseIncompleteLUExt.jl b/ext/ExtendableSparseIncompleteLUExt.jl index 6074574..bfeb7ed 100644 --- a/ext/ExtendableSparseIncompleteLUExt.jl +++ b/ext/ExtendableSparseIncompleteLUExt.jl @@ -13,11 +13,17 @@ import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update! # Deprecated from here +warned=false mutable struct ILUTPreconditioner <: AbstractPreconditioner A::ExtendableSparseMatrix factorization::IncompleteLU.ILUFactorization droptol::Float64 function ExtendableSparse.ILUTPreconditioner(; droptol = 1.0e-3) + global warned + if !warned + @warn "ILUTPreconditioner is deprecated. Use LinearSolve with `precs=ILUTPreconBuilder()` instead" + warned=true + end p = new() p.droptol = droptol p diff --git a/src/factorizations/iluzero.jl b/src/factorizations/iluzero.jl index 2f068c5..efc7ef7 100644 --- a/src/factorizations/iluzero.jl +++ b/src/factorizations/iluzero.jl @@ -1,8 +1,14 @@ +iluzerowarned=false mutable struct ILUZeroPreconditioner <: AbstractPreconditioner A::ExtendableSparseMatrix factorization::ILUZero.ILU0Precon phash::UInt64 function ILUZeroPreconditioner() + global iluzerowarned + if !iluzerowarned + @warn "ILUZeroPreconditioner is deprecated. Use LinearSolve with `precs=ILUZeroPreconBuilder()` instead" + iluzerowarned=true + end p = new() p.phash = 0 p @@ -35,12 +41,18 @@ allow_views(::Type{ILUZeroPreconditioner})=true +biluzerowarned=false mutable struct PointBlockILUZeroPreconditioner <: AbstractPreconditioner A::ExtendableSparseMatrix factorization::ILUZero.ILU0Precon phash::UInt64 blocksize::Int function PointBlockILUZeroPreconditioner(;blocksize=1) + global biluzerowarned + if !biluzerowarned + @warn "PointBlockILUZeroPreconditioner is deprecated. Use LinearSolve with `precs=ILUZeroPreconBuilder(; blocksize=$(blocksize))` instead" + biluzerowarned=true + end p = new() p.phash = 0 p.blocksize=blocksize diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index ffa8d26..4496c25 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -48,6 +48,7 @@ function (b::ILUZeroPreconBuilder)(A0,p) end # Harrr!!! ☠ +# We could resolve this piracy by introducing a wrapper type for the block case. function LinearAlgebra.ldiv!(Y::Vector{Tv}, A::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}}, B::Vector{Tv}) where {N,NN,Tv,Ti} diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index bf3618d..aed53f2 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -82,6 +82,7 @@ allprecs=[ AMGCLWrap.AMGPreconBuilder(), AMGCLWrap.RLXPreconBuilder(), ExtendableSparse.ILUZeroPreconBuilder(), + ExtendableSparse.ILUZeroPreconBuilder(;blocksize=2), ExtendableSparse.ILUTPreconBuilder(), ExtendableSparse.SmoothedAggregationPreconBuilder(), ExtendableSparse.RugeStubenPreconBuilder() From d60e3815bde50bcc7cbb1438bd2ce042945b6e41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sat, 9 Nov 2024 23:58:35 +0100 Subject: [PATCH 17/24] Cholesky: PreconBuilder, deprecation --- src/factorizations/cholmod_cholesky.jl | 5 +++++ src/preconbuilders.jl | 16 ++++++++++++++++ test/test_linearsolve.jl | 8 +++++++- 3 files changed, 28 insertions(+), 1 deletion(-) diff --git a/src/factorizations/cholmod_cholesky.jl b/src/factorizations/cholmod_cholesky.jl index 1e3fbc8..578d749 100644 --- a/src/factorizations/cholmod_cholesky.jl +++ b/src/factorizations/cholmod_cholesky.jl @@ -5,6 +5,7 @@ mutable struct CholeskyFactorization <: AbstractLUFactorization A64::Any end +cholwarned=false """ CholeskyFactorization(;valuetype=Float64, indextype=Int64) CholeskyFactorization(matrix) @@ -12,6 +13,10 @@ CholeskyFactorization(matrix) Default Cholesky factorization via cholmod. """ function CholeskyFactorization() + if !cholwarned + @warn "ExtendableSparse.CholeskyFactorization is deprecated. Use LinearSolve.CholeskyFactorization` instead" + cholwarned=true + end CholeskyFactorization(nothing, nothing, 0, nothing) end diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 4496c25..97c5cf3 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -9,6 +9,22 @@ struct UMFPACKPreconBuilder end (::UMFPACKPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) end +""" + CholeskyPreconBuilder() + +Return callable object constructing a formal left preconditioner from a Cholesky factorization using CHOLMOD to be passed +as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. +""" +struct CholeskyPreconBuilder end +@static if USE_GPL_LIBS +(::CholeskyPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(cholesky(Symmetric(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A)))),LinearAlgebra.I) + +# Harrr!!! ☠ +LinearAlgebra.ldiv!(fact::SparseArrays.CHOLMOD.Factor, v) = fact \ v + +end + + """ SparspakPreconBuilder() diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index aed53f2..517c988 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -84,6 +84,7 @@ allprecs=[ ExtendableSparse.ILUZeroPreconBuilder(), ExtendableSparse.ILUZeroPreconBuilder(;blocksize=2), ExtendableSparse.ILUTPreconBuilder(), + ExtendableSparse.JacobiPreconBuilder(), ExtendableSparse.SmoothedAggregationPreconBuilder(), ExtendableSparse.RugeStubenPreconBuilder() ] @@ -102,6 +103,11 @@ allprecs=[ end end + +moreprecs=[ExtendableSparse.UMFPACKPreconBuilder(), + ExtendableSparse.SparspakPreconBuilder(), + ExtendableSparse.CholeskyPreconBuilder()] + @testset "equationblock" begin n=100 A=fdrand(n,n) @@ -109,7 +115,7 @@ end sol0=ones(n^2) b=A*ones(n^2); - for precs in allprecs + for precs in vcat(allprecs, moreprecs) iteration=KrylovJL_CG(precs=BlockPreconBuilder(;precs, partitioning)) p=LinearProblem(A,b) sol=solve(p, KrylovJL_CG(;precs), abstol=1.0e-12) From ffd67ddad2988e5f12b38975dfba267136383df8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 10 Nov 2024 12:11:39 +0100 Subject: [PATCH 18/24] Add LinearSolve extension --- Project.toml | 3 +++ ext/ExtendableSparseLinearSolveExt.jl | 23 +++++++++++++++++++++++ 2 files changed, 26 insertions(+) create mode 100644 ext/ExtendableSparseLinearSolveExt.jl diff --git a/Project.toml b/Project.toml index 20345a0..179515d 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288" AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" [extensions] @@ -26,6 +27,7 @@ ExtendableSparseAMGCLWrapExt = "AMGCLWrap" ExtendableSparseAlgebraicMultigridExt = "AlgebraicMultigrid" ExtendableSparseIncompleteLUExt = "IncompleteLU" ExtendableSparsePardisoExt = "Pardiso" +ExtendableSparseLinearSolveExt = "LinearSolve" [compat] AMGCLWrap = "2" @@ -33,6 +35,7 @@ AlgebraicMultigrid = "0.4,0.5,0.6" DocStringExtensions = "0.8, 0.9" ILUZero = "0.2" IncompleteLU = "^0.2.1" +LinearSolve = "2.36.0" Pardiso = "0.5.1" Sparspak = "0.3.6" StaticArrays = "1.5.24" diff --git a/ext/ExtendableSparseLinearSolveExt.jl b/ext/ExtendableSparseLinearSolveExt.jl new file mode 100644 index 0000000..6da6894 --- /dev/null +++ b/ext/ExtendableSparseLinearSolveExt.jl @@ -0,0 +1,23 @@ +module ExtendableSparseLinearSolveExt +using LinearSolve +import ExtendableSparse: LinearSolvePreconBuilder +import LinearAlgebra +using SparseArrays: AbstractSparseMatrixCSC + +# Harrr! +# Avoid type piracy by adding a wrapper struct +function (method::LinearSolve.AbstractFactorization)(A) + pr = LinearProblem(A, zeros(eltype(A), size(A, 1))) + init(pr, method) +end + +function LinearAlgebra.ldiv!(u, cache::LinearSolve.LinearCache, b) + cache.b = b + sol = solve!(cache) + copyto!(u, sol.u) +end + +(b::LinearSolvePreconBuilder)(A::AbstractSparseMatrixCSC,p) = (b.method(A), LinearAlgebra.I) + +end + From d4d86931f84e389c777eaca476cc1b90d07c5811 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 10 Nov 2024 12:12:29 +0100 Subject: [PATCH 19/24] Introduce LinearSolvePreconBuilder to handle preconditioners from LU factorization Replaces UMFPACKPreconBuilder and the like --- docs/src/linearsolve.md | 9 +++------ src/ExtendableSparse.jl | 2 +- src/preconbuilders.jl | 36 ++++++------------------------------ test/test_linearsolve.jl | 9 ++++----- 4 files changed, 14 insertions(+), 42 deletions(-) diff --git a/docs/src/linearsolve.md b/docs/src/linearsolve.md index f1d4539..22a4066 100644 --- a/docs/src/linearsolve.md +++ b/docs/src/linearsolve.md @@ -69,16 +69,13 @@ ExtendableSparse.JacobiPreconBuilder LU factorizations of matrices from previous iteration steps may be good preconditioners for Krylov solvers called during a nonlinear solve via -Newton's method. For this purpose, ExtendableSparse provides preconditioner constructors -which wrap sparse LU factorizations. - +Newton's method. For this purpose, ExtendableSparse provides a preconditioner constructor +which wraps sparse LU factorizations supported by LinearSolve.jl ```@docs -ExtendableSparse.UMFPACKPreconBuilder -ExtendableSparse.SparspakPreconBuilder +ExtendableSparse.LinearSolvePreconBuilder ``` Block preconditioner constructors are provided as well - ```@docs; canonical=false ExtendableSparse.BlockPreconBuilder ``` diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 70d3a99..e4a67c8 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -77,7 +77,7 @@ include("factorizations/simple_iteration.jl") export simple, simple! include("preconbuilders.jl") -export SparspakPreconBuilder, UMFPACKPreconBuilder, BlockPreconBuilder, JacobiPreconBuilder +export LinearSolvePreconBuilder, BlockPreconBuilder, JacobiPreconBuilder @public ILUZeroPreconBuilder, ILUTPreconBuilder, SmoothedAggregationPreconBuilder, RugeStubenPreconBuilder diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 97c5cf3..a4bc5c7 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -1,38 +1,14 @@ """ - UMFPACKPreconBuilder() + LinearSolvePreconBuilder(; method=UMFPACKFactorization()) -Return callable object constructing a formal left preconditioner from an LU factorization using UMFPACK to be passed -as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. +Return callable object constructing a formal left preconditioner from a sparse LU factorization using LinearSolve +as the `precs` parameter for a [`BlockPreconBuilder`](@ref) or iterative methods wrapped by LinearSolve.jl. """ -struct UMFPACKPreconBuilder end -@static if USE_GPL_LIBS -(::UMFPACKPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) +Base.@kwdef struct LinearSolvePreconBuilder + method=UMFPACKFactorization() end +(::LinearSolvePreconBuilder)(A,p)= error("import LinearSolve in order to use LinearSolvePreconBuilder") -""" - CholeskyPreconBuilder() - -Return callable object constructing a formal left preconditioner from a Cholesky factorization using CHOLMOD to be passed -as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. -""" -struct CholeskyPreconBuilder end -@static if USE_GPL_LIBS -(::CholeskyPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(cholesky(Symmetric(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A)))),LinearAlgebra.I) - -# Harrr!!! ☠ -LinearAlgebra.ldiv!(fact::SparseArrays.CHOLMOD.Factor, v) = fact \ v - -end - - -""" - SparspakPreconBuilder() - -Return callable object constructing a formal left preconditioner from an LU factorization using Sparspak to be passed -as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. -""" -struct SparspakPreconBuilder end -(::SparspakPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I) """ JacobiPreconBuilder() diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index 517c988..57355be 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -57,6 +57,7 @@ end end factorizations=[UMFPACKFactorization(), + SparspakFactorization(), KLUFactorization(reuse_symbolic=false)] if !Sys.isapple() @@ -104,18 +105,16 @@ allprecs=[ end -moreprecs=[ExtendableSparse.UMFPACKPreconBuilder(), - ExtendableSparse.SparspakPreconBuilder(), - ExtendableSparse.CholeskyPreconBuilder()] +luprecs=[ExtendableSparse.LinearSolvePreconBuilder(factorization) for factorization in factorizations] -@testset "equationblock" begin +@testset "block preconditioning" begin n=100 A=fdrand(n,n) partitioning=A->[1:2:size(A,1), 2:2:size(A,1)] sol0=ones(n^2) b=A*ones(n^2); - for precs in vcat(allprecs, moreprecs) + for precs in vcat(allprecs, luprecs) iteration=KrylovJL_CG(precs=BlockPreconBuilder(;precs, partitioning)) p=LinearProblem(A,b) sol=solve(p, KrylovJL_CG(;precs), abstol=1.0e-12) From ccbb2ad144ed50f40603ee9c6d472f545aa43a9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 10 Nov 2024 12:13:43 +0100 Subject: [PATCH 20/24] fix more deprecations --- ext/ExtendableSparsePardisoExt.jl | 2 +- src/factorizations/cholmod_cholesky.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/ExtendableSparsePardisoExt.jl b/ext/ExtendableSparsePardisoExt.jl index 4faaa2f..69b8e29 100644 --- a/ext/ExtendableSparsePardisoExt.jl +++ b/ext/ExtendableSparsePardisoExt.jl @@ -2,7 +2,7 @@ module ExtendableSparsePardisoExt using ExtendableSparse using LinearAlgebra using Pardiso - +# Deprecated import ExtendableSparse: @makefrommatrix, update!, AbstractLUFactorization abstract type AbstractPardisoLU <: AbstractLUFactorization end diff --git a/src/factorizations/cholmod_cholesky.jl b/src/factorizations/cholmod_cholesky.jl index 578d749..66228b4 100644 --- a/src/factorizations/cholmod_cholesky.jl +++ b/src/factorizations/cholmod_cholesky.jl @@ -13,6 +13,7 @@ CholeskyFactorization(matrix) Default Cholesky factorization via cholmod. """ function CholeskyFactorization() + global cholwarned if !cholwarned @warn "ExtendableSparse.CholeskyFactorization is deprecated. Use LinearSolve.CholeskyFactorization` instead" cholwarned=true From 868991cd0be2d0cccfb8214671d4caf27aad2ddf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 10 Nov 2024 12:23:21 +0100 Subject: [PATCH 21/24] fix sentence point --- docs/src/iter.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/iter.md b/docs/src/iter.md index c3b97dd..31030c0 100644 --- a/docs/src/iter.md +++ b/docs/src/iter.md @@ -1,7 +1,7 @@ # Factorizations & Preconditioners This functionality is deprecated as of version 1.6 an will be removed -with version 2.0. Use the [integration with LinearSolve.jl](/linearsolve/#Integration-with-LinearSolve.jl) +with version 2.0. Use the [integration with LinearSolve.jl](/linearsolve/#Integration-with-LinearSolve.jl). ## Factorizations From 659924be1eaa73075e6203bf28f113a7bcc0dc1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 10 Nov 2024 13:03:02 +0100 Subject: [PATCH 22/24] fix newly introduced type piracies --- ext/ExtendableSparseLinearSolveExt.jl | 19 +++++++++++-------- src/preconbuilders.jl | 25 ++++++++++++++----------- 2 files changed, 25 insertions(+), 19 deletions(-) diff --git a/ext/ExtendableSparseLinearSolveExt.jl b/ext/ExtendableSparseLinearSolveExt.jl index 6da6894..2c5335a 100644 --- a/ext/ExtendableSparseLinearSolveExt.jl +++ b/ext/ExtendableSparseLinearSolveExt.jl @@ -4,20 +4,23 @@ import ExtendableSparse: LinearSolvePreconBuilder import LinearAlgebra using SparseArrays: AbstractSparseMatrixCSC -# Harrr! -# Avoid type piracy by adding a wrapper struct -function (method::LinearSolve.AbstractFactorization)(A) + +struct LinearSolvePrecon{T} + cache::T +end + +function LinearSolvePrecon(A,method::LinearSolve.AbstractFactorization) pr = LinearProblem(A, zeros(eltype(A), size(A, 1))) - init(pr, method) + LinearSolvePrecon(init(pr, method)) end -function LinearAlgebra.ldiv!(u, cache::LinearSolve.LinearCache, b) - cache.b = b - sol = solve!(cache) +function LinearAlgebra.ldiv!(u, P::LinearSolvePrecon, b) + P.cache.b = b + sol = solve!(P.cache) copyto!(u, sol.u) end -(b::LinearSolvePreconBuilder)(A::AbstractSparseMatrixCSC,p) = (b.method(A), LinearAlgebra.I) +(b::LinearSolvePreconBuilder)(A::AbstractSparseMatrixCSC,p) = (LinearSolvePrecon(A,b.method), LinearAlgebra.I) end diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index a4bc5c7..dac9127 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -30,26 +30,29 @@ Base.@kwdef struct ILUZeroPreconBuilder blocksize::Int = 1 end -function (b::ILUZeroPreconBuilder)(A0,p) - A=SparseMatrixCSC(size(A0)..., getcolptr(A0), rowvals(A0),nonzeros(A0)) - if b.blocksize==1 - (ILUZero.ilu0(A),LinearAlgebra.I) - else - (ILUZero.ilu0(pointblock(A,b.blocksize),SVector{b.blocksize,eltype(A)}),LinearAlgebra.I) - end +struct ILUBlockPrecon{N,NN,Tv,Ti} + ilu0::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}} end -# Harrr!!! ☠ -# We could resolve this piracy by introducing a wrapper type for the block case. function LinearAlgebra.ldiv!(Y::Vector{Tv}, - A::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}}, + A::ILUBlockPrecon{N,NN,Tv,Ti}, B::Vector{Tv}) where {N,NN,Tv,Ti} BY=reinterpret(SVector{N,Tv},Y) BB=reinterpret(SVector{N,Tv},B) - ldiv!(BY,A,BB) + ldiv!(BY,A.ilu0,BB) Y end +function (b::ILUZeroPreconBuilder)(A0,p) + A=SparseMatrixCSC(size(A0)..., getcolptr(A0), rowvals(A0),nonzeros(A0)) + if b.blocksize==1 + (ILUZero.ilu0(A),LinearAlgebra.I) + else + (ILUBlockPrecon(ILUZero.ilu0(pointblock(A,b.blocksize),SVector{b.blocksize,eltype(A)})),LinearAlgebra.I) + end +end + + """ ILUTPreconBuilder(; droptol=0.1) From d900f7c880ffd960f633d4e9cc4f06f43836141e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 10 Nov 2024 15:06:48 +0100 Subject: [PATCH 23/24] Make block precon builders mutable Blocksizes & partionings sometimes not known when creating the strategy. --- src/factorizations/blockpreconditioner.jl | 2 +- src/preconbuilders.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorizations/blockpreconditioner.jl b/src/factorizations/blockpreconditioner.jl index d488cfa..e9bc1f8 100644 --- a/src/factorizations/blockpreconditioner.jl +++ b/src/factorizations/blockpreconditioner.jl @@ -112,7 +112,7 @@ from partition of unknowns. - `precs(A,p)` shall return a left precondioner for a matrix block. """ -Base.@kwdef struct BlockPreconBuilder +Base.@kwdef mutable struct BlockPreconBuilder precs=UMFPACKPreconBuilder() partitioning= A -> [1:size(A,1)] end diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index dac9127..9cd5782 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -26,7 +26,7 @@ struct JacobiPreconBuilder end Return callable object constructing a left zero fill-in ILU preconditioner using [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl) """ -Base.@kwdef struct ILUZeroPreconBuilder +Base.@kwdef mutable struct ILUZeroPreconBuilder blocksize::Int = 1 end From 02cf5cf0e02c65d3fee232b6f5eda5c5a2c36405 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Sun, 10 Nov 2024 19:14:23 +0100 Subject: [PATCH 24/24] udate docs --- Project.toml | 2 +- docs/Project.toml | 1 + docs/src/linearsolve.md | 90 ++++++++++++++++++++++++++++++++--------- 3 files changed, 72 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 179515d..eabf350 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableSparse" uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" authors = ["Juergen Fuhrmann "] -version = "1.6.0-dev" +version = "1.6.0" [deps] AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288" diff --git a/docs/Project.toml b/docs/Project.toml index 6179c00..5a3ae26 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,6 +9,7 @@ ILUZero = "88f59080-6952-5380-9ea5-54057fb9a43f" IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" [compat] diff --git a/docs/src/linearsolve.md b/docs/src/linearsolve.md index 22a4066..ffaa9de 100644 --- a/docs/src/linearsolve.md +++ b/docs/src/linearsolve.md @@ -1,16 +1,12 @@ -# Integration with LinearSolve.jl +# Linear System solution -Starting with version 0.9.6, ExtendableSparse is compatible -with [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl). -Since version 0.9.7, this is facilitated via the -AbstractSparseMatrixCSC interface. Since version 1.6, -preconditioner constructors can be passed to iterative solvers via the [`precs` -keyword argument](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#prec). - -We can create a test problem and solve it with the `\` operator. +## The `\` operator +The packages overloads `\` for the ExtendableSparseMatrix type. +The following example uses [`fdrand`](@ref) to create a test matrix and solve +the corresponding linear system of equations. ```@example -using ExtendableSparse # hide +using ExtendableSparse A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) x = ones(1000) b = A * x @@ -18,27 +14,59 @@ y = A \ b sum(y) ``` +This works as well for number types besides `Float64` and related, in this case, +by default a LU factorization based on Sparspak ist used. + +```@example +using ExtendableSparse +using MultiFloats +A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix) +x = ones(Float64x2,1000) +b = A * x +y = A \ b +sum(y) +``` + +## Solving with LinearSolve.jl + +Starting with version 0.9.6, ExtendableSparse is compatible +with [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl). +Since version 0.9.7, this is facilitated via the +AbstractSparseMatrixCSC interface. + + The same problem can be solved via `LinearSolve.jl`: ```@example -using ExtendableSparse # hide -using LinearSolve # hide +using ExtendableSparse +using LinearSolve A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) x = ones(1000) b = A * x -y = solve(LinearProblem(A, b), SparspakFactorization()).u +y = solve(LinearProblem(A, b)).u sum(y) ``` +```@example +using ExtendableSparse +using LinearSolve +using MultiFloats +A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix) +x = ones(Float64x2,1000) +b = A * x +y = solve(LinearProblem(A, b), SparspakFactorization()).u +sum(y) +``` -## Preconditioning +## Preconditioned Krylov solvers with LinearSolve.jl -LinearSolve allows to pass preconditioner constructors via the `precs` keyword +Since version 1.6, preconditioner constructors can be passed to iterative solvers via the [`precs` +keyword argument](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#prec) to the iterative solver specification. ```@example -using ExtendableSparse # hide -using LinearSolve # hide +using ExtendableSparse +using LinearSolve using ExtendableSparse: ILUZeroPreconBuilder A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) x = ones(1000) @@ -48,6 +76,7 @@ y = LinearSolve.solve(LinearProblem(A, b), sum(y) ``` +## Available preconditioners ExtendableSparse provides constructors for preconditioners wich can be used as `precs`. These generally return a tuple `(Pl,I)` of a left preconditioner and a trivial right preconditioner. @@ -75,11 +104,32 @@ which wraps sparse LU factorizations supported by LinearSolve.jl ExtendableSparse.LinearSolvePreconBuilder ``` + Block preconditioner constructors are provided as well ```@docs; canonical=false ExtendableSparse.BlockPreconBuilder ``` + +The example beloww shows how to create a block Jacobi preconditioner where the blocks are defined by even and odd +degrees of freedom, and the diagonal blocks are solved using UMFPACK. +```@example +using ExtendableSparse +using LinearSolve +using ExtendableSparse: LinearSolvePreconBuilder, BlockPreconBuilder +A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) +x = ones(1000) +b = A * x +partitioning=A->[1:2:size(A,1), 2:2:size(A,1)] +umfpackprecs=LinearSolvePreconBuilder(UMFPACKFactorization()) +blockprecs=BlockPreconBuilder(;precs=umfpackprecs, partitioning) +y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(; precs=blockprecs)).u +sum(y) +``` +`umpfackpreconbuilder` e.g. could be replaced by `SmoothedAggregationPreconBuilder()`. Moreover, this approach +works for any `AbstractSparseMatrixCSC`. + + ## Deprecated API Passing a preconditioner via the `Pl` or `Pr` keyword arguments will be deprecated in LinearSolve. ExtendableSparse used to @@ -88,9 +138,9 @@ for this purpose. This approach is deprecated as of v1.6 and will be removed with v2.0. ```@example -using ExtendableSparse # hide -using LinearSolve # hide -using SparseArrays # hide +using ExtendableSparse +using LinearSolve +using SparseArray using ILUZero A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix) x = ones(1000)