DGEQRF and DGERQF return wrong optimal dimension LWORK for null column dimension

[andreasvarga]

I came accross this issue working with LAPACK via Julia (see #42762 ). The 7-th argument of the LAPACK subroutines DGEQRF and DGERQF is LWORK, the length of work array WORK.

According to the documentation of DGEQRF, for a normal call for a M x N matrix, LWORK must satisfy LWORK >= max( 1, N ). For an exceptional call with LWORK = -1 and N = 0, the "optimal" dimension LWORK = 0 is returned in WORK(1). Therefore, at the second call, LWORK = 0 is (normally) used, which leads to error -7. To be consistent with the documentation, the returned "optimal" value of LWORK for the call with LWORK = -1 must be at least WORK(1) = 1.

According to the documentation of DGERQF, for a normal call for a M x N matrix, LWORK must satisfy LWORK >= max( 1, M ). For an exceptional call with LWORK = -1 and N = 0, the "optimal" dimension LWORK = 1 is returned in WORK(1), which is a not consistent with the documentation if M > 1. Therefore, at the second call, LWORK = 1 is (normally) used, which leads to error -7 if M > 1. To be consistent with the documentation, the returned "optimal" value of LWORK for the call with LWORK = -1 must be at least WORK(1) = MAX(1,M).

[weslleyspereira]

Hi @andreasvarga. Thanks!

I think the following diff fix the bugs you reported in sgeqrf.f and sgerqf.f.

diff --git a/SRC/sgeqrf.f b/SRC/sgeqrf.f
index 0fc5ba11a..3e79b7d8d 100644
--- a/SRC/sgeqrf.f
+++ b/SRC/sgeqrf.f
@@ -186,7 +186,8 @@
          INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -4
-      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
+      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY .AND.
+     $         M.GT.0 ) THEN
          INFO = -7
       END IF
       IF( INFO.NE.0 ) THEN
diff --git a/SRC/sgerqf.f b/SRC/sgerqf.f
index 24ff42b4c..176b6923d 100644
--- a/SRC/sgerqf.f
+++ b/SRC/sgerqf.f
@@ -176,7 +176,8 @@
          INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -4
-      ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
+      ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY .AND.
+     $         N.GT.0 ) THEN
          INFO = -7
       END IF
 *
@@ -206,6 +207,7 @@
 *     Quick return if possible
 *
       IF( K.EQ.0 ) THEN
+         WORK( 1 ) = 1
          RETURN
       END IF
 *

If you agree, I can make the pull request. If you would like to do it yourself, please go ahead!

[andreasvarga]

I am pleased that you will be able to make a PR for this issue. I would like to recap some facts, in the hope to arrive to the best handling of this issue.

For the special call with LWORK = -1, both routines should return in WORK(1) the "optimal" workspace LWKOPT, to be externally allocated for WORK before a normal call. If MIN(M,N) = 0, this value should be 1. Otherwise, it is should be N*NB for DGEQRF and M*NB for DGERQF.

In the present version, if MIN(M,N) = 0, DGEQRF returns WORK(1) = 0, which is not the desired value and is not correct (see also the documentation), while DGERQF returns WORK(1) = 1, which is the desired value, but not correct for a normal call if M > 1 (see also the documentation). For MIN(M,N) > 0, the returned values are correct.

To fix the error in DGEQRF it is possible to simply replace the 3-rd statement

   `LWKOPT = N*NB`

by

  ` LWKOPT = MAX(1, N*NB)` 

However, this fix is not optimal, becase for M = 0 and N > 0, it still requires to allocate a working place of N*NB, which is not necessary. So LWKOPT should be set according to the following logic: (see also DGERQF)

      IF( MIN(M,N).EQ.0) then
          LWKOPT = 1
      ELSE
         NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
         LWKOPT = N*NB
      END

For the normal call with the "optimal" value LWORK = WORK(1), we assume the above fix was already performed. For DGEQRF, the case M = 0, N > 1 and LWORK = 1 leads to error -7, while for DGERQF, the case M > 1, N = 0 and LWORK = 1 leads to error -7 as well. The fixes you proposed handle these issues corrctly (I believe).

Personally, I would suggest an implementation of the testing logic of input parameters for DGEQRF which is similar to that of DGERQF (with the modification you proposed).

