lapack icon indicating copy to clipboard operation
lapack copied to clipboard

fix workspace size in tgsen

Open mgates3 opened this issue 2 years ago • 1 comments

Description The workspace query & docs for tgsen appear to be off by 1. It claims to need 2*M*(N-M), but then calls tgsyl with LWORK-2*N1*N2 where N1 = M, N2 = N - M, so the tgsyl lwork = 0: https://github.com/Reference-LAPACK/lapack/blob/master/SRC/ctgsen.f#L628 Thus tgsyl raises an error about invalid lwork (parameter 20). Adding 1 solves the problem.

Discovered when doing LAPACK++ wrapper. To see the bug, disable the workaround in LAPACK++ tgsen:

diff --git a/src/tgsen.cc b/src/tgsen.cc
index 830ca4f..f5f57d5 100644
--- a/src/tgsen.cc
+++ b/src/tgsen.cc
@@ -211,7 +211,7 @@ int64_t tgsen(
         throw Error();
     }
     // LAPACK <= 3.11 has query & documentation error in workspace size; add 1.
-    lapack_int lwork_  = real( qry_work[ 0 ] ) + 1;
+    lapack_int lwork_  = real( qry_work[ 0 ] ); // + 1;

Before this PR:

lapackpp/test> ./tester --ijob 1 --dim 100 --type s,d,c,z tgsen
LAPACK++ version 2022.07.00, id a9361a2
input: ./tester --ijob 1 --dim 100 --type 's,d,c,z' tgsen
                                                                              
type   ijob  jobvl  jobvr       n     error   time (s)  ref time (s)  status  
 ** On entry to STGSYL parameter number 20 had an illegal value
   s      1  novec  novec     100  8.22e-02   0.000734      0.000996  FAILED  

 ** On entry to DTGSYL parameter number 20 had an illegal value
   d      1  novec  novec     100  8.66e-02   0.000581      0.000801  FAILED  

 ** On entry to CTGSYL parameter number 20 had an illegal value
   c      1  novec  novec     100  1.55e-01   0.000239      0.000457  FAILED  

 ** On entry to ZTGSYL parameter number 20 had an illegal value
   z      1  novec  novec     100  1.46e-01   0.000341      0.000471  FAILED  
4 tests FAILED for tgsen.

After this PR:

lapackpp/test> ./tester --ijob 1 --dim 100 --type s,d,c,z tgsen
LAPACK++ version 2022.07.00, id a9361a2
input: ./tester --ijob 1 --dim 100 --type 's,d,c,z' tgsen
                                                                              
type   ijob  jobvl  jobvr       n     error   time (s)  ref time (s)  status  
   s      1  novec  novec     100  0.00e+00   0.000916      0.000841  pass    

   d      1  novec  novec     100  0.00e+00   0.000829      0.000942  pass    

   c      1  novec  novec     100  0.00e+00   0.000377      0.000429  pass    

   z      1  novec  novec     100  0.00e+00   0.000470      0.000460  pass    
All tests passed for tgsen.

Routines like gges that call tgsen in most cases don't use a workspace query, and presumably have a large enough lwork that this issue never arose.

Checklist

  • [x] The documentation has been updated.
  • [ ] If the PR solves a specific issue, it is set to be closed on merge.

mgates3 avatar Dec 02 '22 19:12 mgates3

Good catch, I think this change also needs to propagated to the calling routines as they don't always do workspace queries. I.e. https://github.com/Reference-LAPACK/lapack/blob/master/SRC/cggesx.f#L636

thijssteel avatar Mar 14 '23 10:03 thijssteel