lapack
lapack copied to clipboard
fix workspace size in tgsen
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.
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