gemv not handling N=0 according to math. conventions

[ChriS00]

In order to compute a function t = (t₀, t₁,...) ↦ f(t₀ y + ∑ t_i a_i) where y and a_i are vectors, I have used gemv to perform y ← t₀ y + ∑ t_i a_i. However, in case the sum is empty — ∑_{i ∈ ∅} which translates to N=0 when calling gemv — gemv leaves y unchanged instead of performing the mathematically correct behavior: y ← t₀ y.

Is this broken semantics the intended behavior — or must it be considered as a forgotten corner case? If the former, why? And why hasn't this be fixed — as it seems to me that having the mathematically correct behavior is definitely better and can be obtained at zero additional (computational) cost.

[Julien Langou]

Hi,

This dates back from 1988 and the initial BLAS specification of the Level 2
BLAS. The behavior you are describing (which you do not like) is actually
correct and follows the specification from 1988. Interestingly enough, the
behavior was changed for Level 3 BLAS in 1990 and the Level 3 BLAS follows the
behavior you would have liked to see. See below for historical references and
precise quotes.

To address the issue, you could have an IF statement before calling Level 2
BLAS GEMV or you could use the Level 3 BLAS GEMM instead of the Level 2 BLAS
GEMV.

I understand that you probably have found a fix and made this work and your
comment was higher level than just "I want my code to work". Thanks for the
feedback, you are right, this behavior of Level 2 BLAS GEMV may appear quite
surprising / counter-intuitive indeed. I agree with you.

Cheers,
Julien.

The Level 2 BLAS paper from 1988.

J. J. Dongarra, J. Du Croz, S. Hammarling, and R. J. Hanson, An extended set of
FORTRAN Basic Linear Algebra Subprograms, ACM Trans. Math. Soft., 14 (1988),
pp. 1--17.

page 5 - section 4

The size of the matrix is determined by the arguments M and N for an m-by-n
rectangular matrix [...] Note it is permissible to call the routines with M or
N = 0, in which case the routines exit immediately without referencing their
vector or matrix arguments.

The Level 3 BLAS paper from 1990.

J. J. Dongarra, J. Du Croz, I. S. Duff, and S. Hammarling, A set of Level 3
Basic Linear Algebra Subprograms, ACM Trans. Math. Soft., 16 (1990), pp. 1--17.

page 6 - section 4

The sizes of the matrices are determined by the arguments M, N and K. It is
permissible to call the routines with M or N = 0, in which case the routines
exit immediately without referencing their arguments. If M and N > 0, but K =
0, the operation reduces to C := beta * C [ ... ]

[ChriS00]

Thanks for your answer.