From 533165cc63e4cb5a968e59107b446cdd974602cf Mon Sep 17 00:00:00 2001 From: 6eanut Date: Thu, 17 Sep 2026 06:53:35 +0000 Subject: [PATCH] RISC-V: unroll the ZVL256B SGEMM kernel K-loop by two The ZVL256B SGEMM main pass advanced one k-step per iteration: each iteration loaded eight B scalars and one 16-lane A vector, issued eight accumulating FMAs, then bumped the B and A pointers and took the loop backedge. On the tall-and-narrow shapes this kernel is used on, that per-k loop control is a measurable fraction of the iteration. Fold two consecutive k-steps into one iteration. Both A vectors are loaded before the FMA stream (A double-buffering), so the A-load latency is exposed earlier, and the k counter, one of the two B/A pointer bumps and the backedge are halved. Each accumulator still consumes the two k-steps of a pair in ascending order, so the per-column accumulation order is unchanged. The peeled first k-step (the initial vfmul, ahead of the loop) is left alone, so the restructured main pass reads for (k = K; k > 2; k -= 2) { ... two k-steps ... } for (; --k; ) { ... one k-step ... } which consumes exactly K steps for every K: the unrolled loop runs floor((K-1)/2) times and the trailing loop takes the remainder. The unrolled loop is skipped entirely when K <= 2, in which case the trailing loop alone (or, for K = 1, the peeled step alone) finishes the pass. Because the accumulation order itself is preserved rather than merely re-summed, the emitted results are bit-identical to the previous kernel. Measured on one X100 core at 2.2 GHz, cblas_sgemm, warm-up once and the median of three runs per size: 256x256 +5.13%, 512x512 +5.56%, 1024x1024 +3.43%. In a default RISCV64_ZVL256B build (USE_TRMM=1) this source is compiled once, for SGEMM; STRMM is built from its own strmm_kernel_16x8_zvl256b.c. The Makefile's USE_TRMM=0 branch instead reuses this source for the four STRMM variants, so the patched file was also compiled in those four -DTRMMKERNEL configurations: all build with a warning set unchanged from the unpatched tree. Verified further by a differential run against a double-precision reference: 360 packed-kernel invocations and 20400 cblas_sgemm cases, covering the M/N/K tail and edge paths, all four transpose combinations and several alpha/beta values, produce output bit-identical to the unpatched kernel. Co-authored-by: Yuansheng Co-authored-by: Ning Tian Signed-off-by: jiakai xu --- kernel/riscv64/sgemm_kernel_16x8_zvl256b.c | 52 +++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/kernel/riscv64/sgemm_kernel_16x8_zvl256b.c b/kernel/riscv64/sgemm_kernel_16x8_zvl256b.c index 104a1a0ac2..4782270a1d 100644 --- a/kernel/riscv64/sgemm_kernel_16x8_zvl256b.c +++ b/kernel/riscv64/sgemm_kernel_16x8_zvl256b.c @@ -2170,7 +2170,57 @@ int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, FLOAT alpha, FLOAT* A, FLOAT* B, F vfloat32m2_t resultCD = __riscv_vfmul_vf_f32m2( A00, B6, 16 ); vfloat32m2_t resultEF = __riscv_vfmul_vf_f32m2( A00, B7, 16 ); - for (BLASLONG k = K; --k; ) { + /* 2x K-unroll with A double-buffering: two consecutive k-steps are + * folded into one iteration so the per-k loop control (counter, + * B/A pointer bumps, backedge) is halved, and the two A vectors + * are loaded before the FMA stream to expose A-load latency. + * Each accumulator still consumes the k-steps in the original + * ascending order (step p then step p+1), so the per-column FP + * accumulation order is unchanged. */ + BLASLONG k; + for (k = K; k > 2; k -= 2) { + FLOAT B8 = B[8]; + FLOAT B9 = B[9]; + FLOAT BA = B[10]; + FLOAT BB = B[11]; + FLOAT BC = B[12]; + FLOAT BD = B[13]; + FLOAT BE = B[14]; + FLOAT BF = B[15]; + + B0 = B[0]; + B1 = B[1]; + B2 = B[2]; + B3 = B[3]; + B4 = B[4]; + B5 = B[5]; + B6 = B[6]; + B7 = B[7]; + B += 16; + + vfloat32m2_t A0 = __riscv_vle32_v_f32m2( A, 16 ); + vfloat32m2_t A1 = __riscv_vle32_v_f32m2( A + 16, 16 ); + A += 32; + + result01 = __riscv_vfmacc_vf_f32m2( result01, B0, A0, 16 ); + result23 = __riscv_vfmacc_vf_f32m2( result23, B1, A0, 16 ); + result45 = __riscv_vfmacc_vf_f32m2( result45, B2, A0, 16 ); + result67 = __riscv_vfmacc_vf_f32m2( result67, B3, A0, 16 ); + result89 = __riscv_vfmacc_vf_f32m2( result89, B4, A0, 16 ); + resultAB = __riscv_vfmacc_vf_f32m2( resultAB, B5, A0, 16 ); + resultCD = __riscv_vfmacc_vf_f32m2( resultCD, B6, A0, 16 ); + resultEF = __riscv_vfmacc_vf_f32m2( resultEF, B7, A0, 16 ); + + result01 = __riscv_vfmacc_vf_f32m2( result01, B8, A1, 16 ); + result23 = __riscv_vfmacc_vf_f32m2( result23, B9, A1, 16 ); + result45 = __riscv_vfmacc_vf_f32m2( result45, BA, A1, 16 ); + result67 = __riscv_vfmacc_vf_f32m2( result67, BB, A1, 16 ); + result89 = __riscv_vfmacc_vf_f32m2( result89, BC, A1, 16 ); + resultAB = __riscv_vfmacc_vf_f32m2( resultAB, BD, A1, 16 ); + resultCD = __riscv_vfmacc_vf_f32m2( resultCD, BE, A1, 16 ); + resultEF = __riscv_vfmacc_vf_f32m2( resultEF, BF, A1, 16 ); + } + for (; --k; ) { B0 = B[0]; B1 = B[1]; B2 = B[2];