When P=0 and M=Q (U2 is M-by-M, X21 is M-by-M and orthogonal),
or P=M and M=Q (U1 is M-by-M, X11 is M-by-M and orthogonal),
or Q=0, the minimal dimension R = MIN(P, M-P, Q, M-Q) is zero,
so C and S are empty. In these degenerate cases the preceding
code entered the R == P or R == M-P branch and called the
bidiagonalization subroutines (SORBDB2/SORBDB3) which then
accessed elements X21(I,I+1) or X11(I,I+1) out of bounds in
the second DO loop (I = P+1, Q or I = M-P+1, Q).
Fix: handle R = 0 directly in the four driver routines before
the case-dependent bidiagonalization path. The three trivial
CSD cases are:
Q = 0: set U1 = I, U2 = I, return.
P = 0, M = Q: copy X21 to U2, set V1T = I, return.
P = M, M = Q: copy X11 to U1, set V1T = I, return.
All other R=0 configurations are non-standard and still fall
through to the existing path (which may fail).
Fixes#549
The workspace NB queries in all four ?gelsy variants used
?GERQF and ?UNMRQ, but the routines actually called at runtime
are ?TZRZF and ?UNMRZ (real) / ?UNMRZ (complex).
While ?GERQF and ?TZRZF are related (both produce an RQ-like
factorization), their block sizes may differ in ILAENV's tuned
values. If the user relied on LWKOPT from the query to allocate
exactly the right workspace, the actual ?TZRZF/?UNMRZ calls could
underflow or overflow the work array.
Fix the ILAENV calls to match the routines actually invoked:
?GERQF -> ?TZRZF
?UNMRQ -> ?UNMRZ
Closes#676
Per LAPACK docs, in [cz]tgexc, ilst is [in,out]. The LAPACKE C
interface was passing both ifst and ilst by value as input-only
scalars, so callers could never observe the updated ilst value.
Change both ifst and ilst to lapack_int* pointers in the complex
tgexc variants to match the real-precision stgexc/dtgexc interface
and the underlying Fortran semantics.
Fixes#771
When ?HSEQR encounters non-finite input, it can return a negative
INFO. The undo-scaling block at label 50 then computes N-INFO (which
exceeds N) and indexes WR( INFO+1 ) (out of bounds), causing memory
corruption through ?LASCL.
Add INFO.GE.0 to the outer IF( SCALEA ) guard so the undo-scaling
block is skipped entirely when INFO is negative. The inner
IF( INFO.GT.0 ) guard only covered the second pair of ?LASCL calls.
Fixes#1128
When DFLAG=0 and both DD1 and DD2 need rescaling, the first
scaling loop transitions DFLAG from 0 to -1 and correctly
scales DH11 and DH12. But the second scaling loop then hits
the ELSE branch (matching DFLAG=-1) which unconditionally
resets DH21=-1 and DH12=1, overwriting the scaled DH12.
Change ELSE to ELSE IF (DFLAG.EQ.ONE) so the implied-element
initialization (DH21=-1, DH12=1) only fires when DFLAG=1,
not when DFLAG=-1 (where all elements are already explicit).
Fixes#244
Two comments in the DGEHRD/DHSEQR test incorrectly say "Compute
Schur form":
- The DGEHRD call reduces to upper Hessenberg form, not Schur.
- The DHSEQR call computes eigenvalues and the Schur form, so
use a more descriptive label.
Fixes#587
When JOBU='N' or JOBVT='N', the U and VT matrices are not referenced,
so the leading dimension checks should be skipped. Previously the
code rejected ldvt=0 (and ldu=0) even when those matrices were not
used, causing a spurious INFO=12 error on row-major calls with
JOBVT='N' (or INFO=10 for JOBU='N').
Fixes#1090
The actual SLAQZ0 call uses 'S' when ILV is true and 'E' otherwise,
but the query gated that choice on ILVL, undersizing LWORK for
JOBVL='N', JOBVR='V'. Mirror DGGEV3: gate SLAQZ0 and SGGHD3 on ILV,
SORGQR on ILVL.
DLASWLQ tested NB.LT.0 instead of NB.LE.0, accepting NB=0 silently and
falling into the divide-by-(NB-M) path on line 255 (KK = MOD((N-M),(NB-M)))
when M=0. The single/complex/double-complex siblings all use NB.LE.0.
The argument-validation block at the top of ?LAQZ0 sets INFO=-i and
correctly calls XERBLA(name, -INFO), but the second validation block
(LWORK .LT. LWORKREQ) calls XERBLA(name, INFO) with the negative
value. XERBLA prints its argument verbatim, so this would emit
"parameter number -18 had an illegal value" instead of "18".
Affects S/D/C/Z LAQZ0; ?LAQZ2 and the rest of LAPACK already use
the correct -INFO form.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
INFO=-N value reported to XERBLA on a too-small workspace
(or LRWORK) does not match the actual signature position of the
argument being tested.
Sites:
- ?ORBDB4 / ?UNBDB4: LWORK test reported -14 (WORK), should be -15.
- ?GGSVD3: LWORK.LT.1 test reported -24 (the INFO arg in S/D;
IWORK in C/Z), should be -22 (LWORK).
- ?ORCSD / ?UNCSD: LWORK test reported -22 (LDU2), should be -28
(LWORK). Z/C also report LRWORK as -24 (LDV1T), should be -30.
- C/ZLAQZ0: LWORK test reported -19 (RWORK), should be -18.
(S/DLAQZ0 use ALPHAR+ALPHAI instead of ALPHA, so LWORK is at
position 19 there and -19 is already correct.)
- C/ZLAQZ2: LWORK test reported -26 (RWORK), should be -25.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
SGGES body calls XERBLA at line 423 but the EXTERNAL block omits it;
DGGES correctly lists XERBLA. Add it.
ZRSCL's Doxygen header advertised the routine as ZDRSCL (which is a
different routine that ZRSCL calls internally). Fix the brief, the
download links, and the filename references so the generated docs
correspond to ZRSCL.
- dlaqp2.f:171: replace stale DLARF with DLARF1F
(body calls DLARF1F; slaqp2 already declares SLARF1F).
- zhetrf_aa.f:166: drop stale ZGEMV
(body calls only ZGEMM; chetrf_aa correctly omits CGEMV).
- zlahef_aa.f:173: drop stale ZGEMM
(body calls only ZGEMV; clahef_aa correctly omits CGEMM).
- slaqr2.f:314: replace stale SLARF1L with SLARF1F
(body calls SLARF1F).
- sgelqt.f:147: drop stale SGEQRT2 and SGEQRT3
(body calls only SGELQT3 and SLARFB).
In the LSVEC && RSVEC branch of SGEJSV, the IF arm fills the strict
upper triangle of U with U(p,q) = -SIGN(TEMP1, U(q,p)) while the
ELSE arm zeros it via SLASET('U', ...). Both arms then call SGESVJ
on U with JOBA='L', which tells SGESVJ that U is lower triangular
and the strict upper is zero. After the IF arm that is not true:
SGESVJ silently drops the upper-triangular data the IF arm just
deposited, biasing the SVD result. JOBA='G' (general) is correct
for both arms; DGEJSV already uses 'G' at the matching call site.
In the RSVEC && !LSVEC && !ALMORT branch of DGEJSV, the inner DGESVJ
call passes WORK(N+1) as the workspace base (offset by N) but the
full LWORK as the workspace size. DGESVJ may then write up to N
entries past the end of WORK(LWORK). The float sibling SGEJSV
correctly passes LWORK-N at the same call site.
Bug: In ?ROT call LDX11 is passed as the leading dimension of X21.
Fix: use LDX21 as the INCY argument, matching how X21 is addressed
everywhere else in the same loop body.
The extended complex NRM2 tests can compare an expected +Inf norm with
a computed +Inf norm. The relative-error path then evaluates Inf - Inf,
producing NaN and reporting a spurious failure.
Treat exactly equal results as zero test ratio before computing the
relative difference. This matches the existing SNRM2 and DNRM2 test
logic.
Fixes#1047.
langou