Keep the divide and conquer workspace of cLAED0 and zLAED0 out of the caller's Z

[rmlarsen]

Disclaimer: The initial changes were prepared using Claude Code; review and subsequent test fixes used Codex.

Summary

cLAED0 and zLAED0 solve the leaf problems of the divide and conquer recursion into QSTORE, the N by N workspace with leading dimension LDQS, and then merge them with xLAED7 while passing the caller's Q, which has leading dimension LDQ, as the QSIZ*N packed workspace of xLAED7. xLAED7 and xLAED8 address that workspace with leading dimension QSIZ, so whenever LDQ > QSIZ the merge writes into rows QSIZ+1 to LDQ of Q, which belong to the caller. For cSTEDC and zSTEDC with COMPZ = 'V' and N > SMLSIZ that is Z(N+1:LDZ, :) whenever Z is stored in a larger array, for any matrix: the memory is inside the array the caller passed, so nothing crashes, but the rows below N come back overwritten. This PR keeps the eigenvectors in Q with its own leading dimension and gives QSTORE, which is contiguous, to xLAED7 as its workspace. The numerical changes are in {c,z}laed0.f, with regression coverage in {c,z}chkst.f. The real routines keep their merge workspace in WORK and are not affected.

Description

*     I am free to use Q as a valuable working space until Loop 150.
            CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
     $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
     $                   ...
     $                   Q( 1, SUBMAT ), RWORK( IWREM ),
     $                   IWORK( SUBPBS+1 ), INFO )

The Q(1,SUBMAT) actual argument is xLAED7's WORK, documented as "dimension (QSIZ*N)". xLAED7 passes it to xLAED8 as Q2 with LDQ2 = QSIZ and to xLACRM as a QSIZ by K matrix with leading dimension QSIZ, so the workspace is laid out as if Q had leading dimension QSIZ; with LDQ > QSIZ column j of the workspace starts at offset (j-1)*QSIZ of Q, which is inside the padding of an earlier column. The comment is right that Q is free at that point, since the leaf products have moved its contents into QSTORE; the leading dimension was not accounted for.

Fix. After each leaf product xLACRM into QSTORE, the block is copied back into Q( 1, SUBMAT ) with LDQ, so the eigenvectors accumulate in Q; xLAED7 gets Q( 1, SUBMAT ) with LDQ as its eigenvector matrix and QSTORE as its workspace, which holds LDQS*N >= QSIZ*MATSIZ contiguous entries. The final permutation of the columns goes through QSTORE and back into Q. This costs one copy of QSIZ by MATSIZ per leaf and one of QSIZ by N at the end, O(N^2) against the O(N^3) merges, and changes no arithmetic: the eigenvalues and eigenvectors are bit-identical to the parent commit. The documentation of QSTORE and of LDQS >= max(1,QSIZ), which the leaf products already required, is updated with it.

Minimal reproducer

! ZSTEDC('V') on a 26x26 tridiagonal matrix with Z stored in the first 26 rows
! of a 30-row array: rows 27 to 30 belong to the caller.
program minimal
  implicit none
  integer, parameter :: n = 26, ldz = 30
  double precision :: d(n), e(n), rwork(6*n*n)
  complex(8) :: z(ldz,n), work(n*n), marker
  integer :: iwork(20*n*n), info, i, j
  marker = (-7d0, 7d0)
  do i = 1, n
    d(i) = i
    e(i) = 1d0 / (i + 1)
  end do
  z = marker
  do j = 1, n
    z(1:n, j) = (0d0, 0d0)
    z(j, j) = (1d0, 0d0)
  end do
  call zstedc('V', n, d, e, z, ldz, work, size(work), rwork, size(rwork), iwork, size(iwork), info)
  print '(a,i0,a,i0,a,i0)', 'info = ', info, '  entries of Z(27:30, :) overwritten = ', &
        count(z(n+1:ldz, :) /= marker), ' of ', (ldz-n)*n
end program

BEFORE (master):     info = 0  entries of Z(27:30, :) overwritten = 88 of 104
AFTER (this branch): info = 0  entries of Z(27:30, :) overwritten = 0 of 104

