Julien Jerphanion
b47183a713
kernel/wasm/KERNEL is included after KERNEL.WASM128_GENERIC and
unconditionally pointed S/D GEMV at RISC-V scalar gemv_{n,t}.c, so a
target-file SIMD kernel would never run.
Add kernel/wasm/gemv_{n,t}.c for both precisions. GEMV_T keeps 8
(float) / 4 (double) independent column accumulators with IEEE mul+add
and delays the horizontal add until after the inner loop. GEMV_N
streams four columns into y with a four-lane v128 unroll. Non-unit
stride stays scalar (no WASM gather/scatter).
Guard SGEMV{N,T}KERNEL/DGEMV{N,T}KERNEL with ifndef in KERNEL so the
target file wins. CGEMV/ZGEMV stay RISC-V scalar.
Signed-off-by: Julien Jerphanion <git@jjerphan.xyz>
197 lines
6.9 KiB
C
197 lines
6.9 KiB
C
/***************************************************************************
|
|
Copyright (c) 2026, The OpenBLAS Project
|
|
All rights reserved.
|
|
Redistribution and use in source and binary forms, with or without
|
|
modification, are permitted provided that the following conditions are
|
|
met:
|
|
1. Redistributions of source code must retain the above copyright
|
|
notice, this list of conditions and the following disclaimer.
|
|
2. Redistributions in binary form must reproduce the above copyright
|
|
notice, this list of conditions and the following disclaimer in
|
|
the documentation and/or other materials provided with the
|
|
distribution.
|
|
3. Neither the name of the OpenBLAS project nor the names of
|
|
its contributors may be used to endorse or promote products
|
|
derived from this software without specific prior written permission.
|
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
|
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|
ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
|
|
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
|
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
|
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
|
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
|
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|
POSSIBILITY OF SUCH DAMAGE.
|
|
*****************************************************************************/
|
|
|
|
/*
|
|
* WASM SIMD128 GEMV_N: y += alpha * A * x (AXPY each column of A into y).
|
|
*
|
|
* Compiled twice: SGEMV_N (-UDOUBLE) and DGEMV_N (-DDOUBLE).
|
|
*
|
|
* Four columns share one streaming of y. Inner loop unrolls four v128 lanes
|
|
* (16 floats / 8 doubles). IEEE mul+add (not relaxed madd).
|
|
* Non-unit inc_y stays scalar (no scatter).
|
|
*/
|
|
|
|
#include "common.h"
|
|
|
|
#if defined(__wasm_simd128__)
|
|
#include <wasm_simd128.h>
|
|
|
|
#ifdef DOUBLE
|
|
#define GEMV_VLEN 2
|
|
#define GEMV_SPLAT wasm_f64x2_splat
|
|
#define GEMV_MUL wasm_f64x2_mul
|
|
#define GEMV_ADD wasm_f64x2_add
|
|
#else
|
|
#define GEMV_VLEN 4
|
|
#define GEMV_SPLAT wasm_f32x4_splat
|
|
#define GEMV_MUL wasm_f32x4_mul
|
|
#define GEMV_ADD wasm_f32x4_add
|
|
#endif
|
|
|
|
#define GEMV_UNROLL 4
|
|
#define GEMV_CHUNK (GEMV_VLEN * GEMV_UNROLL)
|
|
#define GEMV_LOAD(p) wasm_v128_load((const void *)(p))
|
|
#define GEMV_STORE(p, v) wasm_v128_store((void *)(p), (v))
|
|
#define GEMV_MADD(y, a, x) GEMV_ADD((y), GEMV_MUL((a), (x)))
|
|
|
|
static void gemv_n_axpy4(BLASLONG m, const FLOAT *a0, const FLOAT *a1,
|
|
const FLOAT *a2, const FLOAT *a3, FLOAT *y, FLOAT t0,
|
|
FLOAT t1, FLOAT t2, FLOAT t3) {
|
|
const v128_t v0 = GEMV_SPLAT(t0);
|
|
const v128_t v1 = GEMV_SPLAT(t1);
|
|
const v128_t v2 = GEMV_SPLAT(t2);
|
|
const v128_t v3 = GEMV_SPLAT(t3);
|
|
BLASLONG i = 0;
|
|
const BLASLONG n_main = m & ~(BLASLONG)(GEMV_CHUNK - 1);
|
|
|
|
for (; i < n_main; i += GEMV_CHUNK) {
|
|
v128_t acc;
|
|
acc = GEMV_LOAD(y + i + 0 * GEMV_VLEN);
|
|
acc = GEMV_MADD(acc, v0, GEMV_LOAD(a0 + i + 0 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v1, GEMV_LOAD(a1 + i + 0 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v2, GEMV_LOAD(a2 + i + 0 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v3, GEMV_LOAD(a3 + i + 0 * GEMV_VLEN));
|
|
GEMV_STORE(y + i + 0 * GEMV_VLEN, acc);
|
|
|
|
acc = GEMV_LOAD(y + i + 1 * GEMV_VLEN);
|
|
acc = GEMV_MADD(acc, v0, GEMV_LOAD(a0 + i + 1 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v1, GEMV_LOAD(a1 + i + 1 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v2, GEMV_LOAD(a2 + i + 1 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v3, GEMV_LOAD(a3 + i + 1 * GEMV_VLEN));
|
|
GEMV_STORE(y + i + 1 * GEMV_VLEN, acc);
|
|
|
|
acc = GEMV_LOAD(y + i + 2 * GEMV_VLEN);
|
|
acc = GEMV_MADD(acc, v0, GEMV_LOAD(a0 + i + 2 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v1, GEMV_LOAD(a1 + i + 2 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v2, GEMV_LOAD(a2 + i + 2 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v3, GEMV_LOAD(a3 + i + 2 * GEMV_VLEN));
|
|
GEMV_STORE(y + i + 2 * GEMV_VLEN, acc);
|
|
|
|
acc = GEMV_LOAD(y + i + 3 * GEMV_VLEN);
|
|
acc = GEMV_MADD(acc, v0, GEMV_LOAD(a0 + i + 3 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v1, GEMV_LOAD(a1 + i + 3 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v2, GEMV_LOAD(a2 + i + 3 * GEMV_VLEN));
|
|
acc = GEMV_MADD(acc, v3, GEMV_LOAD(a3 + i + 3 * GEMV_VLEN));
|
|
GEMV_STORE(y + i + 3 * GEMV_VLEN, acc);
|
|
}
|
|
for (; i + GEMV_VLEN <= m; i += GEMV_VLEN) {
|
|
v128_t acc = GEMV_LOAD(y + i);
|
|
acc = GEMV_MADD(acc, v0, GEMV_LOAD(a0 + i));
|
|
acc = GEMV_MADD(acc, v1, GEMV_LOAD(a1 + i));
|
|
acc = GEMV_MADD(acc, v2, GEMV_LOAD(a2 + i));
|
|
acc = GEMV_MADD(acc, v3, GEMV_LOAD(a3 + i));
|
|
GEMV_STORE(y + i, acc);
|
|
}
|
|
for (; i < m; i++)
|
|
y[i] += t0 * a0[i] + t1 * a1[i] + t2 * a2[i] + t3 * a3[i];
|
|
}
|
|
|
|
static void gemv_n_axpy1(BLASLONG m, const FLOAT *a, FLOAT *y, FLOAT t) {
|
|
const v128_t vt = GEMV_SPLAT(t);
|
|
BLASLONG i = 0;
|
|
const BLASLONG n_main = m & ~(BLASLONG)(GEMV_CHUNK - 1);
|
|
|
|
for (; i < n_main; i += GEMV_CHUNK) {
|
|
v128_t acc;
|
|
acc = GEMV_LOAD(y + i + 0 * GEMV_VLEN);
|
|
GEMV_STORE(y + i + 0 * GEMV_VLEN, GEMV_MADD(acc, vt, GEMV_LOAD(a + i + 0 * GEMV_VLEN)));
|
|
acc = GEMV_LOAD(y + i + 1 * GEMV_VLEN);
|
|
GEMV_STORE(y + i + 1 * GEMV_VLEN, GEMV_MADD(acc, vt, GEMV_LOAD(a + i + 1 * GEMV_VLEN)));
|
|
acc = GEMV_LOAD(y + i + 2 * GEMV_VLEN);
|
|
GEMV_STORE(y + i + 2 * GEMV_VLEN, GEMV_MADD(acc, vt, GEMV_LOAD(a + i + 2 * GEMV_VLEN)));
|
|
acc = GEMV_LOAD(y + i + 3 * GEMV_VLEN);
|
|
GEMV_STORE(y + i + 3 * GEMV_VLEN, GEMV_MADD(acc, vt, GEMV_LOAD(a + i + 3 * GEMV_VLEN)));
|
|
}
|
|
for (; i + GEMV_VLEN <= m; i += GEMV_VLEN) {
|
|
v128_t acc = GEMV_LOAD(y + i);
|
|
GEMV_STORE(y + i, GEMV_MADD(acc, vt, GEMV_LOAD(a + i)));
|
|
}
|
|
for (; i < m; i++)
|
|
y[i] += t * a[i];
|
|
}
|
|
|
|
#undef GEMV_MADD
|
|
#undef GEMV_STORE
|
|
#undef GEMV_LOAD
|
|
#undef GEMV_CHUNK
|
|
#undef GEMV_UNROLL
|
|
#undef GEMV_ADD
|
|
#undef GEMV_MUL
|
|
#undef GEMV_SPLAT
|
|
#undef GEMV_VLEN
|
|
#endif
|
|
|
|
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) {
|
|
BLASLONG i, j, ix, iy;
|
|
FLOAT *a_ptr;
|
|
FLOAT temp;
|
|
|
|
(void)dummy1;
|
|
(void)buffer;
|
|
|
|
if (m < 1 || n < 1)
|
|
return 0;
|
|
if (alpha == (FLOAT)0.0)
|
|
return 0;
|
|
|
|
#if defined(__wasm_simd128__)
|
|
if (inc_y == 1) {
|
|
ix = 0;
|
|
j = 0;
|
|
for (; j + 4 <= n; j += 4) {
|
|
gemv_n_axpy4(m, a + j * lda, a + (j + 1) * lda, a + (j + 2) * lda,
|
|
a + (j + 3) * lda, y, alpha * x[ix],
|
|
alpha * x[ix + inc_x], alpha * x[ix + 2 * inc_x],
|
|
alpha * x[ix + 3 * inc_x]);
|
|
ix += 4 * inc_x;
|
|
}
|
|
for (; j < n; j++) {
|
|
gemv_n_axpy1(m, a + j * lda, y, alpha * x[ix]);
|
|
ix += inc_x;
|
|
}
|
|
return 0;
|
|
}
|
|
#endif
|
|
|
|
ix = 0;
|
|
a_ptr = a;
|
|
for (j = 0; j < n; j++) {
|
|
temp = alpha * x[ix];
|
|
iy = 0;
|
|
for (i = 0; i < m; i++) {
|
|
y[iy] += temp * a_ptr[i];
|
|
iy += inc_y;
|
|
}
|
|
a_ptr += lda;
|
|
ix += inc_x;
|
|
}
|
|
return 0;
|
|
}
|