LoongArch64: fixed cscal and zscal

[XiWeiGu]

For the parameters float x[2] = {NaN, NaN} and float alpha[2] = {0.0, 0.0}, the optimized cscal interface does not directly copy 0.0 to x but continues performing complex multiplication, resulting in an output of {NaN, NaN}.
The optimized zscal has the same issue. This problem was detected in LAPACK tests, but the existing OpenBLAS test cases do not cover this scenario. It may be considered for inclusion in future test cases.

[martin-frbg]

I wonder if this will lead us down the same path of adding a special flag for array zeroing vs IEEE compliance as with non-complex SCAL :(

[XiWeiGu]

You reminded me that we also need to add flags for cscal and zscal.
For the following code:

#include <stdio.h>
#include "cblas.h"
#include <complex.h>
int main() {
    // 向量长度
    int N = 1;

    // 复数缩放因子 alpha
    float alpha[2] = {0.0, 0.0}; 

    // 复数向量 X
    float X[2] = {NAN, NAN};

    // 打印缩放前的向量
    printf("Before scaling:\n");
    for (int i = 0; i < N; i++) {
        printf("X[%d] = %f + %f\n", i, X[i], X[i + 1]);
    }

    // 缩放向量
    cblas_cscal(N, alpha, X, 1);

    // 打印缩放后的向量
    printf("\nAfter scaling:\n");
    for (int i = 0; i < N; i++) {
        printf("X[%d] = %f + %f\n", i, X[i], X[i + 1]);
    }

    return 0;
}

Test output using MKL 2024.2 version：

Before scaling:
X[0] = nan + nan

After scaling:
X[0] = nan + nan

The same output as MKL when using reference BLAS Version 3.12.0.
Test output using OpenBLAS V0.3.29(TARGET=HASWELL)：

Before scaling:
X[0] = nan + nan

After scaling:
X[0] = 0.000000 + 0.000000

[XiWeiGu]

This PR will introduce new issues to s/zscal and needs to be revised. (It seems that other platforms also need modifications to avoid the above issues.)

[XiWeiGu]

I submitted a PR #5081 attempting to fix the implementation in C.

[martin-frbg]

Thanks, I'll try to take a stab at the other implementations over the weekend.

[martin-frbg]

Unfortunately this appears to have broken most of the pre-existing special handling of 0*NAN and 0*INF cases (which mostly covered only one of the real and imaginary components being "special") ?

[XiWeiGu]

I removed some special-case handling code because I felt its correctness was questionable.
Test the following code using MKL 2024.02:

#include <stdio.h>
#include <mkl_cblas.h>
#include <complex.h>
#include <string.h>

int main() {
    int i;
    int N = 17;
    int incX = 1;
    int incY = 1;
    float alpha[2] = {0.0, 0.0};
    float beta[2] = {2.0, 0.0};
    char  trans = 'N';
    float A[17 * 17 * 4];
    float X[17 * 2];
    float Y[17 * 2];

    memset(A, 0, sizeof(A));
    memset(X, 0, sizeof(X));
    for (i = 0; i < (2 * N - 2); i += 4)
    {
        Y[i]     = NAN;
        Y[i + 1] = 1.0;

        Y[i + 2] = INFINITY;
        Y[i + 3] = 1.0;
    }
    Y[2 * N - 2] = NAN;
    Y[2 * N - 1] = 1.0;
    cblas_cgemv(CblasColMajor, CblasNoTrans, N, N, alpha, A, N, X, incX, beta, Y, incY);
    for (i = 0; i < 2 * N; i += 2) {
            printf("%f, %f, ", Y[i], Y[i+1]);
    }
}

The output is:
nan, nan, inf, -nan, nan, nan, inf, -nan, nan, nan, inf, -nan, nan, nan, inf, -nan, nan, nan, inf, -nan, nan, nan, inf, -nan, nan, nan, inf, -nan, nan, nan, inf, -nan, nan, nan
The same output as MKL when using reference BLAS Version 3.12.0.
However, when using OpenBLAS v0.3.29 (TARGET=GENERIC, using the C kernel), the output result is:
nan, 2.000000, inf, 2.000000, nan, 2.000000, inf, 2.000000, nan, 2.000000, inf, 2.000000, nan, 2.000000, inf, 2.000000, nan, 2.000000, inf, 2.000000, nan, 2.000000, inf, 2.000000, nan, 2.000000, inf, 2.000000, nan, 2.000000, inf, 2.000000, nan, 2.000000
I think for the upper-level interfaces that call c/zscal_k (excluding c/zscal), we should also no longer reduce computations.

[XiWeiGu]

OpenBLAS/kernel/arm/zscal.c

Lines 61 to 67 in 76db346

	{
	temp = - da_i * x[ip+1] ;
	if (isnan(x[ip]) || isinf(x[ip])) temp = NAN;
	if (!isinf(x[ip+1]))
	x[ip+1] = da_i * x[ip] ;
	else x[ip+1] = NAN;
	}

Adding special value checks for each number is likely to cause a performance drop and make optimization with assembly much more difficult. Should we consider simplifying it?