Y is misscaled in ?GEMV for some cases

[bE554357]

Description

Description and possible solution is presented in #787

[jgpallero]

Hi,

there are also affected the subroutines xGBMV, xSBMV, xSPMV, and xSYMV

[langou]

Hi @jgpallero and @bE554357, @akobotov also gave a thumbs up,

I think it is best to keep the BLAS behavior as it is. I consulted with a few other folks by email and (at least) three other persons agree that it is best to keep the BLAS behavior as it is. Arguments being: "It's essentially impossible to change those documented semantics of BLAS gemv, since any change would have to be agreed to and propagated through all the vendor libraries" and this is a "classical and well-documented behavior".

For now, I am closing the comment and the PR. I do see that there is a disagreement on what to do since there were three initial thumbs up on the initial PR.

All in all, it is great that you are willing to change the library and propose a PR for it. But I think, for this change, it is best to leave the code as is.

cheers,
Julien.

[jgpallero]

Thank you for your answer. Maybe the "wrong" part in this case is not the behavior of the subroutines, but the documentation. In all of them the docs say referred to the matrix dimensions "must be at least zero", also the doc for MKL and ESSL BLAS. The question is then how behave MKL and ESSL in this situation, but I suppose in the same way as the reference BLAS. I've inspected the code in MAGMA and it works in the same way, i.e., if (m <= 0 || n <= 0) return;. OpenBLAS present the same behavior, as well as ATLAS.

[OFFTOPIC]
And talking about ATLAS, anyone knows what is the current status of the project? The last commit in https://github.com/math-atlas/math-atlas was in 2019

[langou]

My understanding was that ATLAS is not active any longer and is now moved to OpenBLAS. But I am not 100% sure. I am surprised that there was a commit in 2019 for ATLAS. I was assuming that the project was long not maintained any more. So I think you want to look at OpenBLAS.

[langou]

The question is then how behave MKL and ESSL in this situation, but I suppose in the same way as the reference BLAS. I've inspected the code in MAGMA and it works in the same way, i.e., if (m <= 0 || n <= 0) return;. OpenBLAS present the same behavior, as well as ATLAS.

Yes. I want to point to: https://github.com/tlapack/testBLAS. In this software, @weslleyspereira is trying to check BLAS libraries. These corner cases are in. Yes, most current libraries are following the convention as set by the reference BLAS.

[mgates3]

Note that gemm has different semantics: if k = 0, it still scales the C matrix. This isn't in the manpage, but it is in the 1990 paper (https://doi.org/10.1145/77626.79170) and reference code, and appears to work with OpenBLAS, MKL, and Accelerate. So a user could use

    gemm( NoTrans, NoTrans, m, 1, 0, ... )

instead of

    gemv( m, 0, ... )

[mgates3]

Testing several BLAS libraries, I found gemv in OpenBLAS and MKL did not scale y when n = 0. With macOS Accelerate, the Fortran API gemv did not scale y, but cblas_Xgemv did scale y (inconsistent with documented behavior).

blaspp/test> ./tester --dim 3x0 --verbose 2 --alpha 2 --beta 4 gemv
...
y    = [    0.0000    0.1315    0.7556 ]';
yout = [    0.0000    0.1315    0.7556 ]';    // BLAS++ gemv calling Fortran gemv
yref = [    0.0000    0.5262    3.0224 ]';    // cblas_gemv

whereas with MKL:

sh leconte test> ./tester --dim 3x0 --verbose 2 --alpha 2 --beta 4 gemv
...
y    = [    0.8402    0.3944    0.7831 ]';
yout = [    0.8402    0.3944    0.7831 ]';
yref = [    0.8402    0.3944    0.7831 ]';

[langou]

Maybe the "wrong" part in this case is not the behavior of the subroutines, but the documentation. In all of them the docs say referred to the matrix dimensions "must be at least zero", also the doc for MKL and ESSL BLAS.

Ah. Thanks. Yes, the documentation could be fixed. Good point.

[martin-frbg]

I'd read that as "must be at least zero to avoid raising an error " ?

[jgpallero]

Maybe the "wrong" part in this case is not the behavior of the subroutines, but the documentation. In all of them the docs say referred to the matrix dimensions "must be at least zero", also the doc for MKL and ESSL BLAS.

Ah. Thanks. Yes, the documentation could be fixed. Good point.

Maybe a note with something like "if M or N are 0, the function does nothing"

[jgpallero]

I'd read that as "must be at least zero to avoid raising an error " ?

I'm sure this was the original idea

[langou]

I agree with @martin-frbg. I think the comments are fine then.

Yes, we can add If M or N are 0, the function performs a quick return.

[bE554357]

Now it's hard to localize error in gemv while debugging. Is it possible to call xerbla for cases of $\mathbf{y}$ misscale?