Skip to content

Remove BlockArrays dep and handle blocked arrays manually #111

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pjaap merged 3 commits into from
Feb 19, 2026
Merged

Remove BlockArrays dep and handle blocked arrays manually #111

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0618e9d
Remove BlockArrays dep and handle blocked arrays manually
pjaap Feb 18, 2026
f826f08
Add more verbosity > 1 info about assembly/solution process
pjaap Feb 18, 2026
60e8563
fix wrong multiplier in mul!
pjaap Feb 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 0 additions & 2 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@ authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.ja

[deps]
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
Expand All @@ -30,7 +29,6 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
[compat]
ADTypes = "1.16.0"
Aqua = "0.8"
BlockArrays = "1.7.0"
ChunkSplitters = "3.1.2"
CommonSolve = "0.2"
DiffResults = "1"
Expand Down
1 change: 0 additions & 1 deletion src/ExtendableFEM.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@ module ExtendableFEM

using ADTypes: ADTypes, KnownJacobianSparsityDetector
using ChunkSplitters: chunks
using BlockArrays: BlockMatrix, BlockVector, Block, blocks, axes, mortar
using CommonSolve: CommonSolve
using DifferentiationInterface: DifferentiationInterface,
AutoSparse, AutoForwardDiff, prepare_jacobian
Expand Down
67 changes: 24 additions & 43 deletions src/solvers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@ function assemble_system!(A, b, sol, PD, SC, timer; kwargs...)
if SC.parameters[:initialized]
for op in PD.operators
@timeit timer "$(op.parameters[:name])" begin
SC.parameters[:verbosity] > 1 && @info "assemble $(op.parameters[:name]) ... "
stats = @timed assemble!(
A, b, sol, op, SC;
time = SC.parameters[:time],
Expand All @@ -152,6 +153,7 @@ function assemble_system!(A, b, sol, PD, SC, timer; kwargs...)
else
for op in PD.operators
@timeit timer "$(op.parameters[:name]) (first)" begin
SC.parameters[:verbosity] > 1 && @info "first assembly of $(op.parameters[:name]) ... "
stats = @timed assemble!(
A, b, sol, op, SC;
time = SC.parameters[:time],
Expand All @@ -163,14 +165,18 @@ function assemble_system!(A, b, sol, PD, SC, timer; kwargs...)
end
## assemble restrictions
for restriction in PD.restrictions
SC.parameters[:verbosity] > 1 && @info "assemble $(restriction.parameters[:name]) ... "
@timeit timer "$(restriction.parameters[:name])" assemble!(restriction, sol, SC; kwargs...)
end
end

SC.parameters[:verbosity] > 1 && @info "flush FEMatrix... "
flush!(A.entries)

# Apply penalties
for op in PD.operators
@timeit timer "$(op.parameters[:name]) (penalties)" begin
SC.parameters[:verbosity] > 1 && @info "apply penalties of $(op.parameters[:name]) ... "
stats = @timed apply_penalties!(
A, b, sol, op, SC;
assemble_matrix = !SC.parameters[:initialized] || !SC.parameters[:constant_matrix],
Expand All @@ -181,6 +187,8 @@ function assemble_system!(A, b, sol, PD, SC, timer; kwargs...)
time_assembly += stats.time
allocs_assembly += stats.bytes
end

SC.parameters[:verbosity] > 1 && @info "flush FEMatrix... "
flush!(A.entries)

return time_assembly, allocs_assembly
Expand Down Expand Up @@ -271,7 +279,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
end

## restrictions only involve the blocks coressponding to the unknowns and not necessarily the full sol.entries
sol_freedofs = mortar([view(sol[u]) for u in unknowns])
sol_freedofs = vcat([view(sol[u]) for u in unknowns]...)

if length(PD.restrictions) == 0
if linsolve_needs_matrix
Expand All @@ -285,7 +293,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,

# block sizes for the block matrix
block_sizes = [length(b_unrestricted), [ length(b) for b in restriction_rhss ]...]
total_size = sum(block_sizes)
block_ends = cumsum(block_sizes)

# total number of additional LM dofs
nLMs = @views sum(block_sizes[2:end])
Expand Down Expand Up @@ -315,14 +323,15 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
# we need the inverse column permutation
ipcol = invperm(pcol)

# the new combined restriction block (add another transpose)
combined_restriction_matrix = R[:, ipcol]'
# the new combined restriction block (add another transpose and convert into CSC matrix)
combined_restriction_matrix = sparse(R[:, ipcol]')

# the new combined restriction rhs
combined_restriction_rhs = Q'combined_restriction_rhs[prow]

# compress the column space
qr_rank = rank(qr_result)
SC.parameters[:verbosity] > 1 && @info "QR decomposition has rank $qr_rank"
@assert norm(combined_restriction_rhs[(qr_rank + 1):end]) ≤ 1.0e-12 * norm(combined_restriction_rhs) "the rhs of the restriction is not in the image"

combined_restriction_matrix = combined_restriction_matrix[:, 1:qr_rank]
Expand All @@ -337,63 +346,36 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,

# update sizes
block_sizes = [length(b_unrestricted), length(combined_restriction_rhs)]
total_size = sum(block_sizes)
block_ends = cumsum(block_sizes)

SC.parameters[:verbosity] > 1 && @info "Created new CompressedRestriction"
end

A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes)
A_block[Block(1, 1)] = A_unrestricted
linsolve_A = spzeros(Tv, block_ends[end], block_ends[end])
linsolve_A[1:block_ends[1], 1:block_ends[1]] = A_unrestricted
end

## we need to add the (initial) solution to the rhs, since we work with the residual equation
for (B, rhs) in zip(restriction_matrices, restriction_rhss)
rhs .-= B'sol_freedofs
# rhs -= B'sol_freedofs
mul!(rhs, B', sol_freedofs, -1.0, 1.0)
end

b_block = BlockVector(zeros(Tv, total_size), block_sizes)
b_block[Block(1)] = b_unrestricted
linsolve_b = zeros(Tv, block_ends[end])
linsolve_b[1:block_ends[1]] = b_unrestricted

for i in eachindex(restriction_matrices)
if linsolve_needs_matrix
A_block[Block(1, i + 1)] = restriction_matrices[i]
A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i])
linsolve_A[1:block_ends[1], (block_ends[i] + 1):block_ends[i + 1]] = restriction_matrices[i]
linsolve_A[(block_ends[i] + 1):block_ends[i + 1], 1:block_ends[1]] = restriction_matrices[i]'
end
b_block[Block(i + 1)] = restriction_rhss[i]
linsolve_b[(block_ends[i] + 1):block_ends[i + 1]] = restriction_rhss[i]

end

# convert to dense vectors
linsolve_b = Vector(b_block)

if linsolve_needs_matrix

# linsolve.A = sparse(A_block) # convert to CSC Matrix is very slow https://github.com/JuliaArrays/BlockArrays.jl/issues/78
# do it manually:
A_flat = spzeros(size(A_block))

lasts_row = [0, axes(A_block)[1].lasts...] # add leading zero
lasts_col = [0, axes(A_block)[2].lasts...] # add leading zero

(n_row, n_col) = size(blocks(A_block))


SC.parameters[:verbosity] > 1 && @info "Assemble flat matrix out of blocks..."
for i in 1:n_row, j in 1:n_col
range_row = (lasts_row[i] + 1):lasts_row[i + 1]
range_col = (lasts_col[j] + 1):lasts_col[j + 1]

# write each block directly in the resulting matrix
A_flat[range_row, range_col] = A_block[Block(i, j)]
end
SC.parameters[:verbosity] > 1 && @info "... matrix assembly complete!"

linsolve_A = A_flat
end
SC.parameters[:verbosity] > 1 && @info "done setting blocked matrix/vector"
end
end


if SC.parameters[:initialized]
# set/update the linear system
SC.linsolver.b = linsolve_b
Expand Down Expand Up @@ -442,7 +424,6 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,

if length(PD.restrictions) > 0
# extract all residuals for the restriction blocks
block_ends = cumsum(block_sizes)
restriction_residuals = [norm(residual_flat[(block_ends[i] + 1):block_ends[i + 1]]) for i in 1:(length(block_sizes) - 1) ]
push!(stats[:restriction_residuals], restriction_residuals)

Expand Down
Loading