From 045557086f080d0b0edd4e3bd1b6e8027e87579f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Fri, 10 Jul 2026 21:31:29 +0200
Subject: [PATCH 1/7] point block jacobi

---
 src/preconbuilders.jl    | 123 ++++++++++++++++++---------------------
 test/test_block.jl       |  28 ++++++---
 test/test_linearsolve.jl |   7 ++-
 3 files changed, 83 insertions(+), 75 deletions(-)

diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl
index 79dbb86..f201051 100644
--- a/src/preconbuilders.jl
+++ b/src/preconbuilders.jl
@@ -10,16 +10,6 @@ end
 (::LinearSolvePreconBuilder)(A, p) = error("import LinearSolve in order to use LinearSolvePreconBuilder")
 
 
-"""
-    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)
 
@@ -58,7 +48,7 @@ function pointblock(A0::ExtendableSparseMatrixCSC{Tv, Ti}, blocksize) where {Tv,
     n = A.n
     block = zeros(Tv, blocksize, blocksize)
     nblock = n ÷ blocksize
-    b = SMatrix{blocksize, blocksize}(block)
+    b = SMatrix{blocksize, blocksize, Tv, blocksize^2}(block)
     Tb = typeof(b)
     Ab = ExtendableSparseMatrixCSC{Tb, Ti}(nblock, nblock)
 
@@ -154,10 +144,11 @@ function update!(precon::BlockPreconditioner)
 
     np = length(precon.partitioning)
     precon.facts = Vector{Any}(undef, np)
-    return Threads.@threads for ipart in 1:np
+    Threads.@threads for ipart in 1:np
         AP = precon.A[precon.partitioning[ipart], precon.partitioning[ipart]]
         precon.facts[ipart] = precon.factorizations[ipart](AP)
     end
+    return
 end
 
 
@@ -182,8 +173,8 @@ 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])
+    Threads.@threads  for ipart in 1:np
+        if allow_views(p.facts[ipart])
             ldiv!(view(u, partitioning[ipart]), facts[ipart], view(v, partitioning[ipart]))
         else
             uu = u[partitioning[ipart]]
@@ -228,77 +219,79 @@ function (blockprecs::BlockPreconBuilder)(A, p)
     return (bp, LinearAlgebra.I)
 end
 
-"""
-    Allow array for precs => different precoms
-"""
-
 
-mutable struct _JacobiPreconditioner{Tv}
-    invdiag::Vector{Tv}
+mutable struct JacobiPreconditioner{Tb, N}
+    invdiag::Vector{Tb}
 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]
+"""
+    JacobiPreconditioner(A; blocksize=1)
+"""
+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
+        Tb = SMatrix{blocksize, blocksize, eltype(A), blocksize^2}
+        nblock = n ÷ blocksize
+        invdiag = Array{Tb, 1}(undef, nblock)
+        block = zeros(eltype(A), blocksize, blocksize)
+        for iblock in 1:nblock
+            for i in 1:blocksize
+                for j in 1:blocksize
+                    ii = (iblock - 1) * blocksize + i
+                    jj = (iblock - 1) * blocksize + j
+                    block[i, j] = A[ii, jj]
+                    sblock = SMatrix{blocksize, blocksize}(block)
+                    invdiag[iblock] = inv(sblock)
+                end
+            end
+        end
+        return JacobiPreconditioner{Tb, blocksize}(invdiag)
     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)
+function LinearAlgebra.ldiv!(u::AbstractVector{Tu}, p::JacobiPreconditioner{Tb, N}, v::AbstractVector{Tv}) where {Tu, Tv, Tb, N}
     n = length(p.invdiag)
-    for i in 1:n
-        @inbounds u[i] = p.invdiag[i] * v[i]
+    if N == 1
+        for i in 1:n
+            @inbounds u[i] = p.invdiag[i] * v[i]
+        end
+    else
+        bu = reinterpret(SVector{N, Tu}, u)
+        bv = reinterpret(SVector{N, 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)
+LinearAlgebra.ldiv!(p::JacobiPreconditioner, v) = ldiv!(v, p, v)
 
+allow_views(::JacobiPreconditioner{T, 1}) where {T} = true
+allow_views(::JacobiPreconditioner{T, N}) where {T, N} = true
 
-"""
-```
-JacobiPreconditioner()
-JacobiPreconditioner(matrix)
-```
 
