[jorgensd]

Extend MatrixCSR in Python to support non-square block sizes matrixes and block sizes above 9.

The insertion kernels follow the pattern #4341 introduced for
spmv/spmvT in the same header: one loop body takes a BlockSizeArg
block size, either a std::integral_constant (compile-time, unrolled)
or an int (run-time). This replaces the template <int BS0, int BS1>
parameters on impl::insert_*. The public templates
(mat_add_values<BS0, BS1>() etc.) are unchanged and forward with
std::integral_constant.

Changes

• la::MatrixCSR: add/set/mat_add_values/mat_set_values
  handle every data/matrix block-size combination. Before, combinations
  such as data (1, 3) into a (2, 3) matrix passed the guard, then
  inserted wrong values in Release. The guard is removed.
• New overloads mat_add_values(bs0, bs1), mat_set_values(bs0, bs1),
  add(x, rows, cols, bs0, bs1) and set(...) take the block size as an
  int or a std::integral_constant. The template versions forward to
  them, so the compile-time path is unchanged.
• Wrappers: assemble_matrix and interpolation_matrix share a
  visit_block_size dispatch. It uses compile-time sizes (k, k) for
  k ∈ {1, 2, 3, 4, 6, 9}, plus (1, 2), (2, 1), (1, 3) and (3, 1), and
  falls back to int otherwise. assemble_matrix now dispatches on the
  DOF map bs(), which sets the data layout, rather than
  index_map_bs(). The two differ for a sub-space, which is a view into
  the parent space and assembles into the parent matrix.
• Python MatrixCSR.add/set accept any block size, not only 1–3.

Benchmarks

assemble_matrix into an existing MatrixCSR, P1 mass matrices and
couplings on 82,944 tetrahedra (Release, serial, aarch64):

Block size	compact: main → PR	expanded: main → PR
1, 2, 3, 4, 6, 9	within ±2%	within ±2%
5	0.0200 → 0.0204 s	0.0413 → 0.0449 s
7	0.0278 → 0.0307 s	0.0605 → 0.0670 s
8	0.0323 → 0.0376 s	0.0713 → 0.0809 s
(4, 4) tensor (16)	raises → 0.150 s	raises → 0.252 s
(1, 3) / (3, 1)	raises → 0.0060 / 0.0063 s	raises → 0.0076 / 0.0141 s

• Block sizes 5, 7 and 8 are 2–16% slower than on main because
  they are no longer compiled; restoring them removes the difference. In
  C++, int block sizes are up to 34% slower than compile-time ones
  across block sizes, hence the compiled set.
• interpolation_matrix (P1 → P2) matches main for value shapes
  (2,) and (3,) (±2%), and is about 5% slower for scalars
  (0.0339 → 0.0356 s). Block sizes 4 and 5 now work.

Tests

• Mixed block sizes: scalar–vector couplings with 2–4 components,
  with and without boundary conditions.
• Sub-spaces: W.sub(0) arguments, where bs() != index_map_bs(),
  assembled into the parent matrix.
• Large block size: a block-size-16 mass matrix.
• Interpolation: blocked interpolation_matrix with 2, 4 and 5
  components.
• MatrixCSR.add/set: block-size-4 data into matrices of block
  size 1, 2 and 4.

AI assistance: I used Claude Code (Claude Opus 5.5) to draft parts of this
PR. I reviewed, edited, tested, and take responsibility for the final
contribution.

🤖 Generated with Claude Code

[chrisrichardson]

Is the logic here correct? Why is it only negated for bs={1,1}?

[chrisrichardson]

OK, sorry the not is around everything, I see. Not very easy to read, however.

[chrisrichardson]

Maybe get rid of the ! and just return if true?