[velonica0]

Description

In kernel/riscv64/zgemm_kernel_8x4_zvl256b.c, the main 8x4 tile updates C (C += alpha * ACC) with stride-16-byte accesses: per tile, 16 vlse64 (one for the real parts and one for the imaginary parts of each of the 8 C segments) and 16 vsse64. Each cache line of C is therefore walked twice on load and twice on store.

C is stored interleaved (re, im). A 2-field segment access (vlseg2e64 / vsseg2e64) moves gvl contiguous complex elements in one instruction and de-interleaves them into the same real and imaginary vectors. #6051 already does this for the A panel in this file. This PR does the same for C.

Why it matters: on a single-thread 512³ ZGEMM on the X100, the C update block takes 24.2% of the samples in zgemm_kernel_n. The K loop has no spills and is not touched. With ZGEMM_DEFAULT_Q=64, the C update runs once every 64 K iterations, so at K=512 the whole of C is streamed 8 times.

Change

For each of the 8 C segments of the main 8x4 tile:

vfloat64m1x2_t Cseg0 = __riscv_vlseg2e64_v_f64m1x2( &C[ci*2], gvl );
vfloat64m1_t C0r = __riscv_vget_v_f64m1x2_f64m1( Cseg0, 0 );
vfloat64m1_t C0i = __riscv_vget_v_f64m1x2_f64m1( Cseg0, 1 );
C0r = __riscv_vfmacc ( C0r, alphar, ACC0r, gvl );
C0i = __riscv_vfmacc ( C0i, alphar, ACC0i, gvl );
C0r = __riscv_vfnmsac( C0r, alphai, ACC0i, gvl );
C0i = __riscv_vfmacc ( C0i, alphai, ACC0r, gvl );
__riscv_vsseg2e64_v_f64m1x2( &C[ci*2], __riscv_vcreate_v_f64m1x2( C0r, C0i ), gvl );

Each segment is loaded, updated and stored before the next one is loaded. My first version loaded all 8 C tuples up front, as the old code did. With the 16 accumulators also live, that ran out of vector registers: GCC emitted 5 spill/reload pairs (vs1r/vl1re64), and zgemm_kernel_n grew from 1470 to 1933 instructions. Interleaving the updates keeps one C tuple live at a time.

Only the main 8x4 path changes. The M/N edge tiles keep their strided accesses and are left for a follow-up.

Why the results are bit-identical

• The 8 C segments are disjoint: two per column, at ci and ci+gvl with gvl = 4, and the columns are ldc >= M apart. Moving each store before the next segment's load therefore cannot change any value that is read.
• vlseg2e64 field 0 / field 1 give exactly the elements that the old vlse64 at &C[2ci] / &C[2ci+1] with stride 16 gave, and vsseg2e64 writes them back to the same addresses.
• Each element goes through the same 4 FMAs (vfmacc by alphar, then vfnmsac/vfmacc by alphai), on the same operands and in the same order. No reassociation is involved.

Validation

Differential test. I wrote a standalone program, linked once against the baseline libopenblas.a and once against the patched one, and compared their raw output bytes with cmp. It covers 33040 cases:

• cblas_zgemm: both orders, all 16 (trans, conj) combinations for A and B, M, N ∈ {1,3,4,7,8,9,12,15,16,17,33,65}, K ∈ {1,5,8,63,64,65,129}, padded lda/ldb/ldc, three alpha and three beta values including 0 and pure imaginary;
• cblas_ztrsm: all side/uplo/trans;
• cblas_zherk: both uplo, both trans.

Every A/B/C buffer is mmaped so that it ends exactly at a PROT_NONE guard page, so any read or write past the end faults. The outputs are bit-identical on both X100 and A100 (178 MB, sha256 251ed73a50fcd0a8…), and nothing faulted.

Tests (CMake build, ctest, on each core type):

test	X100	A100
blas-smoke (openblas_utest, openblas_utest_ext)	pass	same 4 failures as the unpatched baseline (see below)
xzcblat3 (default threads)	pass	pass
xzcblat3 (1 thread)	pass	pass
utest zgemm	94/94	94/94

On the A100 the unpatched baseline already fails matrix_packing:ssymm, matrix_packing:dsymm, kernel_regress:skx_avx and dgemmt:upper_M_257_K_64_a_notrans_b_notrans. The patched build fails the same four, with identical messages. All four are real precision and none goes through this kernel. LAPACK tests were not run (the build has no Fortran).

Disassembly of zgemm_kernel_n (GCC 15.2.0, -march=rv64imafdcv_zvfh_zfh_zvl256b):

	insns	vlse64	vsse64	vlseg2e64	vsseg2e64	spills	vmv1r
baseline	1470	100	56	6	0	0	21
this PR	1464	84	40	14	8	0	44

The extra vmv1r come from GCC copying the tuple into and out of a second register pair around the FMAs: 4 per segment, once per tile. GCC 15.2 emits the same copies whether the tuple is rebuilt with vcreate, updated with vset, or updated with nested vset(vget(...)), so I left the source in its most readable form. Their cost is included in the numbers below.

Performance

SpacemiT K3: X100 core (VLEN 256, 2.2 GHz) and A100 core (VLEN 1024, 1.8 GHz). Single thread, pinned, TARGET=RISCV64_ZVL256B.

Size sweep, benchmark/gemm.c (zgemm NN, M=N=K), median of 5 paired rounds in alternating order, MFlops:

core	n	baseline	this PR	gain
X100	256	11890	12100	+1.8%
X100	512	10184	11974	+17.6%
X100	1024	10593	12105	+14.3%
A100	256	1168	1221	+4.5%
A100	512	1208	1237	+2.5%
A100	1024	1215	1247	+2.6%

The X100 gain is small at n=256 and 14–18% from n=512 up. I have not profiled why it is size-dependent. The A100 gains are small but consistent. Every cell's run-to-run spread is within 1%, except the A100 n=256 baseline (2.7%).

Other kernels that use the same tile path (median of paired rounds, bootstrap 95% CI):

case	X100	A100
zgemm 512³ (kernel_n)	+14.9% [14.8, 15.0]	+2.4% [2.4, 2.4]
zgemm conj-A variants (kernel_l)	+14.7%	+2.7%
zgemm conj-B variants (kernel_r)	+15.1%	+2.5%
zgemm conj-both variants (kernel_b)	+14.9%	+2.8%
zherk (via syrk → GEMM_KERNEL)	+13.1%	+2.5%
ztrsm (via trsm_kernel → GEMM_KERNEL)	+9.1%	+2.8%
dgemm, dtrsm (controls, unchanged)	±0.03%	±0.05%

The 8-thread zgemm results match the single-thread ones (X100 +14.9%, A100 +2.4%).