Regression test. cCHKST and zCHKST call xSTEDC('V') on a tridiagonal matrix of order N = max(2, SMLSIZ+1) (26 by default). The case owns its diagonal and off-diagonal arrays, eigenvalues, and a (N+1) by N eigenvector matrix. The extra row contains a marker. This reaches divide and conquer independently of the input sizes and respects the caller's array capacities. An overwritten marker, nonzero INFO, or excessive residual/orthogonality ratio from xSTT21 is a failure. The parent overwrites 25 markers in both precisions; the fix preserves all markers and passes both ratios.

Validation

• Current tests: gfortran 13.3 on x86-64, reference BLAS, -O2 -fcheck=all; all 4 focused sep.in/se2.in driver runs pass.
• 16 AddressSanitizer cases pass with small caller arrays and compact/padded leading dimensions. The storage checks also use SMLSIZ = 25 and 40, and verify that the dedicated regression leaves caller matrix and eigenvalue arrays untouched.
• The updated regression still fails against the parent numerical kernels in every affected precision.
• The supplied minimal reproducer was independently checked against the parent and the numerical fix and matches the output above. Additional rectangular ZLAED0 checks at N=26,51,100,131 preserve all padding with the same normalized residuals as the parent.
• The full suite and cross-architecture/compiler matrix have not been rerun for this test-only follow-up. Earlier validation of the unchanged numerical kernels included the full LAPACK suite and the special-value/finite sweeps described in the original submission.

Found while auditing the symmetric tridiagonal eigensolvers for the NaN and overflow handling of #1377-#1391: the sweep's canaries caught the write for a NaN input, and the reproducer shows it is independent of the input.

Update: Codex review identified the regression test's reliance on LDU as array capacity. Dedicated arrays now preserve the padding, residual, and orthogonality checks without accessing caller storage beyond its documented bounds.

[codecov]

Codecov Report

❌ Patch coverage is 80.95238% with 16 lines in your changes missing coverage. Please review.
✅ Project coverage is 69.57%. Comparing base (a6c6e74) to head (973813c).
✅ All tests successful. No failed tests found.

Files with missing lines	Patch %	Lines
TESTING/EIG/cchkst.f	79.48%	8 Missing ⚠️
TESTING/EIG/zchkst.f	79.48%	8 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1406      +/-   ##
==========================================
+ Coverage   69.36%   69.57%   +0.21%     
==========================================
  Files        6122     6122              
  Lines      486711   486766      +55     
  Branches    23268    23268              
==========================================
+ Hits       337584   338680    +1096     
+ Misses     148689   147648    -1041     
  Partials      438      438              

Components	Coverage Δ
BLAS	97.94% <ø> (ø)
CBLAS	96.98% <ø> (ø)
LAPACK	82.94% <100.00%> (+0.55%)	⬆️
LAPACKE	2.17% <ø> (ø)
TMGLIB	55.69% <ø> (ø)
BLAS testing	88.32% <ø> (-0.02%)	⬇️
CBLAS testing	89.63% <ø> (ø)
LAPACK testing	82.20% <79.48%> (-0.01%)	⬇️
LAPACKE testing	∅ <ø> (∅)

Files with missing lines	Coverage Δ
SRC/claed0.f	84.90% <100.00%> (+84.90%)	⬆️
SRC/zlaed0.f	84.90% <100.00%> (+84.90%)	⬆️
TESTING/EIG/cchkst.f	64.21% <79.48%> (+1.06%)	⬆️
TESTING/EIG/zchkst.f	64.21% <79.48%> (+1.06%)	⬆️

... and 13 files with indirect coverage changes

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a6c6e74...973813c. Read the comment docs.

[ACSimon33]

Hi @rmlarsen,
could you please merge the current master back into all your PRs? We fixed numerical issues and errors on all platforms that are tested in CI, and the CI jobs are now configured to fail if new errors occur. That will help us to quickly determine if the changes in the PRs introduce new errors or can be merged.

[rmlarsen]

Hi @rmlarsen, could you please merge the current master back into all your PRs? We fixed numerical issues and errors on all platforms that are tested in CI, and the CI jobs are now configured to fail if new errors occur. That will help us to quickly determine if the changes in the PRs introduce new errors or can be merged.

@ACSimon33 done!