-Jacobi preconditioner.
 """
-function JacobiPreconditioner end
+    JacobiPreconBuilder()
 
-function update!(p::JacobiPreconditioner)
-    flush!(p.A)
-    Tv = eltype(p.A)
-    p.factorization = jacobi(SparseMatrixCSC(p.A))
-    return p
+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 JacobiPreconBuilder
+    blocksize::Int = 1
 end
+(b::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC, p) = (JacobiPreconditioner(A; blocksize = b.blocksize), LinearAlgebra.I)
 
-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
diff --git a/test/test_block.jl b/test/test_block.jl
index c43dbec..f443683 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, factorization = JacobiPreconditioner), log = true)
+    @test sol ≈ sol0
+
+    sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = 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, factorization = 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, factorization = 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)
 

From 7ac9a7990f3946c1366c257e79300539015c48ee Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sat, 11 Jul 2026 22:12:40 +0200
Subject: [PATCH 2/7] reorg: separate preconditioners and preconbuilders

---
 src/ExtendableSparse.jl |   1 +
 src/preconbuilders.jl   | 249 +---------------------------------------
 src/preconditioners.jl  | 244 +++++++++++++++++++++++++++++++++++++++
 3 files changed, 247 insertions(+), 247 deletions(-)
 create mode 100644 src/preconditioners.jl

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 f201051..c06a757 100644
--- a/src/preconbuilders.jl
+++ b/src/preconbuilders.jl
@@ -20,64 +20,15 @@ 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, 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
-
 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)
 
@@ -90,104 +41,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)
-    Threads.@threads for ipart in 1:np
-        AP = precon.A[precon.partitioning[ipart], precon.partitioning[ipart]]
-        precon.facts[ipart] = precon.factorizations[ipart](AP)
-    end
-    return
-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.facts[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)]
@@ -220,67 +73,6 @@ function (blockprecs::BlockPreconBuilder)(A, p)
 end
 
 
-mutable struct JacobiPreconditioner{Tb, N}
-    invdiag::Vector{Tb}
-end
-
-
-"""
-    JacobiPreconditioner(A; blocksize=1)
-"""
-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
-        Tb = SMatrix{blocksize, blocksize, eltype(A), blocksize^2}
-        nblock = n ÷ blocksize
-        invdiag = Array{Tb, 1}(undef, nblock)
-        block = zeros(eltype(A), blocksize, blocksize)
-        for iblock in 1:nblock
-            for i in 1:blocksize
-                for j in 1:blocksize
-                    ii = (iblock - 1) * blocksize + i
-                    jj = (iblock - 1) * blocksize + j
-                    block[i, j] = A[ii, jj]
-                    sblock = SMatrix{blocksize, blocksize}(block)
-                    invdiag[iblock] = inv(sblock)
-                end
-            end
-        end
-        return JacobiPreconditioner{Tb, blocksize}(invdiag)
-    end
-end
-
-
-function LinearAlgebra.ldiv!(u::AbstractVector{Tu}, p::JacobiPreconditioner{Tb, N}, v::AbstractVector{Tv}) where {Tu, Tv, Tb, N}
-    n = length(p.invdiag)
-    if N == 1
-        for i in 1:n
-            @inbounds u[i] = p.invdiag[i] * v[i]
-        end
-    else
-        bu = reinterpret(SVector{N, Tu}, u)
-        bv = reinterpret(SVector{N, 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, N}) where {T, N} = true
-
-
 """
     JacobiPreconBuilder()
 
@@ -290,45 +82,8 @@ to be passed as the `precs` parameter to iterative methods wrapped by LinearSolv
 Base.@kwdef struct JacobiPreconBuilder
     blocksize::Int = 1
 end
-(b::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC, p) = (JacobiPreconditioner(A; blocksize = b.blocksize), LinearAlgebra.I)
 
