Implement S/C/D/Z AXPBY

[grisuthedragon]

Description

This PR implements the AXPBY operation

$$y \leftarrow \alpha x + \beta y$$

which extends the axpy operation by the second scaling factor, just like in gemm or gemv .

This is required to reduce the memory transfers in algorithms like the CG algorithm, where one step is

$$\mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k$$

Until now, this needs to be implemented in one scal and one axpy step. The introduction of the axpby routine allows to read and write p_{k+1} only once from the memory. In other iterative algorithms, like BiCGStab, the subroutine can be used as well.

The routine already exists, for example, in

• https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2024-2/axpby.html
• https://github.com/OpenMathLib/OpenBLAS/blob/develop/interface/axpby.c
• https://www.ibm.com/docs/en/essl/6.2?topic=vss-saxpby-daxpby-scales-vector-vector-by-two-scalars-adds-them-one-another-stores-in-vector

Checklist

• The documentation has been updated.
• Tests for Fortran
• Tests for CBLAS

[langou]

There are some changes related to TFSM (both in SRC and LAPACKE). I do not think this should be part of this commit.

[martin-frbg]

please rebase, it looks like your fork is about two weeks out of date...

[grisuthedragon]

Rebased to current master's state.

[angsch]

Thank you taking on the task of adding another longstanding BLAS extension to reference-BLAS. I really appreciate that.

[angsch]

nit: constanttimes -> constant times

[angsch]

Copy & paste?

[angsch]

Add "implicit none"?

[angsch]

Manual unrolling feels a tad out of time. I would just do the scalar loop. Note that [c,z]axpby do not have the unrolled loop.