The row-major branch transposed the full m-by-n matrix through
temporaries and copied the never-written complementary region of the
destination temporary back over B, corrupting the triangle uplo requires
to be preserved and reading uninitialized memory. Call the Fortran kernel
directly with uplo swapped (U<->L) and m/n swapped, exactly as the
correct column-major path does. Applies identically to s/d/c/z.
Offset computations like j*lda overflow 32-bit lapack_int once
n = lda >= 46341, e.g. LAPACKE_dpotrf crashes in its NaN check before
the factorization runs. Cast to size_t (or int64_t for the tz offsets)
the same way ge/gb/tp/tf already do.
The Golub-Reinsch-SVD-style iteration in xBBCSD used an absolute
convergence threshold THRESH ~ 90*EPS to decide when off-diagonal
bulges are negligible and can be skipped. For rows where the
diagonal entries are small (tiny singular values in the ratio),
this absolute threshold is too large relative to the diagonal,
causing the algorithm to stop chasing bulges prematurely and
producing inaccurate singular vectors (||X21 - U2 D2 V^*|| up to
50*EPS ||X21||).
Fix: change all bulge-convergence checks (RESTART flags in the
inner loop, initial bulge-chase decisions, and the IMAX-1 cleanup)
from absolute to relative by scaling THRESH by the adjacent diagonal
entries:
RESTART11 = |B11E|^2+|BULGE|^2 <= (THRESH * MAX(|B11D(I-1)|,|B11D(I)|,UNFL))^2
This mirrors DBDSQR's relative convergence check |E| <= TOL*|D|.
Applied to all 4 bidiagonal blocks (B11/B21/B12/B22) at all 3
check-points in the iteration.
Fixes#965
The workspace query formula LWKOPT = NW*NB + TSIZE used TSIZE = LDT*NBMAX
= 65*64 = 4160, a hardcoded constant. The blocked algorithm only stores
one T matrix block at a time (reused across loop iterations), so the
per-iteration workspace is LDT*NB, not LDT*NBMAX. For tiny M,N,K where
a single block suffices (NB >= N or NB >= K), LWKOPT was always >= 4160
regardless of problem size.
Fix: change LWKOPT = NW*NB + TSIZE -> LWKOPT = NW*NB + LDT*NB.
The NB adjustment formula when LWORK is limited must consistently use
the per-iteration T storage instead of TSIZE:
NB = (LWORK - TSIZE) / LDWORK -> NB = LWORK / (LDWORK + LDT)
Applied to all 16 routines (s,d,c,z x {orm,unm}{qr,rq,lq,ql}).
Closes#546
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
The row-major path hardcoded `lda_t = MAX(1,k)`, which is only correct
when SIDE='L' (A is K-by-N). When SIDE='R', A is M-by-K and the
transposition dimensions, allocation, and leading-dimension checks were
all wrong.
- Introduce nrowsA/ncolsA/nrowsV determined by SIDE (following the
same pattern as the LAPACKE_?tpmqrt_work fix in PR #540).
- lda_t now uses nrowsA instead of always k.
- lda check compares against ncolsA instead of always m.
- A transpose uses (nrowsA, ncolsA) instead of (k, m).
- LDV check depends on STOREV: ldv >= nrowsV for STOREV='C',
ldv >= k for STOREV='R'.
- V transpose uses nrowsV rows instead of ldv.
Closes#734
Fix DBDSDC/SBDSDC returning non-orthogonal U/V for bidiagonal matrices
with many nearly-equal singular values (e.g., all singular values ≈ 1).
The divide-and-conquer bidiagonal SVD has two code paths:
- Full vector path (DLASD2/SLASD2): deflation tolerance = 8 * EPS
- Compact path (DLASD7/SLASD7): deflation tolerance = 64 * EPS
With the weaker tolerance (8*EPS), singular values differing by only
~10*EPS (e.g., ~2e-15 for double precision) escape deflation. The
subsequent Z computation in DLASD3 then suffers catastrophic
cancellation from denominators (σ_i - σ_j) that are tiny, polluting the
singular vectors and causing loss of orthogonality.
Raise DLASD2 and SLASD2 to 64*EPS, matching DLASD7 and SLASD7, so more
close singular values are deflated and the singular-vector computation
remains stable.
Closes#255
Problem: When NCVT = NRU = NCC = 0 (no singular vectors requested),
BDSQR calls the dqds algorithm. If dqds failed with INFO = 1 or 3,
the routine returned immediately with that error, even though the
standard QR algorithm could still compute the singular values. Only
INFO = 2 triggered the fallback.
Change `IF( INFO .NE. 2 ) RETURN` to `IF( INFO .EQ. 0 ) RETURN`
after the dqds call, so the standard QR fallback runs on any
dqds failure (INFO=1,2,3) and not only INFO=2.
Closes#242
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.
meng004