SymmetricBlockMatrix

A SymmetricBlockMatrix is a matrices that is symmetric and (at least) sparse at the block level. This can, for example, occur when the matrix originates from a fast multipole method (FMM) or any other algorithm that naturally groups rows and columns into interacting blocks (e.g. near‑field interactions).

Below you will find a compact, step‑by‑step guide that shows how to

1. Build a symmetric block‑sparse matrix from two lists of dense blocks and their index maps.
2. Multiply it with vectors (including transposes).
3. Convert it to a standard SparseMatrixCSC for comparison.
4. Enable multi‑threaded construction with OhMyThreads.

Constructing a SymmetricBlockMatrix

using CompScienceMeshes, BEAST, H2Trees
using UnicodePlots
using BlockSparseMatrices

# make sure that only one of the two symmetric  blocks is stored
struct GalerkinSymmetricIsNearFunctor{N}
    isnear::N
end

function (f::GalerkinSymmetricIsNearFunctor)(tree, nodea, nodeb)
    if H2Trees.isleaf(tree, nodeb)
        return f.isnear(tree, nodea, nodeb) && nodea > nodeb
    else
        return f.isnear(tree, nodea, nodeb)
    end
end

m = meshsphere(1.0, 0.1)
X = raviartthomas(m)
sizeS = (numfunctions(X), numfunctions(X))
tree = TwoNTree(X, 0.2)

selfvalues, colvalues, rowvalues = H2Trees.nearinteractions(
    tree; extractselfvalues=true, isnear=GalerkinSymmetricIsNearFunctor(H2Trees.isnear)
)

diagonals = Matrix{Float64}[]
for i in eachindex(selfvalues)
    push!(diagonals, rand(Float64, length(selfvalues[i]), length(selfvalues[i])))
    diagonals[end] = (diagonals[end] + transpose(diagonals[end])) / 2  # diagonals are symmetric
end

offdiagonals = Matrix{Float64}[]
for i in eachindex(colvalues)
    push!(offdiagonals, rand(Float64, length(colvalues[i]), length(rowvalues[i])))
end

S = SymmetricBlockMatrix(diagonals, selfvalues, offdiagonals, colvalues, rowvalues, sizeS)

