Support negative strides in BLAS.gemv!

[sostock]

Fixes #30765.

This calculates the correct pointers for the BLAS call in the case of vectors with negative stride. I also added an error if stride(A,2) doesn’t meet the requirements.

I can also go through the other BLAS functions to see if some of them need to be fixed as well, but I wanted to get feedback on this PR first.

[ViralBShah]

Do all BLAS implementations handle this correctly?

[sostock]

I tested it with OpenBLAS and MKL (via MKL.jl) on 64-bit Linux, both handle it correctly (at least for the test cases in this PR).

[ViralBShah]

@dkarrasch Does this look good? I like the tests and wondering if we can merge.

[martinholters]

This broke the case where the output length happened to be zero:

julia> BLAS.gemv!('N', 1.0, zeros(0, 10), rand(10), 1., zeros(0))
Float64[] # before
ERROR: `stride(A,2)` must be at least `max(1, size(A,1))` # now

Was that unsafe (e.g. could cause some memory corruption or the like) to do before? Or is the check too strict, i.e. would lda >= size(A,1) suffice?

[sostock]

I haven’t checked whether lda = size(A,1) = 0 causes any problems, but the references (netlib, MKL) state that lda must be at least max(1, m). If we make sure that lda = size(A,1) = 0 works correctly, I’m in favor of relaxing the check to lda >= size(A,1).

[martinholters]

Ah, well, with LBT enabling simple switching of the used BLAS library, it's probably better not to expect anything from it not stated in the specs. Still, it's a bit unfortunate this potentially causes breakage. Alternatively, is it safe to set lda = 1 if size(A,1) == 0? (Looking closer, I think this is what was effectively used before, rather than lda==0.)

[sostock]

Alternatively, is it safe to set lda = 1 if size(A,1) == 0?

I would guess so (I would be surprised if any elements of A are accessed if size(A,1)==0). I will check and prepare a PR.

[martinholters]

If I read the specs correctly, the condition is that lda*n <= length(A), from which it would not be allowed to arbitrarily set lda=1 for n>0 just because m==0. However, the reference code does have an early exit for m==0 which would allow this...

[ViralBShah]

@martinholters if this is an MKL bug (or something they should handle better), even while we need to address as discussed here, let's also file as a separate issue here and I can report to Intel.

[martinholters]

There is no bug in MKL AFAICT. Prior to this change, one could have size(A,1)==size(Y)==0 and gemv! would happily do nothing, whether using MKL or not. And that is also what one would expect when looking at the Fortran code. However, from the docs I think it follows that this is actually not allowed. And this PR replaces the "doesn't do anything" with a "throws an error". So I guess the question boils down to: May we expect other BLAS libraries to behave like the Fortran code, or could they instead require what the original documentation said? As an important data point, MKL seems to be inline with the Fortran code in this regard.

[sostock]

Setting lda = 1 in the case size(A,1) == 0 works with both OpenBLAS and MKL.

However, in the case size(A,2) == 0, both implementations seem to do nothing, which one could consider incorrect:

julia> BLAS.gemv!('N', 1.0, zeros(5,0), zeros(0), 2.0, ones(5)) # non-empty output vector should be scaled by 2
5-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0

Here, beta = 2, so the output vector should be multiplied by 2 instead of doing nothing. Should we error in the case size(A,2) == 0 && size(A,1) > 0? Should we report this as a bug in OpenBLAS/MKL (they both do the same)?

[martinholters]

However, in the case size(A,2) == 0, both implementations seem to do nothing, which one could consider incorrect:

Yes, that's also something I had noticed a while ago and wondered whether the old Fortran code (which also does this) or its accompanying documentation were to be regarded as the reference...

[sostock]

#42012 fixes the case where size(A,1) == 0.

[sostock]

I just noticed that we actually don’t have to error for negative stride(A,2). We can just do something like

if stride(A,2) < 0
    A = @view A[:, end:-1:begin]
    x = @view x[end:-1:begin]
end

Would that be fine? For better performance, one could also calculate the resulting pointer and stride without creating the SubArrays.