GitHub

BlockSparseMatrix

A BlockSparseMatrix is a matrices that is (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 block‑sparse matrix from a list 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 BlockSparseMatrix

using CompScienceMeshes, BEAST, H2Trees
using UnicodePlots
using BlockSparseMatrices

m = meshsphere(1.0, 0.1)
X = raviartthomas(m)
tree = TwoNTree(X, 0.2)
colvalues, rowvalues = H2Trees.nearinteractions(tree)

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

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

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

Matrix–Vector Products

y = randn(Float64, size(B, 1))
@time B * y
@time B' * y
@time transpose(B) * y
  0.335012 seconds (635.10 k allocations: 32.649 MiB, 3.41% gc time, 99.38% compilation time)
  0.067908 seconds (171.44 k allocations: 8.632 MiB, 96.82% compilation time)
  0.006797 seconds (30 allocations: 38.898 KiB, 69.53% compilation time)

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
Bsparse = sparse(B)

@show Base.format_bytes(Base.summarysize(B));
@show Base.format_bytes(Base.summarysize(Bsparse));
Base.format_bytes(Base.summarysize(B)) = "26.868 MiB"
Base.format_bytes(Base.summarysize(Bsparse)) = "46.581 MiB"

In specific examples BlockSparseMatrix 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
B = BlockSparseMatrix(
    blocks,
    colvalues,
    rowvalues,
    (numfunctions(X), numfunctions(X));
    scheduler=DynamicScheduler(),
)