-
-"""
-    struct ProductPreconditioner
-
-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*}
-```
-"""
-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
+(b::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC, p) = (JacobiPreconditioner(A; blocksize = b.blocksize), LinearAlgebra.I)
 
 
 """
diff --git a/src/preconditioners.jl b/src/preconditioners.jl
new file mode 100644
index 0000000..5e110a4
--- /dev/null
+++ b/src/preconditioners.jl
@@ -0,0 +1,244 @@
+"""
+    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
+
+mutable struct JacobiPreconditioner{Tb, N}
+    invdiag::Vector{Tb}
+end
+
+
+"""
+    JacobiPreconditioner(A; blocksize=1)
+"""
+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
+        Tb = SMatrix{blocksize, blocksize, eltype(A), blocksize^2}
+        nblock = n ÷ blocksize
+        invdiag = Array{Tb, 1}(undef, nblock)
+        block = zeros(eltype(A), blocksize, blocksize)
+        for iblock in 1:nblock
+            for i in 1:blocksize
+                for j in 1:blocksize
+                    ii = (iblock - 1) * blocksize + i
+                    jj = (iblock - 1) * blocksize + j
+                    block[i, j] = A[ii, jj]
+                    sblock = SMatrix{blocksize, blocksize}(block)
+                    invdiag[iblock] = inv(sblock)
+                end
+            end
+        end
+        return JacobiPreconditioner{Tb, blocksize}(invdiag)
+    end
+end
+
+
+function LinearAlgebra.ldiv!(u::AbstractVector{Tu}, p::JacobiPreconditioner{Tb, N}, v::AbstractVector{Tv}) where {Tu, Tv, Tb, N}
+    n = length(p.invdiag)
+    if N == 1
+        for i in 1:n
+            @inbounds u[i] = p.invdiag[i] * v[i]
+        end
+    else
+        bu = reinterpret(SVector{N, Tu}, u)
+        bv = reinterpret(SVector{N, 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, N}) where {T, N} = true
+
+
+"""
+    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, 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
+
+struct ILU0BlockPrecon{N, NN, Tv, Ti}
+    ilu0::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}}
+end
+
+function LinearAlgebra.ldiv!(
+        Y::Vector{Tv},
+        A::ILU0BlockPrecon{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
+
+######
+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
+
+
+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)
+    Threads.@threads for ipart in 1:np
+        AP = precon.A[precon.partitioning[ipart], precon.partitioning[ipart]]
+        precon.facts[ipart] = precon.factorizations[ipart](AP)
+    end
+    return
+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.facts[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])
+
+
+"""
+    struct ProductPreconditioner
+
+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*}
+```
+"""
+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

From 7708a9fea6f6b1a412fa1f9412aaea030d26c2a7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sat, 11 Jul 2026 23:15:50 +0200
Subject: [PATCH 3/7] Rework BlockPreconditioner

The way it was structured still used old ideas from ExtendableSparse
preconditioning. Without breaking (due to being not exported anymore),
we can implement this in a much clearer and simple way.
---
 src/preconbuilders.jl  |  17 +++--
 src/preconditioners.jl | 144 ++++++++++++++++++++---------------------
 test/test_block.jl     |   8 +--
 3 files changed, 87 insertions(+), 82 deletions(-)

diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl
index c06a757..09945e9 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 iterativer solvers as they are handeled in LinearSolve.jl.
+# See https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#Specifying-Preconditioners
+#
+
 """
         LinearSolvePreconBuilder(; method=UMFPACKFactorization())
 
@@ -11,7 +17,7 @@ end
 
 
 """
-    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,7 +36,7 @@ function (b::ILUZeroPreconBuilder)(A0, p)
 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)
@@ -74,7 +80,7 @@ end
 
 
 """
-    JacobiPreconBuilder()
+    JacobiPreconBuilder(; blocksize = 1)
 
 Return callable object constructing a left Jacobi preconditioner
 to be passed as the `precs` parameter to iterative methods wrapped by LinearSolve.jl.
@@ -87,9 +93,10 @@ end
 
 
 """
-    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
index 5e110a4..3aaaf46 100644
--- a/src/preconditioners.jl
+++ b/src/preconditioners.jl
@@ -1,23 +1,31 @@
+#
+# 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.
+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, N}
+#################################################################################
+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)
@@ -36,8 +44,7 @@ function JacobiPreconditioner(A::AbstractSparseMatrixCSC; blocksize = 1)
                     ii = (iblock - 1) * blocksize + i
                     jj = (iblock - 1) * blocksize + j
                     block[i, j] = A[ii, jj]
-                    sblock = SMatrix{blocksize, blocksize}(block)
-                    invdiag[iblock] = inv(sblock)
+                    invdiag[iblock] = inv(SMatrix{blocksize, blocksize}(block))
                 end
             end
         end
@@ -45,16 +52,15 @@ function JacobiPreconditioner(A::AbstractSparseMatrixCSC; blocksize = 1)
     end
 end
 
-
-function LinearAlgebra.ldiv!(u::AbstractVector{Tu}, p::JacobiPreconditioner{Tb, N}, v::AbstractVector{Tv}) where {Tu, Tv, Tb, N}
+function LinearAlgebra.ldiv!(u::AbstractVector{Tu}, p::JacobiPreconditioner{Tb, BSize}, v::AbstractVector{Tv}) where {Tu, Tv, Tb, BSize}
     n = length(p.invdiag)
-    if N == 1
+    if BSize == 1
         for i in 1:n
             @inbounds u[i] = p.invdiag[i] * v[i]
         end
     else
-        bu = reinterpret(SVector{N, Tu}, u)
-        bv = reinterpret(SVector{N, Tv}, v)
+        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
@@ -65,13 +71,13 @@ end
 LinearAlgebra.ldiv!(p::JacobiPreconditioner, v) = ldiv!(v, p, v)
 
 allow_views(::JacobiPreconditioner{T, 1}) where {T} = true
-allow_views(::JacobiPreconditioner{T, N}) where {T, N} = true
-
+allow_views(::JacobiPreconditioner{T, BSize}) where {T, BSize} = false
 
+############################################################################
 """
-    pointblock(matrix,blocksize)
+        pointblock(A,blocksize)
 
-Create a pointblock matrix.
+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)
@@ -84,8 +90,6 @@ function pointblock(A0::ExtendableSparseMatrixCSC{Tv, Ti}, blocksize) where {Tv,
     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]
@@ -101,86 +105,77 @@ function pointblock(A0::ExtendableSparseMatrixCSC{Tv, Ti}, blocksize) where {Tv,
     return flush!(Ab)
 end
 
-struct ILU0BlockPrecon{N, NN, Tv, Ti}
-    ilu0::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}}
+"""
+    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{N, NN, Tv, Ti},
+        A::ILU0BlockPrecon{BSize, BSquare, Tv, Ti},
         B::Vector{Tv}
-    ) where {N, NN, Tv, Ti}
-    BY = reinterpret(SVector{N, Tv}, Y)
-    BB = reinterpret(SVector{N, Tv}, B)
+    ) 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
 
-######
-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
+#######################################################################################
+struct BlockPreconditioner
+    partitioning::Vector{AbstractVector}
+    factorizations::Vector{Any}
 end
 
-
 """
-     BlockPreconditioner(;partitioning, factorization)
+     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].
-Factorization is a callable (Function or struct) which allows to create a factorization (with `ldiv!` methods) from a submatrix of A.
-"""
-function BlockPreconditioner end
 
-
-function update!(precon::BlockPreconditioner)
-    flush!(precon.A)
-    nall = sum(length, precon.partitioning)
-    n = size(precon.A, 1)
+`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
-        @warn "sum(length,partitioning)=$(nall) but n=$(n)"
+        @error "BlockPreconditioner: sum(length,partitioning)=$(nall) but n=$(n)"
     end
