DGEQRF and DGERQF return wrong optimal dimension LWORK for null column dimension #633
Comments
andreasvarga
changed the title
|
Hi @andreasvarga. Thanks! I think the following diff fix the bugs you reported in 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! |
|
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 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 by 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) 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. |
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.
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. |
|
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 with The documentation should say that LWORK >= 1, if MIN(M,N) = 0, and LWORK >= N, otherwise. Regarding DGERQF, the only change in code should be to replace with 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. |
|
Thank you @andreasvarga and @VasileSima4 for reporting and for the proposal for the fix. I created #638 that reproduce the changes you mentioned in this thread. Please, feel free to review it or propose changes. |
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).
The text was updated successfully, but these errors were encountered: