Replace the GEMV-loop ?GEMMT with a blocked level-3 driver

[jgillis]

(written by Claude Opus 5)

interface/gemmt.c has implemented ?GEMMT / ?GEMMTR as a loop over the columns of
the triangle with one ?GEMV per column ever since the routine was added in 0.3.22.
That is correct but, as #4921 notes, far off the performance of GEMM, and above the
GEMV multithreading threshold every single column becomes its own thread dispatch, so
on a many-core machine it gains almost nothing from extra cores and ends up far slower
than computing the whole product with GEMM and throwing half of it away.

This replaces it with a real level-3 driver. No new microkernels are needed: the blocks
that cross the diagonal are handled by the existing SYRK kernel, which already clips a
block against the diagonal at register-block granularity.

What is added

• driver/level3/level3_gemmt.c — the driver, derived from level3_syrk.c. The
  blocking, the diagonal handling and the loop structure are the same; the difference is
  that the packed B buffer is filled from a second matrix, so the buffer sharing that
  SYRK can do (shared) is not available, and op(A) / op(B) are selected by the same
  NN/NT/…/CC defines that level3.c uses.
• driver/level3/gemmt_k.c — the per-variant wrapper: the triangular beta scaling and
  the choice of triangular kernel. Built once per uplo × op(A) × op(B), exactly as
  syrk_k.c is built per uplo × trans.
• For complex, syrk_kernel.c is additionally built with CONJA / CONJB
  (?gemmt_kernel_UCN, _UNC, _UCC and the L counterparts); the file already
  supported those defines but nothing instantiated them.

interface/gemmt.c now sets up blas_arg_t and dispatches through a
[uplo][op(B)][op(A)] table like interface/gemm.c. Above the existing SMP threshold it
splits the triangle into column ranges of equal area with syrk_thread() — the
partitioner SYRK and SYR2K already use — so every thread runs the driver on a disjoint
set of columns of C.

Only the blocks that intersect the requested triangle are packed and computed, and only
the triangle is scaled by beta, so no work is wasted on the half of C that is not
referenced.

Side effect: complex GEMMT no longer writes to B

For transb = 'C' / 'R' the old code conjugated the caller's B in place
(IMATCOPY_K_CNC) and conjugated it back at the end. That made a complex GEMMT unsafe
against a B that is shared between threads or that lives in read-only memory. Conjugation
is now part of the kernel selection, so all input matrices are left untouched — the tests
below check that explicitly.

Numbers

dgemmt, 8-core/16-thread i7-11800H, GCC 13, TARGET=SKYLAKEX, time for one call
(best of ≥5), and dsyrk on the same data for reference:

n	k	threads	before	after	dsyrk
400	128	1	1.05 ms	0.36 ms	0.36 ms
800	256	1	9.30 ms	2.60 ms	2.66 ms
1600	128	1	18.7 ms	5.14 ms	5.12 ms
1600	256	1	47.0 ms	10.7 ms	10.7 ms
400	128	4	0.82 ms	0.14 ms	0.17 ms
1600	128	4	6.92 ms	1.54 ms	1.53 ms
400	128	16	1.31 ms	0.10 ms	0.11 ms
800	256	16	6.92 ms	0.66 ms	0.64 ms
1600	128	16	11.7 ms	1.01 ms	1.01 ms
1600	256	16	21.2 ms	2.06 ms	2.05 ms

Measurements are interleaved old/new/syrk, each the best of three rounds of a separate
process (each round the best of >= 5 calls), on an otherwise idle machine.

Single-threaded the new driver is 3–4× faster than the loop and matches SYRK. The old
implementation barely benefits from more threads (18.7 → 11.7 ms from 1 to 16 threads at
n=1600), while the new one scales like SYRK (18.7 → 1.01 ms), so the gap grows to ~10–13×
on 16 threads. Under load, or with several GEMMT calls in flight, the old per-column
dispatch degrades much further — that is how it was noticed downstream, where MUMPS calls
GEMMT for every symmetric front update when the BLAS provides it
(casadi/casadi#4417).

Testing

• utest, ctest and test pass (make -C utest: 1539/1539, including four new
  cases, see below). Checked with 1, 8 and 16 threads.
• A randomized differential harness against a NumPy reference: every precision × uplo ×
  op(A) × op(B) (including the CBLAS-only ConjNoTrans) × Fortran/CBLAS ×
  column/row-major, with padded leading dimensions and alpha/beta ∈ {0, 1, −1,
  random, pure imaginary}, over n ∈ {0…200} and k ∈ {0…129} — 39 644 cases, plus 6 360
  large cases (n up to 700) at 2, 3, 5, 7 and 16 threads. It also checks that the
  opposite triangle, the ldc padding, the memory around every buffer and the input
  matrices A and B are left untouched.
• Verified across kernel families with a DYNAMIC_ARCH build and OPENBLAS_CORETYPE =
  PRESCOTT, NEHALEM, SANDYBRIDGE, HASWELL, SKYLAKEX, ZEN (their register blocking, and
  hence the diagonal clipping, differ).
• valgrind: clean (0 errors, no leaks) over all variants at 1 and 2 threads.
• The drivers called directly with explicit row and column ranges (which the interface
  never passes), for both triangles and beta != 1, and the ?gemmtr entry points
  through the same harness.
• Builds checked: Makefile default, DYNAMIC_ARCH=1, USE_THREAD=0, and CMake.

The four new utest cases (?gemmt.upper_M_257_K_64_a_notrans_b_notrans) use a size
that is not a multiple of the register blocking and is large enough to be split over
threads, which is the case the existing tests (capped at 100) did not reach; they compare
against GEMM and assert that the lower triangle is bit-identical to its input.

Note on level3_syrk.c

For a row band that lies entirely below the column block, level3_syrk.c bounds its column
loop with for (jjs = js; jjs < min_j; …), mixing an absolute index with a length: when
js > 0 it packs nothing and the kernel then reads an unfilled buffer. The GEMMT driver
uses jjs < js + min_j. As far as I can see nothing in the tree calls the SYRK driver with a
row range, so the branch is unreachable today; I left level3_syrk.c alone to keep this PR
to GEMMT, but it is a one-line fix if you would like it here.

🤖 Generated with Claude Code