TESTING: scale the xGEEV eigenvalue-consistency check by n ulp |W|

[ACSimon33]

Summary

xDRVEV and xGET23 test 5 require xGEEV to return bit-identical eigenvalues whether or not eigenvectors were asked for, setting the ratio to 1/ulp on any difference at all. That is stronger than anything the routine promises: the two JOB values run different code paths, and nothing obliges those to round alike. This PR measures the difference against n |W| ulp instead, the same normalisation tests 1 and 2 of the same drivers already use.

Closes #732.

Why bit equality is not available

xGEEV with JOBVL = JOBVR = 'N' calls xHSEQR with JOB = 'E', which is free to skip work the Schur form would need. In xLAHQR, xLAQR3 and xLAQR5 that shows up as narrower update ranges - I1, I2, JTOP and JBOT all depend on WANTT - so the two JOB values reach the same eigenvalues through a different sequence of BLAS calls. Nothing in xGEEV's specification promises that two such sequences round alike, and under LLVM Flang on aarch64 they do not.

The divergence is FMA contraction inside the reference ZGEMM, and it can be isolated to that one file. Recompiling BLAS/SRC/zgemm.f alone, with the rest of the library untouched, clears every failing row:

zgemm.f compiled with	.2d ops	vector FMA	scalar FMA	ZEV test 5
flang 21.1.8 default, -ffp-contract=fast	129	26	20	fails, 1/ulp
flang -ffp-contract=on	140	0	0	passes
flang -ffp-contract=off	140	0	0	passes
flang -fno-vectorize	77	4	13	passes
gfortran 15.2 or 16.1	0	0	30	passes

What is not established is which individual ZGEMM call in the sweep is the first to differ. ZGEMM itself is deterministic. The two paths simply perform different sequences of fused arithmetic. The measured gap is 2 to 14.5 ulp, and every residual test passes, so nothing here suggests a wrong answer - only that the last bits are not reproducible across code paths.

Note finally that any optimised BLAS is free to do the same, and a threaded one certainly will, so the test's premise does not hold for a real deployment either.

The fix

Each of the three comparisons in xDRVEV, and the three in xGET23, becomes

               WNRM = ZERO
               WDIF = ZERO
               DO 150 J = 1, N
                  WNRM = MAX( WNRM, ABS( W( J ) ), ABS( W1( J ) ) )
                  WDIF = MAX( WDIF, ABS( W( J )-W1( J ) ) )
  150          CONTINUE
               RESULT( 5 ) = MAX( RESULT( 5 ), WDIF /
     $                       MAX( UNFL, ULP*REAL( N )*
     $                       MAX( WNRM, WDIF ) ) )

with ABS( WR )+ABS( WI ) in the real drivers, SMLNUM for UNFL in xGET23, and DBLE( N ) in double precision. Including WDIF in the denominator caps the ratio at 1/(n ulp), so a genuinely wrong eigenvalue still reports about 2e14 and the test keeps its full power. The printed description in xDRVEV and xDRVVX changes to match.

Results

Max ZEV ratio over all 1102 cases of zed.in on aarch64:

test	gfortran 15.2 before	after	flang 21.1.8 before	after
1, right-eigenvector residual	5.35	5.35	4.60	4.60
2, left-eigenvector residual	3.89	3.89	4.21	4.21
3, right-eigenvector normalisation	4.49	4.49	5.71	5.71
4, left-eigenvector normalisation	2.00	2.00	2.00	2.00
5, eigenvalue consistency	0.00	0.00	4.5e+15	0.725
6, VR independent of JOBVL	0.00	0.00	0.00	0.00
7, VL independent of JOBVR	0.00	0.00	0.00	0.00

Only test 5 moves, and only under flang, where 4.5e+15 is the 1/ulp sentinel the old code assigned on any difference at all. After the change it reads 0.725. Every other ratio is unchanged to three figures on both compilers, which is the point: nothing but the one over-strict comparison is affected.

Note also that flang's residuals in tests 1 and 2 are better than gfortran's, 4.60 against 5.35 and 4.21 against 3.89 respectively, so the eigenvalues and eigenvectors flang computes are not in question. Only the demand that two code paths agree bit-for-bit was.

Notes

• znep_64 still reports ZHS: 1 out of 2016 under ATfL 22.1. That is xCHKHS test 8, a different check, and it is left alone deliberately: it belongs with SEP test threshold #1350 rather than here.
• Tests 6 and 7 still demand bit equality, of the eigenvector matrices rather than the eigenvalues. They are not currently hit and are left untouched; the same argument would apply to them.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 69.37%. Comparing base (9eaccc1) to head (d1ac8d7).
✅ All tests successful. No failed tests found.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1401   +/-   ##
=======================================
  Coverage   69.36%   69.37%           
=======================================
  Files        6122     6122           
  Lines      486337   486409   +72     
  Branches    23268    23268           
=======================================
+ Hits       337330   337426   +96     
+ Misses     148569   148545   -24     
  Partials      438      438           

Components	Coverage Δ
BLAS	97.94% <ø> (ø)
CBLAS	96.98% <ø> (ø)
LAPACK	82.38% <ø> (ø)
LAPACKE	2.17% <ø> (ø)
TMGLIB	55.69% <ø> (ø)
BLAS testing	88.33% <ø> (ø)
CBLAS testing	89.63% <ø> (ø)
LAPACK testing	82.27% <100.00%> (+0.02%)	⬆️
LAPACKE testing	∅ <ø> (∅)

Files with missing lines	Coverage Δ
TESTING/EIG/cdrvev.f	75.20% <100.00%> (+2.17%)	⬆️
TESTING/EIG/cdrvvx.f	75.53% <ø> (ø)
TESTING/EIG/cget23.f	75.08% <100.00%> (+1.92%)	⬆️
TESTING/EIG/ddrvev.f	76.24% <100.00%> (+2.03%)	⬆️
TESTING/EIG/ddrvvx.f	75.66% <ø> (ø)
TESTING/EIG/dget23.f	75.95% <100.00%> (+1.85%)	⬆️
TESTING/EIG/sdrvev.f	76.24% <100.00%> (+2.03%)	⬆️
TESTING/EIG/sdrvvx.f	75.66% <ø> (ø)
TESTING/EIG/sget23.f	75.60% <100.00%> (+1.86%)	⬆️
TESTING/EIG/zdrvev.f	75.20% <100.00%> (+2.17%)	⬆️
... and 2 more

---

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 9eaccc1...d1ac8d7. Read the comment docs.