From 7b1ec93c473c2766e31fc99b46dbc19404f13e0f Mon Sep 17 00:00:00 2001 From: 6eanut Date: Fri, 18 Sep 2026 09:08:27 +0000 Subject: [PATCH] RISC-V: register-block the ZVL256B GEMV-N kernel The ZVL256B GEMV-N kernel for unit y stride walked one column of A at a time. For each column it loaded y, loaded the column, issued one vfmacc, stored y back, and then repeated the whole sequence for the next column. Every column therefore paid a y load and a y store round-trip, and because VL was recomputed from the remaining row count on each pass, the vsetvli and the pointer bumps also sat inside the hot loop. Process GEMV_N_BLOCK (4) consecutive columns per pass over y. A single y slice stays resident in vector registers while all four columns are accumulated into it, so the y round-trip is paid once per block instead of once per column, and the main loop's VL is pinned to the hardware maximum, which lets the compiler hoist the vsetvli and strength-reduce the pointer arithmetic out of the loop. The columns are still consumed in ascending j and the four contributions are folded into the accumulator in that same order, so each y[i] sees the same sequence of addends as the original column-at-a-time loop; the intermediate y vector merely stays in a register instead of being stored and reloaded. The row remainder of each block keeps a variable VL, and the columns left over after the last full block are finished by the original two-column and single-column loops, unchanged. The result is bit-identical. A differential run of the baseline and patched kernels against each other, and of both against a double-precision reference, over 21375 single-precision and 21375 double-precision cases covering m = 1..129, n = 1..33, several leading dimensions and x increments, and alpha values including 0, with guard pages around every buffer, shows no mismatch and no out-of-bounds access. The same sweep for a non-unit y stride, which the patch leaves untouched, is likewise bit-identical. The library built from the patched kernel passes 125/125 utests and 1473/1473 extension tests. Measured on one X100 core at 2.2 GHz, warm-up once and the median of three runs per size: sgemv 256 -1.15%, 512 +30.01%, 1024 +13.57%; dgemv 256 +31.04%, 512 +22.14%, 1024 +2.40%. Co-authored-by: Yuansheng Co-authored-by: Ning Tian Signed-off-by: jiakai xu --- kernel/riscv64/gemv_n_vector.c | 111 ++++++++++++++++++++++++++------- 1 file changed, 87 insertions(+), 24 deletions(-) diff --git a/kernel/riscv64/gemv_n_vector.c b/kernel/riscv64/gemv_n_vector.c index 845e4de296..678f6b3c93 100644 --- a/kernel/riscv64/gemv_n_vector.c +++ b/kernel/riscv64/gemv_n_vector.c @@ -29,6 +29,7 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #if !defined(DOUBLE) #define VSETVL(n) RISCV_RVV(vsetvl_e32m8)(n) +#define VSETVLMAX RISCV_RVV(vsetvlmax_e32m8)() #define FLOAT_V_T vfloat32m8_t #define VLEV_FLOAT RISCV_RVV(vle32_v_f32m8) #define VLSEV_FLOAT RISCV_RVV(vlse32_v_f32m8) @@ -37,6 +38,7 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #define VFMACCVF_FLOAT RISCV_RVV(vfmacc_vf_f32m8) #else #define VSETVL(n) RISCV_RVV(vsetvl_e64m8)(n) +#define VSETVLMAX RISCV_RVV(vsetvlmax_e64m8)() #define FLOAT_V_T vfloat64m8_t #define VLEV_FLOAT RISCV_RVV(vle64_v_f64m8) #define VLSEV_FLOAT RISCV_RVV(vlse64_v_f64m8) @@ -45,6 +47,13 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #define VFMACCVF_FLOAT RISCV_RVV(vfmacc_vf_f64m8) #endif +/* Number of consecutive columns whose contribution is accumulated while a + single y vector slice stays resident in vector registers. Keeping the y + slice in registers removes the per-column y load/store round-trip, while + the fixed VL used by the main loop lets the compiler hoist the vsetvli and + strength-reduce the pointer increments out of the hot loop. */ +#define GEMV_N_BLOCK 4 + int CNAME(BLASLONG m, BLASLONG n, BLASLONG dummy1, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y, FLOAT *buffer) { if (n < 0) return(0); @@ -54,39 +63,93 @@ int CNAME(BLASLONG m, BLASLONG n, BLASLONG dummy1, FLOAT alpha, FLOAT *a, BLASLO FLOAT_V_T va, vy, va2; if (inc_y == 1) { - for (j = 0; j < (n >> 1); j++) { - temp = alpha * x[0]; - temp2 = alpha * x[inc_x]; - y_ptr = y; - a_ptr = a; - a2_ptr = a + lda; - for (i = m; i > 0; i -= vl) { - vl = VSETVL(i); + const BLASLONG epr = VSETVLMAX; + + /* Register-blocked main loop. A single y vector slice is kept + resident in registers while GEMV_N_BLOCK consecutive columns are + accumulated into it, so the y load/store round-trip happens once + per block instead of once per column. VL is fixed to the hardware + maximum for the whole main loop, which lets the compiler hoist the + vsetvli and strength-reduce the pointer increments out of the hot + loop. Columns are still consumed in ascending order, therefore the + FP64 accumulation sequence of every y[i] is bit-identical to the + original column-at-a-time loop: the intermediate y vector is only + held in a register rather than stored and reloaded. */ + for (j = 0; j + GEMV_N_BLOCK <= n; j += GEMV_N_BLOCK) { + FLOAT *ap0 = a + (j + 0) * lda; + FLOAT *ap1 = a + (j + 1) * lda; + FLOAT *ap2 = a + (j + 2) * lda; + FLOAT *ap3 = a + (j + 3) * lda; + FLOAT t0 = alpha * x[(j + 0) * inc_x]; + FLOAT t1 = alpha * x[(j + 1) * inc_x]; + FLOAT t2 = alpha * x[(j + 2) * inc_x]; + FLOAT t3 = alpha * x[(j + 3) * inc_x]; + + for (i = 0; i + epr <= m; i += epr) { + y_ptr = y + i; + vy = VLEV_FLOAT(y_ptr, epr); + vy = VFMACCVF_FLOAT(vy, t0, VLEV_FLOAT(ap0 + i, epr), epr); + vy = VFMACCVF_FLOAT(vy, t1, VLEV_FLOAT(ap1 + i, epr), epr); + vy = VFMACCVF_FLOAT(vy, t2, VLEV_FLOAT(ap2 + i, epr), epr); + vy = VFMACCVF_FLOAT(vy, t3, VLEV_FLOAT(ap3 + i, epr), epr); + VSEV_FLOAT(y_ptr, vy, epr); + } + if (i < m) { + vl = VSETVL(m - i); + y_ptr = y + i; vy = VLEV_FLOAT(y_ptr, vl); - va = VLEV_FLOAT(a_ptr, vl); - va2 = VLEV_FLOAT(a2_ptr, vl); + vy = VFMACCVF_FLOAT(vy, t0, VLEV_FLOAT(ap0 + i, vl), vl); + vy = VFMACCVF_FLOAT(vy, t1, VLEV_FLOAT(ap1 + i, vl), vl); + vy = VFMACCVF_FLOAT(vy, t2, VLEV_FLOAT(ap2 + i, vl), vl); + vy = VFMACCVF_FLOAT(vy, t3, VLEV_FLOAT(ap3 + i, vl), vl); + VSEV_FLOAT(y_ptr, vy, vl); + } + } + + /* Remaining columns keep the original two-column (n >> 1) structure. */ + for (; j + 2 <= n; j += 2) { + temp = alpha * x[(j + 0) * inc_x]; + temp2 = alpha * x[(j + 1) * inc_x]; + a_ptr = a + (j + 0) * lda; + a2_ptr = a + (j + 1) * lda; + for (i = 0; i + epr <= m; i += epr) { + y_ptr = y + i; + vy = VLEV_FLOAT(y_ptr, epr); + va = VLEV_FLOAT(a_ptr + i, epr); + va2 = VLEV_FLOAT(a2_ptr + i, epr); + vy = VFMACCVF_FLOAT(vy, temp, va, epr); + vy = VFMACCVF_FLOAT(vy, temp2, va2, epr); + VSEV_FLOAT(y_ptr, vy, epr); + } + if (i < m) { + vl = VSETVL(m - i); + y_ptr = y + i; + vy = VLEV_FLOAT(y_ptr, vl); + va = VLEV_FLOAT(a_ptr + i, vl); + va2 = VLEV_FLOAT(a2_ptr + i, vl); vy = VFMACCVF_FLOAT(vy, temp, va, vl); vy = VFMACCVF_FLOAT(vy, temp2, va2, vl); VSEV_FLOAT(y_ptr, vy, vl); - y_ptr += vl; - a_ptr += vl; - a2_ptr += vl; } - x += inc_x * 2; - a += lda * 2; } - if (n & 1) { - temp = alpha * x[0]; - y_ptr = y; - a_ptr = a; - for (i = m; i > 0; i -= vl) { - vl = VSETVL(i); + /* Remaining odd column keeps the original single-column structure. */ + if (j < n) { + temp = alpha * x[j * inc_x]; + a_ptr = a + j * lda; + for (i = 0; i + epr <= m; i += epr) { + y_ptr = y + i; + vy = VLEV_FLOAT(y_ptr, epr); + va = VLEV_FLOAT(a_ptr + i, epr); + vy = VFMACCVF_FLOAT(vy, temp, va, epr); + VSEV_FLOAT(y_ptr, vy, epr); + } + if (i < m) { + vl = VSETVL(m - i); + y_ptr = y + i; vy = VLEV_FLOAT(y_ptr, vl); - va = VLEV_FLOAT(a_ptr, vl); + va = VLEV_FLOAT(a_ptr + i, vl); vy = VFMACCVF_FLOAT(vy, temp, va, vl); VSEV_FLOAT(y_ptr, vy, vl); - y_ptr += vl; - a_ptr += vl; } } } else {