Aliasing in geadd

[david-cortes]

Fortran has a function logic in which arrays are assumed not to alias with each other, but many BLAS functions like ddot nevertheless still work correctly with aliasing in the const inputs - e.g. one can pass the same array under both x and y:

double result = cblas_ddot(n, x, 1, x, 1);

OpenBLAS recently introduced functions geadd:
https://github.com/OpenMathLib/OpenBLAS/wiki/OpenBLAS-Extensions

which performs matrix addition:

C := alpha*A + beta*B

Compared to functions like gemm, in these each value of the inputs is involved in exactly one value of the outputs, so if it were to be implemented as a simple C loop, it should work correctly even if C is aliased with A or B.

From some quick experiments, it seems to work as expected when B and C refer to the same array, but want to ask nevertheless: is this function guaranteed to produce correct output when there is aliasing in the inputs?

i.e. can it be used to perform an operation like this?

B := alpha*A + beta*B

If so, I think it'd be a nice addition to the documentation, given that function signatures do not use restrict.

[martin-frbg]

As recently as 2015 it seems, back when I was just a hapless user... At first glance, the implementation uses AXPY so I am not sure if aliasing is safe to use on all platforms (given the various optimized assembler routines for that).

[martin-frbg]

This is actually coded as B := alpha *A + beta *B without a standalone C ever coming into play. (Presumably this is how ATLAS did it - I have not checked there yet, and it certainly does not help that IBM's ESSL has another, more elaborate API for ?GEADD that even supports transposition of the matrices)

[david-cortes]

There's a similar one from intel too (omatadd). Their docs do however specify in which cases is aliasing allowed:
https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-dpcpp/2024-0/omatadd.html

In general, A, B, and C must not overlap in memory, with the exception of the following in-place operations:

• A and C may point to the same memory if op(A) is non-transpose and lda = ldc;
• B and C may point to the same memory if op(B) is non-transpose and ldb = ldc.

[martin-frbg]

That reads a lot like ESSL's version of GEADD... and there's already a request to implement this (#4236)