(Analysis by Claude)
Summary
cheevd (single-precision complex Hermitian eigensolver, jobz='V') produces
substantially less accurate eigenvectors for uplo='L' than for uplo='U'
on the same input matrix, starting somewhere between OpenBLAS v0.3.33
(good) and v0.3.34 (bad). We bisected this to a single commit:
db20661747657af702011ff2fb94912bf6d2dd63
"Increase deflation tolerance to match recent ?SLASD changes (Reference-LAPACK PR 1317)"
Author: Martin Kroeker
Date: Sun Jul 5 17:37:17 2026 +0200
Reverting only this commit on top of current develop fixes the failure.
Restoring it reproduces the failure again. See Bisection below.
Symptom
Running cheevd(jobz='V', uplo=..., n=100, ...) on a random Hermitian test
matrix (fixed seed, so fully reproducible), then comparing the returned
eigenpairs against the known-truth eigenvalues/eigenvectors used to
construct the matrix (phase-correcting each eigenvector, since A*v = lambda*v is equally satisfied by v' = alpha*v for any unit-modulus
complex alpha):
_EIGEN_VAL_SPACING (gap) = 2.00000009E-03
thresh_evalmax = 9.99999975E-06
thresh_ecolmax = 1.19209290E-03 ( max of floor 5.00000024E-04 and 20.0000000 * eps/gap )
thresh_ealphamax = 1.19209280E-02 ( max of floor 2.00000009E-03 and 200.000000 * eps/gap )
uplo=U evalmax= 9.53674316E-07 ecolmax= 4.86821118E-05 ealphamax= 7.85589218E-05
uplo=U Passed
uplo=L evalmax= 8.34465027E-07 ecolmax= 6.33120025E-03 ealphamax= 2.75581479E-02
uplo=L Failed
uplo='U' and uplo='L' are just two different ways of reading the same
fully-populated Hermitian input matrix (our generator fills both
triangles), so this isn't a difference in test data -- it's cheevd
itself recovering the eigenvectors far less accurately for the
lower-triangular storage/reduction path, on this input.
The input matrix used in our test is from a generator derived from this one, with
an additional logic to enforce a minimum eigenvalue spacing.
The threshold value for the eigenvector verification is based on Davis-Kahan.
See comments in reproducer for details.
Environment
- OpenBLAS
develop @ 632ef8743 (also reproduces at tagged v0.3.34)
- x86_64 (sapphirerapids), CMake build,
USE_THREAD default
- Reproduces identically with two different compilers/toolchains we tried
(nvfortran two different releases, and gfortran 11.5.0), so it isn't
compiler-specific.
v0.3.33 does not reproduce the failure (both uplo cases pass,
~1e-4 accuracy).
Root cause
It is found that the root cause is the change in #5885, confirmed by a revert/rebuild/rerun/restore round-trip:
--- a/lapack-netlib/SRC/slaed2.f
+++ b/lapack-netlib/SRC/slaed2.f
@@ -312,7 +312,7 @@
EPS = SLAMCH( 'Epsilon' )
- TOL = EIGHT*EPS*MAX( ABS( D( JMAX ) ), ABS( Z( IMAX ) ) )
+ TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( JMAX ) ), ABS( Z( IMAX ) ) )
--- a/lapack-netlib/SRC/slaed8.f
+++ b/lapack-netlib/SRC/slaed8.f
@@ -355,7 +355,7 @@
EPS = SLAMCH( 'Epsilon' )
- TOL = EIGHT*EPS*ABS( D( JMAX ) )
+ TOL = EIGHT*EIGHT*EPS*ABS( D( JMAX ) )(commit db2066174 also makes the equivalent change in dlaed2.f/
dlaed8.f for double precision). SLAED2/SLAED8 decide, during
SSTEDC/CSTEDC's divide-and-conquer merge, which eigenvalue pairs are
"close enough" to deflate -- skip solving the secular equation for them
and copy/rotate their eigenvectors directly instead. This commit widens
that tolerance 8x (8*EPS -> 64*EPS). CSTEDC (used by CHEEVD)
and SSTEDC (used by SSYEVD) both route their real-valued secular-
equation merge through this same SLAED2/SLAED8 code, so we'd expect
SSYEVD to be susceptible to the same failure mode, though we haven't
built a matrix that demonstrates it there yet.
Revert/rebuild/rerun round-trip (n=100, iseed=4321, _EIGEN_VAL_SPACING=2.d-3):
| Build |
uplo='U' |
uplo='L' |
develop HEAD (db2066174 present) |
Passed (ecolmax 4.87E-05) |
Failed (ecolmax 6.33E-03, ealphamax 2.76E-02) |
develop HEAD with db2066174 reverted |
Passed (ecolmax 4.82E-05) |
Passed (ecolmax 4.82E-05, ealphamax 2.02E-04) |
develop HEAD, commit restored |
Passed (ecolmax 4.87E-05) |
Failed (ecolmax 6.33E-03, matches row 1 exactly) |
_EIGEN_VAL_SPACING (gap) = 2.00000009E-03
thresh_evalmax = 9.99999975E-06
thresh_ecolmax = 1.19209290E-03 ( max of floor 5.00000024E-04 and 20.0000000 * eps/gap )
thresh_ealphamax = 1.19209280E-02 ( max of floor 2.00000009E-03 and 200.000000 * eps/gap )
uplo=U evalmax= 9.53674316E-07 ecolmax= 4.81969910E-05 ealphamax= 7.83205032E-05
uplo=U Passed
uplo=L evalmax= 8.34465027E-07 ecolmax= 4.81597272E-05 ealphamax= 2.02119350E-04
uplo=L Passed
Sensitivity to eigenvalue spacing
Our matrix generator enforces a floor on the gap between generated
eigenvalues (_EIGEN_VAL_SPACING, build-time constant) to avoid
unrelated single-precision resolution issues with genuinely
near-degenerate eigenvalues. Sweeping this floor against the same
buggy develop build shows the failure is specifically about eigenvalue
proximity, and disappears once the floor is wide enough:
_EIGEN_VAL_SPACING |
uplo='L' |
2.d-3 |
Failed (ecolmax 6.33E-03) |
5.d-3 |
Failed (ecolmax 5.81E-03) |
6.d-3 |
Passed (ecolmax 1.07E-04) |
7.d-3 -- 1.d-1 |
Passed |
This is consistent with the root cause: the new deflation tolerance
(TOL = 64*EPS*max(|D|,|Z|), ~1e-5 for float32 at this matrix's
eigenvalue scale) doesn't act directly on the globally enforced minimum
gap between the eigenvalues we feed in at generation time -- it acts on
gaps that appear locally, between intermediate eigenvalues inside the
divide-and-conquer merge tree (~7 levels for n=100), which can end up
much smaller than the global floor. A 2e-3 global floor still lets some
such local gap fall inside the widened deflation window on the uplo='L'
reduction path for this matrix (but not uplo='U'); a 6e-3+ floor
apparently gives enough margin to avoid it. This is a data point for
diagnosis, not a fix -- widening the test's spacing floor would only
mask the regression, not correct it.
Reproducing
Run the whole thing end to end, from a fresh clone, in one shot, with the attached reproducer reproducer.tgz:
(Analysis by Claude)
Summary
cheevd(single-precision complex Hermitian eigensolver,jobz='V') producessubstantially less accurate eigenvectors for
uplo='L'than foruplo='U'on the same input matrix, starting somewhere between OpenBLAS
v0.3.33(good) and
v0.3.34(bad). We bisected this to a single commit:Reverting only this commit on top of current
developfixes the failure.Restoring it reproduces the failure again. See Bisection below.
Symptom
Running
cheevd(jobz='V', uplo=..., n=100, ...)on a random Hermitian testmatrix (fixed seed, so fully reproducible), then comparing the returned
eigenpairs against the known-truth eigenvalues/eigenvectors used to
construct the matrix (phase-correcting each eigenvector, since
A*v = lambda*vis equally satisfied byv' = alpha*vfor any unit-moduluscomplex
alpha):uplo='U'anduplo='L'are just two different ways of reading the samefully-populated Hermitian input matrix (our generator fills both
triangles), so this isn't a difference in test data -- it's
cheevditself recovering the eigenvectors far less accurately for the
lower-triangular storage/reduction path, on this input.
Environment
develop@632ef8743(also reproduces at taggedv0.3.34)USE_THREADdefault(
nvfortrantwo different releases, andgfortran11.5.0), so it isn'tcompiler-specific.
v0.3.33does not reproduce the failure (bothuplocases pass,~1e-4 accuracy).
Root cause
It is found that the root cause is the change in #5885, confirmed by a revert/rebuild/rerun/restore round-trip:
(commit
db2066174also makes the equivalent change indlaed2.f/dlaed8.ffor double precision).SLAED2/SLAED8decide, duringSSTEDC/CSTEDC's divide-and-conquer merge, which eigenvalue pairs are"close enough" to deflate -- skip solving the secular equation for them
and copy/rotate their eigenvectors directly instead. This commit widens
that tolerance 8x (
8*EPS->64*EPS).CSTEDC(used byCHEEVD)and
SSTEDC(used bySSYEVD) both route their real-valued secular-equation merge through this same
SLAED2/SLAED8code, so we'd expectSSYEVDto be susceptible to the same failure mode, though we haven'tbuilt a matrix that demonstrates it there yet.
Revert/rebuild/rerun round-trip (
n=100,iseed=4321,_EIGEN_VAL_SPACING=2.d-3):uplo='U'uplo='L'developHEAD (db2066174present)developHEAD withdb2066174reverteddevelopHEAD, commit restoredSensitivity to eigenvalue spacing
Our matrix generator enforces a floor on the gap between generated
eigenvalues (
_EIGEN_VAL_SPACING, build-time constant) to avoidunrelated single-precision resolution issues with genuinely
near-degenerate eigenvalues. Sweeping this floor against the same
buggy
developbuild shows the failure is specifically about eigenvalueproximity, and disappears once the floor is wide enough:
_EIGEN_VAL_SPACINGuplo='L'2.d-35.d-36.d-37.d-3--1.d-1This is consistent with the root cause: the new deflation tolerance
(
TOL = 64*EPS*max(|D|,|Z|), ~1e-5 forfloat32at this matrix'seigenvalue scale) doesn't act directly on the globally enforced minimum
gap between the eigenvalues we feed in at generation time -- it acts on
gaps that appear locally, between intermediate eigenvalues inside the
divide-and-conquer merge tree (~7 levels for
n=100), which can end upmuch smaller than the global floor. A
2e-3global floor still lets somesuch local gap fall inside the widened deflation window on the
uplo='L'reduction path for this matrix (but not
uplo='U'); a6e-3+ floorapparently gives enough margin to avoid it. This is a data point for
diagnosis, not a fix -- widening the test's spacing floor would only
mask the regression, not correct it.
Reproducing
Run the whole thing end to end, from a fresh clone, in one shot, with the attached reproducer reproducer.tgz: