diff --git a/.codespellrc b/.codespellrc index 7505bec..7b63a75 100644 --- a/.codespellrc +++ b/.codespellrc @@ -1,2 +1,2 @@ [codespell] -ignore-words-list = missings,rcall,linke,fo,coo,alledges +ignore-words-list = missings,rcall,linke,fo,coo,alledges,bu diff --git a/CHANGELOG.md b/CHANGELOG.md index 5ca57a8..530572e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,10 @@ # Changelog +## [2.4.0] - 2026-07-12 +- JacobiPreconBuilder now allows for blocksize +- More tests +- Internal reorganization (simplified BlockPreconditioner etc.) + ## [2.3.0] - 2026-06-16 - Product and Identity preconditioners and PreconBuilders - Block preconditioner and PreconBuilder now allow for vector @@ -13,7 +18,6 @@ - fix allocations in lnk+csc matrix addition ## [2.0.0] - 2026-01-06 - - Remove solver + precon API which is not based on precs or directly overloading `\`. Fully rely on LinearSolve (besides `\`) - Move AMGBuilder etc to corresponding packages (depending on the PRs) diff --git a/Project.toml b/Project.toml index 4295113..39c6514 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ExtendableSparse" uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" -version = "2.3.1" +version = "2.4.0" authors = ["Juergen Fuhrmann ", "Daniel Runge"] [deps] diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index cc55a8e..06efec9 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -41,6 +41,7 @@ export flush!, updateindex!, rawupdateindex!, reset!, nnznew export eliminate_dirichlet, eliminate_dirichlet!, mark_dirichlet +include("preconditioners.jl") include("preconbuilders.jl") export LinearSolvePreconBuilder, BlockPreconBuilder, JacobiPreconBuilder export ProductPreconBuilder, IdentityPreconBuilder diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 79dbb86..71854a4 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -1,3 +1,9 @@ +# +# This file defines `preconbuilders` which allow to specify `precs` preconditioner constructors +# for iterative solvers as they are handled in LinearSolve.jl. +# See https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#Specifying-Preconditioners +# + """ LinearSolvePreconBuilder(; method=UMFPACKFactorization()) @@ -11,17 +17,7 @@ end """ - 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) - - -""" - ILUZeroPreconBuilder(;blocksize=1) + ILUZeroPreconBuilder(; blocksize = 1) Return callable object constructing a left zero fill-in ILU preconditioner using [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl) @@ -30,66 +26,17 @@ Base.@kwdef mutable struct ILUZeroPreconBuilder blocksize::Int = 1 end -struct ILUBlockPrecon{N, NN, Tv, Ti} - ilu0::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}} -end - -function LinearAlgebra.ldiv!( - Y::Vector{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.ilu0, BB) - return Y -end - -""" - pointblock(matrix,blocksize) - -Create a pointblock matrix. -""" -function pointblock(A0::ExtendableSparseMatrixCSC{Tv, Ti}, blocksize) where {Tv, Ti} - A = SparseMatrixCSC(A0) - colptr = A.colptr - rowval = A.rowval - nzval = A.nzval - n = A.n - block = zeros(Tv, blocksize, blocksize) - nblock = n ÷ blocksize - b = SMatrix{blocksize, blocksize}(block) - Tb = typeof(b) - Ab = ExtendableSparseMatrixCSC{Tb, Ti}(nblock, nblock) - - - for i in 1:n - for k in colptr[i]:(colptr[i + 1] - 1) - j = rowval[k] - iblock = (i - 1) ÷ blocksize + 1 - jblock = (j - 1) ÷ blocksize + 1 - ii = (i - 1) % blocksize + 1 - jj = (j - 1) % blocksize + 1 - block[ii, jj] = nzval[k] - rawupdateindex!(Ab, +, SMatrix{blocksize, blocksize}(block), iblock, jblock) - block[ii, jj] = zero(Tv) - end - end - return flush!(Ab) -end - function (b::ILUZeroPreconBuilder)(A0, p) A = SparseMatrixCSC(size(A0)..., getcolptr(A0), rowvals(A0), nonzeros(A0)) return if b.blocksize == 1 (ILUZero.ilu0(A), LinearAlgebra.I) else - (ILUBlockPrecon(ILUZero.ilu0(pointblock(A, b.blocksize), SVector{b.blocksize, eltype(A)})), LinearAlgebra.I) + (ILU0BlockPrecon(ILUZero.ilu0(pointblock(A, b.blocksize), SVector{b.blocksize, eltype(A)})), LinearAlgebra.I) end end - """ - ILUTPreconBuilder(; droptol=0.1) + ILUTPreconBuilder(; droptol = 0.1) Return callable object constructing a left ILUT preconditioner using [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl) @@ -100,103 +47,6 @@ end (::ILUTPreconBuilder)(A, p) = error("import IncompleteLU.jl in order to use ILUTBuilder") -mutable struct BlockPreconditioner - A::AbstractMatrix - factorizations - partitioning::Union{Nothing, Vector{AbstractVector}} - facts::Vector - function BlockPreconditioner(A; partitioning = nothing, factorization = nothing, factorizations = nothing) - p = new() - p.A = A - p.partitioning = partitioning - if !isnothing(factorization) - p.factorizations = [factorization for i in 1:length(partitioning)] - else - p.factorizations = factorizations - end - update!(p) - return p - end -end - - -""" - BlockPreconditioner(;partitioning, factorization) - -Create a block preconditioner from partition of unknowns given by `partitioning`, 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]. -Factorization is a callable (Function or struct) which allows to create a factorization (with `ldiv!` methods) from a submatrix of A. -""" -function BlockPreconditioner end - -""" - allow_views(::preconditioner_type) - -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. -""" -allow_views(::Any) = false - - -function update!(precon::BlockPreconditioner) - flush!(precon.A) - nall = sum(length, precon.partitioning) - n = size(precon.A, 1) - if nall != n - @warn "sum(length,partitioning)=$(nall) but n=$(n)" - end - - if isnothing(precon.partitioning) - partitioning = [1:n] - end - - np = length(precon.partitioning) - precon.facts = Vector{Any}(undef, np) - return Threads.@threads for ipart in 1:np - AP = precon.A[precon.partitioning[ipart], precon.partitioning[ipart]] - precon.facts[ipart] = precon.factorizations[ipart](AP) - end -end - - -function LinearAlgebra.ldiv!(p::BlockPreconditioner, v) - partitioning = p.partitioning - facts = p.facts - np = length(partitioning) - - Threads.@threads for ipart in 1:np - if allow_views(p.factorizations[ipart]) - ldiv!(facts[ipart], view(v, partitioning[ipart])) - else - vv = v[partitioning[ipart]] - ldiv!(facts[ipart], vv) - view(v, partitioning[ipart]) .= vv - end - end - return v -end - -function LinearAlgebra.ldiv!(u, p::BlockPreconditioner, v) - partitioning = p.partitioning - facts = p.facts - np = length(partitioning) - Threads.@threads for ipart in 1:np - if allow_views(p.factorizations[ipart]) - ldiv!(view(u, partitioning[ipart]), facts[ipart], view(v, partitioning[ipart])) - else - uu = u[partitioning[ipart]] - ldiv!(uu, facts[ipart], v[partitioning[ipart]]) - view(u, partitioning[ipart]) .= uu - end - end - return u -end - -Base.eltype(p::BlockPreconditioner) = eltype(p.facts[1]) - - """ BlockPreconBuilder(;precs=UMFPACKPreconBuilder(), partitioning = A -> [1:size(A,1)] @@ -228,120 +78,25 @@ function (blockprecs::BlockPreconBuilder)(A, p) return (bp, LinearAlgebra.I) end -""" - Allow array for precs => different precoms -""" - - -mutable struct _JacobiPreconditioner{Tv} - invdiag::Vector{Tv} -end - -function jacobi(A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti} - invdiag = Array{Tv, 1}(undef, A.n) - n = A.n - @inbounds for i in 1:n - invdiag[i] = one(Tv) / A[i, i] - end - return _JacobiPreconditioner(invdiag) -end - -function jacobi!(p::_JacobiPreconditioner{Tv}, A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti} - n = A.n - @inbounds for i in 1:n - p.invdiag[i] = one(Tv) / A[i, i] - end - return p -end - -mutable struct JacobiPreconditioner - A::AbstractMatrix - factorization::Union{_JacobiPreconditioner, Nothing} - function JacobiPreconditioner(A) - p = new() - p.A = A - p.factorization = nothing - update!(p) - return p - end -end - -function LinearAlgebra.ldiv!(u, p::_JacobiPreconditioner, v) - n = length(p.invdiag) - for i in 1:n - @inbounds u[i] = p.invdiag[i] * v[i] - end - return u -end - -LinearAlgebra.ldiv!(p::_JacobiPreconditioner, v) = ldiv!(v, p, v) - - -""" -``` -JacobiPreconditioner() -JacobiPreconditioner(matrix) -``` - -Jacobi preconditioner. -""" -function JacobiPreconditioner end - -function update!(p::JacobiPreconditioner) - flush!(p.A) - Tv = eltype(p.A) - p.factorization = jacobi(SparseMatrixCSC(p.A)) - return p -end - -LinearAlgebra.ldiv!(u, fact::JacobiPreconditioner, v) = ldiv!(u, fact.factorization, v) -LinearAlgebra.ldiv!(fact::JacobiPreconditioner, v) = ldiv!(fact.factorization, v) - -allow_views(::JacobiPreconditioner) = true -allow_views(::Type{JacobiPreconditioner}) = true """ - struct ProductPreconditioner + JacobiPreconBuilder(; blocksize = 1) -Product of two left preconditioning steps. -The operation ``u=M^{-1}v`` is defined by two simple iteration steps -with two different preconditioners ``M_1`` and ``M_2``: -Let ``u_0=0``. Then calculate -```math - \\begin{align*} - u_1&= u_0 - M_1^{-1}(Au_0 - v) = M_1^{-1}v\\\\ - u &= u_1 - M_2^{-1}(Au_1 - v) - \\end{align*} -``` +Return callable object constructing a left Jacobi preconditioner +to be passed as the `precs` parameter to iterative methods wrapped by LinearSolve.jl. """ -Base.@kwdef struct ProductPreconditioner{TA, TM1, TM2} - A::TA - M1::TM1 - M2::TM2 -end - -function LinearAlgebra.ldiv!(u, p::ProductPreconditioner, v) - (; A, M1, M2) = p - u1 = similar(u) - u2 = similar(u) - ldiv!(u1, M1, v) - mul!(u2, A, u1) - ldiv!(u, M2, v - u2) - u .+= u1 - return u +Base.@kwdef struct JacobiPreconBuilder + blocksize::Int = 1 end -function LinearAlgebra.ldiv!(p::ProductPreconditioner, v) - u = ldiv!(copy(v), p, v) - v .= u - return v -end +(b::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC, p) = (JacobiPreconditioner(A; blocksize = b.blocksize), LinearAlgebra.I) """ - struct ProductPreconBuilder + ProductPreconBuilder(precs1, precs2) -LinearSolve `precs` compatible preonditioner constructor for [`ProductPreconditioner`](@ref) +Return LinearSolve `precs` compatible callable object which constructs +a [`ProductPreconditioner`](@ref) from the `precs1` and `precs2` preconbuilders. """ Base.@kwdef mutable struct ProductPreconBuilder precs1 = JacobiPreconBuilder() diff --git a/src/preconditioners.jl b/src/preconditioners.jl new file mode 100644 index 0000000..34c4c3f --- /dev/null +++ b/src/preconditioners.jl @@ -0,0 +1,241 @@ +# +# This file defines a number of preconditioners uses by the preconbuilders. +# + + +""" + allow_views(::preconditioner) + +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. +""" +allow_views(::Any) = false +allow_views(::ILUZero.ILU0Precon) = true + +################################################################################# +mutable struct JacobiPreconditioner{Tb, BSize} + invdiag::Vector{Tb} +end + +""" + JacobiPreconditioner(A; blocksize=1) + +Create a Jacobi preconditioner. If `blocksize>1` it is a point block preconditioner with blocks of +type of `StaticArrays.SMatrix` +""" +function JacobiPreconditioner(A::AbstractSparseMatrixCSC; blocksize = 1) + n = size(A, 1) + if blocksize == 1 + Tv = eltype(A) + invdiag = Array{Tv, 1}(undef, n) + @inbounds for i in 1:n + invdiag[i] = one(Tv) / A[i, i] + end + return JacobiPreconditioner{Tv, blocksize}(invdiag) + else + n % blocksize == 0 || throw(ArgumentError("JacobiPreconditioner: size(A, 1)=$(n) must be divisible by blocksize=$(blocksize)")) + Tb = SMatrix{blocksize, blocksize, eltype(A), blocksize^2} + nblock = n ÷ blocksize + invdiag = Vector{Tb}(undef, nblock) + block = zeros(eltype(A), blocksize, blocksize) + for iblock in 1:nblock + @inbounds for i in 1:blocksize, j in 1:blocksize + ii = (iblock - 1) * blocksize + i + jj = (iblock - 1) * blocksize + j + block[i, j] = A[ii, jj] + end + invdiag[iblock] = inv(SMatrix{blocksize, blocksize}(block)) + end + return JacobiPreconditioner{Tb, blocksize}(invdiag) + end +end + +function LinearAlgebra.ldiv!(u::AbstractVector{Tu}, p::JacobiPreconditioner{Tb, BSize}, v::AbstractVector{Tv}) where {Tu, Tv, Tb, BSize} + n = length(p.invdiag) + if BSize == 1 + for i in 1:n + @inbounds u[i] = p.invdiag[i] * v[i] + end + else + bu = reinterpret(SVector{BSize, Tu}, u) + bv = reinterpret(SVector{BSize, Tv}, v) + for i in 1:n + @inbounds bu[i] = p.invdiag[i] * bv[i] + end + end + return u +end + +LinearAlgebra.ldiv!(p::JacobiPreconditioner, v) = ldiv!(v, p, v) + +allow_views(::JacobiPreconditioner{T, 1}) where {T} = true +allow_views(::JacobiPreconditioner{T, BSize}) where {T, BSize} = false + +############################################################################ +""" + pointblock(A,blocksize) + +Create a pointblock matrix with entries of type `StaticArrays.SMatrix` of size `blocksize x blocksize` from A. +""" +function pointblock(A0::ExtendableSparseMatrixCSC{Tv, Ti}, blocksize) where {Tv, Ti} + A = SparseMatrixCSC(A0) + colptr = A.colptr + rowval = A.rowval + nzval = A.nzval + n = A.n + block = zeros(Tv, blocksize, blocksize) + nblock = n ÷ blocksize + b = SMatrix{blocksize, blocksize, Tv, blocksize^2}(block) + Tb = typeof(b) + Ab = ExtendableSparseMatrixCSC{Tb, Ti}(nblock, nblock) + for i in 1:n + for k in colptr[i]:(colptr[i + 1] - 1) + j = rowval[k] + iblock = (i - 1) ÷ blocksize + 1 + jblock = (j - 1) ÷ blocksize + 1 + ii = (i - 1) % blocksize + 1 + jj = (j - 1) % blocksize + 1 + block[ii, jj] = nzval[k] + rawupdateindex!(Ab, +, SMatrix{blocksize, blocksize}(block), iblock, jblock) + block[ii, jj] = zero(Tv) + end + end + return flush!(Ab) +end + +""" + ILU0BlockPrecon + +Point-block preconditioner based on ILUZero.ilu0 +""" +struct ILU0BlockPrecon{BSize, BSquare, Tv, Ti} + ilu0::ILUZero.ILU0Precon{SMatrix{BSize, BSize, Tv, BSquare}, Ti, SVector{BSize, Tv}} +end + +function LinearAlgebra.ldiv!( + Y::Vector{Tv}, + A::ILU0BlockPrecon{BSize, BSquare, Tv, Ti}, + B::Vector{Tv} + ) where {BSize, BSquare, Tv, Ti} + BY = reinterpret(SVector{BSize, Tv}, Y) + BB = reinterpret(SVector{BSize, Tv}, B) + ldiv!(BY, A.ilu0, BB) + return Y +end + +####################################################################################### +struct BlockPreconditioner + partitioning::Vector{AbstractVector} + factorizations::Vector{Any} +end + +""" + BlockPreconditioner(A;partitioning, factorizations) + +Create a block preconditioner from partition of unknowns given by `partitioning`, 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]. + +`factorizations` is a thread safe callable `factorizations(A)` (Function or struct) which allows to create a factorization (with `ldiv!` methods) +from a submatrix of A, or a vector thereof. +""" +function BlockPreconditioner( + A::AbstractSparseMatrixCSC; + partitioning = [1:size(A, 1)], + factorizations = (A) -> nothing + ) + nall = sum(length, partitioning) + n = size(A, 1) + if nall != n + @error "BlockPreconditioner: sum(length,partitioning)=$(nall) but n=$(n)" + end + np = length(partitioning) + + if !isa(factorizations, Vector) + factorizations = fill(factorizations, np) + end + + facts = Vector{Any}(undef, np) + Threads.@threads for ipart in 1:np + AP = A[partitioning[ipart], partitioning[ipart]] + facts[ipart] = factorizations[ipart](AP) + end + + return BlockPreconditioner(partitioning, facts) +end + +function LinearAlgebra.ldiv!(p::BlockPreconditioner, v) + (; factorizations, partitioning) = p + np = length(partitioning) + + Threads.@threads for ipart in 1:np + if allow_views(factorizations[ipart]) + ldiv!(factorizations[ipart], view(v, partitioning[ipart])) + else + vv = v[partitioning[ipart]] + ldiv!(factorizations[ipart], vv) + view(v, partitioning[ipart]) .= vv + end + end + return v +end + +function LinearAlgebra.ldiv!(u, p::BlockPreconditioner, v) + (; factorizations, partitioning) = p + np = length(partitioning) + Threads.@threads for ipart in 1:np + if allow_views(factorizations[ipart]) + ldiv!(view(u, partitioning[ipart]), factorizations[ipart], view(v, partitioning[ipart])) + else + uu = u[partitioning[ipart]] + ldiv!(uu, factorizations[ipart], v[partitioning[ipart]]) + view(u, partitioning[ipart]) .= uu + end + end + return u +end + +Base.eltype(p::BlockPreconditioner) = eltype(p.factorizations[1]) + + +####################################################################################### + +""" + ProductPreconditioner(A,M1,M2) + +Product M of two left preconditioning steps using `M1` and `M2`, respectively. + +The operation ``u=M^{-1}v`` is defined by two simple iteration steps +with two different preconditioners ``M_1`` and ``M_2``: + +Let ``u_0=0``. Then calculate +```math + \\begin{align*} + u_1&= u_0 - M_1^{-1}(Au_0 - v) = M_1^{-1}v\\\\ + u &= u_1 - M_2^{-1}(Au_1 - v) + \\end{align*} +``` +""" +Base.@kwdef struct ProductPreconditioner{TA, TM1, TM2} + A::TA + M1::TM1 + M2::TM2 +end + +function LinearAlgebra.ldiv!(u, p::ProductPreconditioner, v) + (; A, M1, M2) = p + u1 = similar(u) + u2 = similar(u) + ldiv!(u1, M1, v) + mul!(u2, A, u1) + ldiv!(u, M2, v - u2) + u .+= u1 + return u +end + +function LinearAlgebra.ldiv!(p::ProductPreconditioner, v) + u = ldiv!(copy(v), p, v) + v .= u + return v +end diff --git a/test/test_block.jl b/test/test_block.jl index c43dbec..dad59b8 100644 --- a/test/test_block.jl +++ b/test/test_block.jl @@ -1,7 +1,7 @@ module test_block using Test using ExtendableSparse -using ExtendableSparse: BlockPreconditioner, jacobi +using ExtendableSparse: BlockPreconditioner, JacobiPreconditioner using ILUZero using IterativeSolvers using LinearAlgebra @@ -10,25 +10,37 @@ using AMGCLWrap ExtendableSparse.allow_views(::typeof(ilu0)) = true -function main(; n = 100) +function main(; n = 100, blocksize = 4) A = fdrand(n, n) - partitioning = [1:2:(n^2), 2:2:(n^2)] - sol0 = ones(n^2) - b = A * ones(n^2) + N = n^2 + + partitioning = [i:blocksize:(n^2) for i in 1:blocksize ] + sol0 = ones(N) + + b = A * ones(N) sol = cg(A, b, Pl = ilu0(A)) @test sol ≈ sol0 - sol = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = ilu0)) + sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorizations = JacobiPreconditioner), log = true) + @test sol ≈ sol0 + + sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorizations = ilu0), log = true) @test sol ≈ sol0 - sol = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = jacobi)) + sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorizations = sparspaklu), log = true) @test sol ≈ sol0 - sol = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = sparspaklu)) + partitioning = [i:(i + blocksize - 1) for i in 1:blocksize:N] + sol, hist_eq = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorizations = sparspaklu), log = true) @test sol ≈ sol0 + sol, hist_pt = cg(A, b, Pl = JacobiPreconditioner(A; blocksize), log = true) + @test sol ≈ sol0 + @test hist_pt.iters == hist_eq.iters + + return end diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl index 8756567..c7a5faf 100644 --- a/test/test_linearsolve.jl +++ b/test/test_linearsolve.jl @@ -83,7 +83,8 @@ allprecs = [ ExtendableSparse.ILUZeroPreconBuilder(), ExtendableSparse.ILUZeroPreconBuilder(; blocksize = 2), ExtendableSparse.ILUTPreconBuilder(), - # ExtendableSparse.JacobiPreconBuilder(), + ExtendableSparse.JacobiPreconBuilder(), + ExtendableSparse.JacobiPreconBuilder(; blocksize = 2), SmoothedAggregationPreconBuilder(), RugeStubenPreconBuilder(), ] @@ -107,8 +108,10 @@ luprecs = [ExtendableSparse.LinearSolvePreconBuilder(factorization) for factori @testset "block preconditioning" begin n = 100 + blocksize = 4 A = fdrand(n, n) - partitioning = A -> [1:2:size(A, 1), 2:2:size(A, 1)] + N = n^2 + partitioning = [i:blocksize:(n^2) for i in 1:blocksize ] sol0 = ones(n^2) b = A * ones(n^2)