+    np = length(partitioning)
 
-    if isnothing(precon.partitioning)
-        partitioning = [1:n]
+    if !isa(factorizations, Vector)
+        factorizations = fill(factorizations, np)
     end
 
-    np = length(precon.partitioning)
-    precon.facts = Vector{Any}(undef, np)
+    facts = Vector{Any}(undef, np)
     Threads.@threads for ipart in 1:np
-        AP = precon.A[precon.partitioning[ipart], precon.partitioning[ipart]]
-        precon.facts[ipart] = precon.factorizations[ipart](AP)
+        AP = A[partitioning[ipart], partitioning[ipart]]
+        facts[ipart] = factorizations[ipart](AP)
     end
-    return
-end
 
+    return BlockPreconditioner(partitioning, facts)
+end
 
 function LinearAlgebra.ldiv!(p::BlockPreconditioner, v)
-    partitioning = p.partitioning
-    facts = p.facts
+    (; factorizations, partitioning) = p
     np = length(partitioning)
 
     Threads.@threads for ipart in 1:np
-        if allow_views(p.factorizations[ipart])
-            ldiv!(facts[ipart], view(v, partitioning[ipart]))
+        if allow_views(factorizations[ipart])
+            ldiv!(factorizations[ipart], view(v, partitioning[ipart]))
         else
             vv = v[partitioning[ipart]]
-            ldiv!(facts[ipart], vv)
+            ldiv!(factorizations[ipart], vv)
             view(v, partitioning[ipart]) .= vv
         end
     end