Sorry, if I was not clear enough.

I think the documentation should be also updated.

[VasileSima4]

Hi @andreasvarga. Thanks!

I think the following diff fix the bugs you reported in sgeqrf.f and sgerqf.f.

diff --git a/SRC/sgeqrf.f b/SRC/sgeqrf.f
index 0fc5ba11a..3e79b7d8d 100644
--- a/SRC/sgeqrf.f
+++ b/SRC/sgeqrf.f
@@ -186,7 +186,8 @@
          INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -4
-      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
+      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY .AND.
+     $         M.GT.0 ) THEN
          INFO = -7
       END IF
       IF( INFO.NE.0 ) THEN
diff --git a/SRC/sgerqf.f b/SRC/sgerqf.f
index 24ff42b4c..176b6923d 100644
--- a/SRC/sgerqf.f
+++ b/SRC/sgerqf.f
@@ -176,7 +176,8 @@
          INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -4
-      ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
+      ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY .AND.
+     $         N.GT.0 ) THEN
          INFO = -7
       END IF
 *
@@ -206,6 +207,7 @@
 *     Quick return if possible
 *
       IF( K.EQ.0 ) THEN
+         WORK( 1 ) = 1
          RETURN
       END IF
 *

If you agree, I can make the pull request. If you would like to do it yourself, please go ahead!

Hi,

The fix above seems incomplete, since the optimal LWKOPT is not set for the case N = 0. I include below the full sequence that should replace the existing one in the reference dgeqrf code. This way is easier to check mentally its correctness.

  K = MIN( M, N )
  INFO = 0
  NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  LQUERY = ( LWORK.EQ.-1 )
  IF( M.LT.0 ) THEN
     INFO = -1
  ELSE IF( N.LT.0 ) THEN
     INFO = -2
  ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
     INFO = -4
  ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
     INFO = -7
  END IF
  IF( INFO.NE.0 ) THEN
     CALL XERBLA( 'DGEQRF', -INFO )
     RETURN
  ELSE IF( LQUERY ) THEN
     IF( K.EQ.0 ) THEN
        LWKOPT = 1
     ELSE
        LWKOPT = N*NB
     END IF
     WORK( 1 ) = LWKOPT
     RETURN
  END IF

• Quick return if possible
  

• IF( K.EQ.0 ) THEN
     WORK( 1 ) = 1
     RETURN
  END IF
  

Moreover, the documentation should say that LWORK = 1, if N = 0, and LWORK >= N, otherwise.

[andreasvarga]

Moreover, the documentation should say that LWORK = 1, if N = 0, and LWORK >= N, otherwise.

A small correction: the documentation should say that LWORK = 1, if MIN(M,N) = 0, and LWORK >= N, otherwise.

[andreasvarga]

After consulting Vasile Sima, we arrived to the following test, which, we hope, covers all possible calls of DGEQRF. In the above sequence proposed by Vasile Sima, replace

       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN

with

      ELSE IF ( .NOT.LQUERY ) THEN
         IF( LWORK.LE.0 .OR. ( M.GT.0 .AND. LWORK.LT.MAX( 1, N ) ) )  INFO = -7

The documentation should say that LWORK >= 1, if MIN(M,N) = 0, and LWORK >= N, otherwise.
Similar modifications applies to SGEQRF, CGEQRF and ZGEQRF.

Regarding DGERQF, the only change in code should be to replace

         IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
            INFO = -7
         END IF

with

        IF ( .NOT.LQUERY ) THEN
           IF( LWORK.LE.0 .OR. ( N.GT.0 .AND. LWORK.LT.MAX( 1, M ) ) )  INFO = -7
        END IF

The documentation should say that LWORK >= 1, if MIN(M,N) = 0, and LWORK >= M, otherwise.

Note: I just observed that the logic of the code for SGERQF differs, unexpectedly, from that of DGERQF in testing input parameters. This is why, probably for SGERQF another fix is necessary (similar to DGEQRF). For CGERQF and ZGERQF the above fix is sufficient.

I hope we are done with our proposals for these two routines.

[weslleyspereira]

Thank you @andreasvarga and @VasileSima4 for reporting and for the proposal for the fix.

I created JuliaLang/julia#638 that reproduce the changes you mentioned in this thread. Please, feel free to review it or propose changes.