[SmaCow666]

Summary

Fix incorrect zero-output-stride accumulation in the RISC-V RVV AXPY kernel.

For incy == 0, the current RVV implementation broadcasts the initial output value to all active vector lanes, accumulates the corresponding input values, and then reduces all lanes. This causes the initial output value to be included once per active lane rather than once for the whole operation.

This change handles the zero-output-stride path with an ordered scalar fused accumulation, while leaving the regular RVV AXPY paths unchanged.

Minimal reproducer

For example:

n = 2;
alpha = 1;
x = {1, 1};
incx = 1;
y[0] = 2;
incy = 0;

The sequential update is:

y = 2 + 1 = 3
y = 3 + 1 = 4

The generic RISC-V implementation produces 4, while the affected RVV path produces 6.

For one RVV chunk with active length L, the affected implementation effectively computes

L * initial_y + alpha * sum(x)

instead of

initial_y + alpha * sum(x)

so the extra contribution is

(L - 1) * initial_y

The n = 2 reproducer uses exactly representable values and does not depend on floating-point rounding or multiple vector-length iterations.

Root cause

The incy == 0 RVV path initializes every active vector lane with the same initial y[0] value before reducing the lanes.

The initial output value is not a neutral reduction element and must only contribute once. Broadcasting it to all active lanes therefore duplicates the reduction seed.

Changing LMUL changes the active vector length for some inputs but does not address the underlying issue.

Fix

For incy == 0, perform the accumulation in program order using scalar fused multiply-add operations.

Other stride cases continue to use the existing implementation.

The change is limited to the RISC-V ZVL256B AXPY implementation and does not change the regular unit-stride RVV path.

Tests

Added two regression tests in utest/test_axpy.c, covering 240 cases per precision with varying lengths, input strides, initial output values, and alpha values.

The regression tests fail with the previous RVV implementation and pass with this fix.

Tested with:

• single- and double-precision AXPY;
• zero and nonzero strides across the existing extended test set;
• VLEN=256 and VLEN=512;
• OpenBLAS unit tests and BLAS Level-1 tests;
• guard-page and input-preservation checks.

The current upstream m4 kernel file was rebuilt with the fixed support objects and validated in QEMU at VLEN=256 and VLEN=512.

Native validation on Milk-V Jupiter covered the OpenBLAS 0.3.34-based version of this correction.

The generic implementation and the fixed RVV implementation both return 4 for the minimal reproducer. Broader zero-output-stride tests were checked against an ordered scalar fused-accumulation reference.

Scope

This PR addresses the behavior of the existing OpenBLAS zero-output-stride path. It does not change the normal incy != 0 AXPY path or introduce a performance optimization.

I have intentionally kept this correctness fix separate from unrelated short-vector AXPY performance work.

incy == 0 is outside the usual positive-stride BLAS usage, so I would also appreciate confirmation from the maintainers on the intended OpenBLAS compatibility contract for zero output stride.

[martin-frbg]

Thank you - zero output stride looks unusual but is allowed by the standard, so OpenBLAS should behave. (It's unfortunate that the tests taken from Reference-BLAS do not cover this case, will be interesting :( to see if other kernels are similarly flawed)

[martin-frbg]

Looks like it's only MIPS64/LOONGARCH64 axpy kernels that have this flaw, and the non-ZVL256B branch of the file your PR patches. (Not entirely sure if the old Xuantie-gcc supports builtin_fma, but it probably does).
I'll try to fix the other platforms before merging your PR, so as not to cause a broken build via the added test)