@@ -188,30 +183,33 @@ function LinearAlgebra.ldiv!(p::BlockPreconditioner, v)
 end
 
 function LinearAlgebra.ldiv!(u, p::BlockPreconditioner, v)
-    partitioning = p.partitioning
-    facts = p.facts
+    (; factorizations, partitioning) = p
     np = length(partitioning)
     Threads.@threads  for ipart in 1:np
-        if allow_views(p.facts[ipart])
-            ldiv!(view(u, partitioning[ipart]), facts[ipart], view(v, partitioning[ipart]))
+        if allow_views(factorizations[ipart])
+            ldiv!(view(u, partitioning[ipart]), factorizations[ipart], view(v, partitioning[ipart]))
         else
             uu = u[partitioning[ipart]]
-            ldiv!(uu, facts[ipart], v[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.facts[1])
+Base.eltype(p::BlockPreconditioner) = eltype(p.factorizations[1])
 
 
+#######################################################################################
+
 """
-    struct ProductPreconditioner
+    ProductPreconditioner(A,M1,M2)
+
+Product M of two left preconditioning steps using `M1` and `M2`, respectively.
 
-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*}
diff --git a/test/test_block.jl b/test/test_block.jl
index f443683..dad59b8 100644
--- a/test/test_block.jl
+++ b/test/test_block.jl
@@ -23,17 +23,17 @@ function main(; n = 100, blocksize = 4)
 
     @test sol ≈ sol0
 
-    sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = JacobiPreconditioner), log = true)
+    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, factorization = ilu0), log = true)
+    sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorizations = ilu0), log = true)
     @test sol ≈ sol0
 
-    sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = sparspaklu), log = true)
+    sol, hist = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorizations = sparspaklu), log = true)
     @test sol ≈ sol0
 
     partitioning = [i:(i + blocksize - 1) for i in 1:blocksize:N]
-    sol, hist_eq = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = sparspaklu), log = true)
+    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)

From e2aec088b0f8d6e606c6d6b5abed1e3c4aba1aa9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sat, 11 Jul 2026 23:56:18 +0200
Subject: [PATCH 4/7] update project+changelog

---
 CHANGELOG.md | 6 +++++-
 Project.toml | 2 +-
 2 files changed, 6 insertions(+), 2 deletions(-)

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 <juergen.fuhrmann@wias-berlin.de>", "Daniel Runge"]
 
 [deps]

From 6b626b837a238accfa8a590503f589d98f864d14 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sun, 12 Jul 2026 00:09:39 +0200
Subject: [PATCH 5/7] take `inv` out of the block filling loop for Jacobi

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
---
 src/preconditioners.jl | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/src/preconditioners.jl b/src/preconditioners.jl
index 3aaaf46..34c4c3f 100644
--- a/src/preconditioners.jl
+++ b/src/preconditioners.jl
@@ -34,19 +34,18 @@ function JacobiPreconditioner(A::AbstractSparseMatrixCSC; blocksize = 1)
         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 = Array{Tb, 1}(undef, nblock)
+        invdiag = Vector{Tb}(undef, nblock)
         block = zeros(eltype(A), blocksize, blocksize)
         for iblock in 1:nblock
-            for i in 1:blocksize
-                for j in 1:blocksize
-                    ii = (iblock - 1) * blocksize + i
-                    jj = (iblock - 1) * blocksize + j
-                    block[i, j] = A[ii, jj]
-                    invdiag[iblock] = inv(SMatrix{blocksize, blocksize}(block))
-                end
+            @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

From 845e841e8b24f82a340dc6d0a01f8630548d683e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sun, 12 Jul 2026 00:11:05 +0200
Subject: [PATCH 6/7] typos

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
---
 src/preconbuilders.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl
index 09945e9..71854a4 100644
--- a/src/preconbuilders.jl
+++ b/src/preconbuilders.jl
@@ -1,6 +1,6 @@
 #
 # This file defines `preconbuilders` which allow to specify `precs` preconditioner constructors
-# for iterativer solvers as they are handeled in LinearSolve.jl.
+# for iterative solvers as they are handled in LinearSolve.jl.
 # See https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#Specifying-Preconditioners
 #
 

From 9658f11944841e656ac3f8874fab335d84889242 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sun, 12 Jul 2026 00:25:31 +0200
Subject: [PATCH 7/7] update .codespellrc

---
 .codespellrc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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