diff --git a/kernel/wasm/KERNEL b/kernel/wasm/KERNEL index 0c8dcfcc29..332db5322a 100644 --- a/kernel/wasm/KERNEL +++ b/kernel/wasm/KERNEL @@ -40,8 +40,12 @@ DSUMKERNEL = ../arm/sum.c CSUMKERNEL = ../arm/zsum.c ZSUMKERNEL = ../arm/zsum.c +ifndef SAXPYKERNEL SAXPYKERNEL = ../riscv64/axpy.c +endif +ifndef DAXPYKERNEL DAXPYKERNEL = ../riscv64/axpy.c +endif CAXPYKERNEL = ../riscv64/zaxpy.c ZAXPYKERNEL = ../riscv64/zaxpy.c diff --git a/kernel/wasm/KERNEL.WASM128_GENERIC b/kernel/wasm/KERNEL.WASM128_GENERIC index 0265085b1c..0b873ef3f7 100644 --- a/kernel/wasm/KERNEL.WASM128_GENERIC +++ b/kernel/wasm/KERNEL.WASM128_GENERIC @@ -40,8 +40,8 @@ DSUMKERNEL = ../arm/sum.c CSUMKERNEL = ../arm/zsum.c ZSUMKERNEL = ../arm/zsum.c -SAXPYKERNEL = ../x86_64/saxpy.c -DAXPYKERNEL = ../x86_64/daxpy.c +SAXPYKERNEL = axpy.c +DAXPYKERNEL = axpy.c CAXPYKERNEL = ../riscv64/zaxpy.c ZAXPYKERNEL = ../riscv64/zaxpy.c diff --git a/kernel/wasm/axpy.c b/kernel/wasm/axpy.c new file mode 100644 index 0000000000..dc2726d134 --- /dev/null +++ b/kernel/wasm/axpy.c @@ -0,0 +1,175 @@ +/*************************************************************************** +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 AXPY: y += da * x + * + * Compiled twice: SAXPY (-UDOUBLE) and DAXPY (-DDOUBLE). Unit-stride uses + * eight independent v128 lanes (32 floats / 16 doubles) so load/mul/add/store + * can overlap; remainder is one vector then scalar. IEEE mul+add (not + * relaxed madd): AXPY is checked to machine epsilon, and relaxed madd in + * the generic V_SIMD path previously slowed L1. + * + * Non-unit stride stays scalar (WASM SIMD128 has no gather). inc==0 is the + * scalar path so y[0] += n * da * x[0] still holds. + */ + +#include "common.h" + +#if defined(__wasm_simd128__) +#include + +#ifdef DOUBLE +#define AXPY_VLEN 2 +#define AXPY_SPLAT wasm_f64x2_splat +#define AXPY_MUL wasm_f64x2_mul +#define AXPY_ADD wasm_f64x2_add +#else +#define AXPY_VLEN 4 +#define AXPY_SPLAT wasm_f32x4_splat +#define AXPY_MUL wasm_f32x4_mul +#define AXPY_ADD wasm_f32x4_add +#endif + +#define AXPY_UNROLL 8 +#define AXPY_CHUNK (AXPY_VLEN * AXPY_UNROLL) + +#define AXPY_LOAD(p) wasm_v128_load((const void *)(p)) +#define AXPY_STORE(p, v) wasm_v128_store((void *)(p), (v)) +#define AXPY_MADD(y, a, x) AXPY_ADD((y), AXPY_MUL((a), (x))) + +static void axpy_kernel_unit(BLASLONG n, const FLOAT *x, FLOAT *y, FLOAT da) { + const v128_t va = AXPY_SPLAT(da); + BLASLONG i = 0; + const BLASLONG n_main = n & ~(BLASLONG)(AXPY_CHUNK - 1); + + for (; i < n_main; i += AXPY_CHUNK) { + v128_t x0 = AXPY_LOAD(x + i + 0 * AXPY_VLEN); + v128_t x1 = AXPY_LOAD(x + i + 1 * AXPY_VLEN); + v128_t x2 = AXPY_LOAD(x + i + 2 * AXPY_VLEN); + v128_t x3 = AXPY_LOAD(x + i + 3 * AXPY_VLEN); + v128_t x4 = AXPY_LOAD(x + i + 4 * AXPY_VLEN); + v128_t x5 = AXPY_LOAD(x + i + 5 * AXPY_VLEN); + v128_t x6 = AXPY_LOAD(x + i + 6 * AXPY_VLEN); + v128_t x7 = AXPY_LOAD(x + i + 7 * AXPY_VLEN); + + v128_t y0 = AXPY_LOAD(y + i + 0 * AXPY_VLEN); + v128_t y1 = AXPY_LOAD(y + i + 1 * AXPY_VLEN); + v128_t y2 = AXPY_LOAD(y + i + 2 * AXPY_VLEN); + v128_t y3 = AXPY_LOAD(y + i + 3 * AXPY_VLEN); + v128_t y4 = AXPY_LOAD(y + i + 4 * AXPY_VLEN); + v128_t y5 = AXPY_LOAD(y + i + 5 * AXPY_VLEN); + v128_t y6 = AXPY_LOAD(y + i + 6 * AXPY_VLEN); + v128_t y7 = AXPY_LOAD(y + i + 7 * AXPY_VLEN); + + AXPY_STORE(y + i + 0 * AXPY_VLEN, AXPY_MADD(y0, va, x0)); + AXPY_STORE(y + i + 1 * AXPY_VLEN, AXPY_MADD(y1, va, x1)); + AXPY_STORE(y + i + 2 * AXPY_VLEN, AXPY_MADD(y2, va, x2)); + AXPY_STORE(y + i + 3 * AXPY_VLEN, AXPY_MADD(y3, va, x3)); + AXPY_STORE(y + i + 4 * AXPY_VLEN, AXPY_MADD(y4, va, x4)); + AXPY_STORE(y + i + 5 * AXPY_VLEN, AXPY_MADD(y5, va, x5)); + AXPY_STORE(y + i + 6 * AXPY_VLEN, AXPY_MADD(y6, va, x6)); + AXPY_STORE(y + i + 7 * AXPY_VLEN, AXPY_MADD(y7, va, x7)); + } + + for (; i + AXPY_VLEN <= n; i += AXPY_VLEN) { + v128_t yi = AXPY_LOAD(y + i); + v128_t xi = AXPY_LOAD(x + i); + AXPY_STORE(y + i, AXPY_MADD(yi, va, xi)); + } + + for (; i < n; i++) + y[i] += da * x[i]; +} + +#undef AXPY_LOAD +#undef AXPY_STORE +#undef AXPY_MADD +#undef AXPY_CHUNK +#undef AXPY_UNROLL +#undef AXPY_VLEN +#undef AXPY_SPLAT +#undef AXPY_MUL +#undef AXPY_ADD +#endif + +int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, + BLASLONG inc_x, FLOAT *y, BLASLONG inc_y, FLOAT *dummy, + BLASLONG dummy2) { + BLASLONG i = 0; + BLASLONG ix = 0, iy = 0; + + (void)dummy0; + (void)dummy1; + (void)dummy; + (void)dummy2; + + if (n <= 0) + return 0; + if (da == 0.0) + return 0; + + if ((inc_x == 1) && (inc_y == 1)) { +#if defined(__wasm_simd128__) + axpy_kernel_unit(n, x, y, da); +#else + while (i < n) { + y[i] += da * x[i]; + i++; + } +#endif + return 0; + } + + { + BLASLONG n1 = n & ~(BLASLONG)3; + while (i < n1) { + FLOAT m1 = da * x[ix]; + FLOAT m2 = da * x[ix + inc_x]; + FLOAT m3 = da * x[ix + 2 * inc_x]; + FLOAT m4 = da * x[ix + 3 * inc_x]; + + y[iy] += m1; + y[iy + inc_y] += m2; + y[iy + 2 * inc_y] += m3; + y[iy + 3 * inc_y] += m4; + + ix += inc_x * 4; + iy += inc_y * 4; + i += 4; + } + } + + while (i < n) { + y[iy] += da * x[ix]; + ix += inc_x; + iy += inc_y; + i++; + } + return 0; +}