4821×4821 SymmetricBlockMatrix{Float64} with 3050335 non-zero entries
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⣿⣿⣿⣿⣿⡇⡇⠀⠀⡇⠃⠀⠀⠀⠀⠀⣿⣿⡇⠀⠀⢰⣦⢠⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⡄⡄⣤⣶⣶⣦⠀
⠀⣿⣿⣿⣿⣿⣷⣷⠀⠀⣷⣦⠀⠀⠀⠀⠀⣤⣤⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠶⠗⠛⠻⠿⠿⠟⠀
⠀⠿⠿⢿⣿⣿⣿⣿⣤⣤⣿⡏⠀⠀⠀⠀⠀⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⢀⣀⡀⠀⢀⡀⠀⠀⠀
⠀⠉⠉⠙⠛⠛⣿⣿⣿⣾⣿⣆⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠸⠿⠽⠿⠯⠍⠘⠋⠉⠁⠀
⠀⠤⠤⢤⣤⣤⣿⣾⣿⣿⣿⣿⣿⣿⠀⠀⠀⢤⣤⠄⠀⠀⠀⠀⠀⠀⠀⠀⠤⠿⠇⠸⠿⠼⠿⠦⠄⠀⠀⠀⠀⠀
⠀⠉⠀⠈⠛⠋⠉⠈⢹⣿⣿⣿⣿⣿⣿⣿⡟⠛⠛⠛⠓⠂⠀⠀⠀⠀⢄⣤⡤⣤⡄⢠⣤⠤⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠘⠛⠛⣿⣿⣿⣿⣿⣧⢀⣠⠀⣤⣤⠀⢀⠀⣀⣴⣶⣦⣶⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣿⠿⠿⣿⣿⣿⣿⣿⣿⣿⡏⠉⠈⠀⠉⠉⠉⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⣿⣿⠀⣿⡇⠀⠀⠀⠀⣷⣿⠀⠀⣰⣿⣿⣿⣿⣿⣿⣁⣀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡀⠀
⠀⠉⠉⠀⠉⠁⠀⠀⠀⠀⠁⢿⠀⠀⣤⣿⣿⣿⣿⣿⣿⣿⣿⡿⠸⠿⢼⣤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢹⠁⠀
⠀⢀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠛⡏⠉⠁⢸⣿⣿⣿⣿⣾⣾⣿⣿⠛⠀⠀⡗⠂⠀⠀⠀⠀⠀⡀⠀⣰⣀⠀
⠀⠈⣛⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠂⠀⠀⠈⣛⡋⣺⣿⣻⣾⣿⣶⠒⠀⠀⡶⡶⠀⠀⠀⠀⢰⣶⠀⣻⣿⠀
⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⢄⢀⣼⡇⠀⠀⠀⣛⣇⣿⣿⢻⣿⣿⣿⣿⣤⣤⣿⡇⠀⠀⠀⠀⠈⠉⠀⠙⠉⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡄⠀⡿⠸⣿⡇⠀⠀⠀⠀⠛⠛⠀⠘⠀⠛⣿⣿⣿⣿⣿⣧⣤⣤⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠿⠇⠀⠿⠸⠿⠇⠀⠀⠀⠀⠀⢤⠤⢠⡤⣤⣿⣿⣿⣿⣿⣿⢿⣿⣀⣀⣠⣤⡀⠠⠤⠀
⠀⠀⠀⠐⠀⠀⠠⣶⡆⣶⡆⠀⣶⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠘⠋⠉⠉⠉⣿⣿⣟⣻⣾⣿⣿⣿⡟⠛⠋⠉⠉⠀
⠀⠀⠀⢠⡄⠀⢰⣷⡇⣶⡇⠀⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠛⠛⢻⣿⣿⣿⣿⣿⣧⣤⡄⠀⠀⠀
⠀⠀⠭⣽⠁⠀⠈⡏⠇⠈⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡀⠀⠀⠀⠀⣸⣿⠿⠿⣿⣿⣿⣿⣯⣭⣭⠀
⠀⢠⣿⣿⡆⠀⠰⡶⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠘⠛⠃⠀⠀⠀⠀⠻⡿⠀⠀⠿⡿⣿⣿⣿⣿⣿⠀
⠀⠸⣿⣿⠇⠀⠀⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠸⠗⠒⠐⢺⣿⣾⡗⠀⠀⠀⠀⡆⡇⠀⠀⠀⡇⣿⣿⣿⣿⣿⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀

S now behaves like a regular matrix: you can query its size, inspect its blocks, etc.

Matrix–Vector Products

y = randn(Float64, numfunctions(X))
@time S * y
@time S' * y
@time transpose(S) * y
@show maximum(abs, S * y - transpose(S) * y)

  0.394126 seconds (1.36 M allocations: 72.425 MiB, 99.41% compilation time)
  0.248337 seconds (702.90 k allocations: 36.145 MiB, 5.96% gc time, 99.07% compilation time)
  0.007767 seconds (57 allocations: 44.039 KiB, 71.88% compilation time)
maximum(abs, S * y - transpose(S) * y) = 5.684341886080802e-14

All three operations are implemented in pure Julia and respect the block‑sparsity, so they are typically much faster than converting the matrix to a dense format first.

Converting to a Classical Sparse Matrix

Sometimes you need a SparseMatrixCSC. The conversion is straightforward. You can compare the memory footprints:

using SparseArrays
Ssparse = sparse(S)
display(Ssparse);

@show Base.format_bytes(Base.summarysize(S))
@show Base.format_bytes(Base.summarysize(Ssparse))

Base.format_bytes(Base.summarysize(S)) = "16.242 MiB"
Base.format_bytes(Base.summarysize(Ssparse)) = "46.581 MiB"

In specific examples SymmetricBlockMatrix consumes less memory.

Multi‑Threaded Construction (Optional)

Depending on the example, the block-coloring step can be a bottleneck, thus multithreading is switched-off by default. To enable multithreading you can determine a scheduler from OhMyThreads.

using OhMyThreads
S = SymmetricBlockMatrix(
    diagonals,
    selfvalues,
    offdiagonals,
    colvalues,
    rowvalues,
    sizeS;
    scheduler=DynamicScheduler(),
)