#694625 libblas3gf: DGEMV does not handle properly N=0

Package:
libblas3
Source:
lapack
Description:
Basic Linear Algebra Reference implementations, shared library
Submitter:
Christophe TROESTLER
Date:
2015-12-23 23:03:12 UTC
Severity:
normal
#694625#5
Date:
2012-11-28 14:05:45 UTC
From:
To:
Dear Maintainer,

The routine DGEMV (with TRANS = 'N') is supposed to perform the
operation

    y ← beta * y + alpha * A * x

where the matrix A has dimensions M×N.  In particular, if N=0 (and M >
0), A * x is the null vector (all elements are given by sums on an
empty set of indices).  Consequently, if N = 0, the operation should
be

    y ← beta * y

and NOT leave y untouched.  For your convenience, I have attached a
program demonstrating the latter (wrong) behavior.
#694625#10
Date:
2012-11-28 14:20:33 UTC
From:
To:
severity 694625 normal
thanks

Bonjour Christophe!
   1.0000000000000000
   2.0000000000000000
What should I get ?
2 & 2 ?

Thanks
Sylvestre
#694625#17
Date:
2012-11-28 14:47:08 UTC
From:
To:
Christophe TROESTLER <Christophe.Troestler@umons.ac.be> writes:

I have tested the three BLAS implementations available in Debian (netlib
BLAS, ATLAS and OpenBLAS), and they all give the same result (bad from
your point of view).

So I am not sure this is a bug. Maybe this is on purpose. If one
considers this operation to be invalid (as I do), then it makes sense to
do nothing.

Does the documentation clearly states that BLAS supports operations on
objects with one dimension equal to zero?
#694625#22
Date:
2012-11-28 14:49:09 UTC
From:
To:
Sylvestre Ledru <sylvestre@debian.org> writes:

Christophe is expecting 2 & 4 (Y=[1,2] and beta=2 in his example)
#694625#27
Date:
2012-11-28 15:05:01 UTC
From:
To:
Salut Sylvestre !

2.
4.

(beta=2 in the code).
#694625#32
Date:
2012-11-28 15:05:09 UTC
From:
To:
If the operation is considered invalid, then it should return an error
instead of doing nothing silently.

However, I insist that the operation is perfectly valid from a
mathematical point of view.  Indeed, in mathematics, ∑_{i ∈ ∅} a_i = 0.
Therefore, A * x = 0 when A has dimensions M×0.  This corner case is
also useful in practice.  For a problem I have, I needed to compute

    t ↦ f(t₀ y + ∑ t_i a_i)

where a_i are vectors and t_i real numbers.  It so happens that in
some circumstances, the sum may be empty leading to a constant
function instead of the correct t ↦ f(t₀ y).  A very bad surprise
IMHO, incoherent with the semantics of the operation.

Not the man pages.  To the contrary, for N, the man page clearly
includes the possibility that N = 0:

       N      - INTEGER.
              On  entry, N specifies the number of columns of the matrix A.  N
              must be at least zero.  Unchanged on exit.
#694625#37
Date:
2012-11-28 15:23:39 UTC
From:
To:
Christophe TROESTLER <Christophe.Troestler@umons.ac.be> writes:

I understand your point, though one could argue that an object of
dimension M×0 is not a matrix. Indeed, Octave for example handles well
these corner cases.

Yes indeed. However it does not tell what it does in that case, though I
agree with you that doing y=beta*y would seem logical.

Do you have examples of other BLAS functions which handle
null-dimensions as you would expect?