GEQRF: document corner case of LWORK

[amigalemming]

I expected that a workspace query would always return a valid LWORK. This is not the case, a workspace query can yield LWORK=0. Should a warning be added to the documentation of the workspace query?

[ilayn]

Can you also add which query leads to LWORK=0?

[amigalemming]

E.g. M=1, N=0.

[ilayn]

If N=0 then it will quick exit and WORK will not be addressed so I think that is not a problem.

[amigalemming]

But INFO is set to -7.

[ilayn]

Ah that can be too. Netlib site is down so I can't look at the source code but either it will raise an error or quick exit due to N=0, depends on which one is checked first. In both cases WORK is not addressed if I remember correctly.

[amigalemming]

Even on quick return the first element of WORK is written: https://github.com/Reference-LAPACK/lapack/blob/master/SRC/sgeqrf.f

[ilayn]

There are two cases in the first one LWORK=-1 and N = 0. In this case, WORK(1) will be set to 0 and calling SGEQRF with this N=0 and LWORK = 0 will lead to INFO=-7. That is the expected outcome since LWORK < MAX(1, N)

In the second case assume we used an arbitrary number X > 1 but N=0 still then in the following block

*
*     Quick return if possible
*
      K = MIN( M, N )
      IF( K.EQ.0 ) THEN
         WORK( 1 ) = 1
         RETURN
     END IF

This will cause the quick exit. That will be the result because A has a shape (LDA, 0) hence nothing is needed to be modified and TAU has the shape 0 again nothing is needed. Initially work size is unspecified ( * ) but its specification is that it will have the shape (MAX(1,LWORK)) hence first element is written so I don't see the problem of being corner case.

[mgates3]

The problem seems to be if a user did something like this:

    // query for lwork
    double dummy[1];
    info = dgeqrf( m, n, ..., dummy, -1 );
    // allocate work
    lwork = dummy[0];
    work = new double[ lwork ];
    // call QR
    info = dgeqrf( m, n, ..., work, lwork );
    delete[] work;
    // handle errors if info != 0

the code will work except when n = 0, when it will yield info = -7. Since there's nothing to do, it would be nicer if that code finished with info = 0. Either the user should set

    lwork = max( 1, dummy[0] );

or LAPACK should set:

    LWKOPT = MAX( 1, N*NB )

instead of

    LWKOPT = N*NB

so that the returned "optimal size of the WORK array" is a legal value. This would potentially affect a large number of LAPACK routines — geqrf, gelqf, gebrd, hetrd, gehrd, to name just a few.

[weslleyspereira]

I believe this issue was fixed by #638. Maybe we want to let this open to remember that the problem still persists for several other routines, e.g., gelqf, gebrd, hetrd and gehrd.

[dklyuchinskiy]

I believe this issue was fixed by #638. Maybe we want to let this open to remember that the problem still persists for several other routines, e.g., gelqf, gebrd, hetrd and gehrd.

This is actual, indeed. I have recently met this problem during gebrd/sytrd usage inside other computational routine. In this case you need to remember, that lwork can be 0 for some corner cases, which will lead to memory problems. And user should always set lwork = max( 1, dummy[0] ); which is inconvenient.