In the LAPACK interface of Accelerate (the default, without ACCELERATE_NEW_LAPACK), zheev with JOBZ='V', UPLO='U' returns eigenvectors that are neither orthonormal nor eigenvectors, while INFO=0. The eigenvalues are correct. zhegv with UPLO='U' fails the same way (tested with B = identity).
It happens when:
- the matrix is complex Hermitian of order N >= 130 (N = 129 is correct, N = 130 already fails, N = 300 fails), and
- its last two or more rows and columns are exactly zero (one zero row is fine). Zero rows at the start or in the middle are fine. Banded and dense matrices both fail.
On the same matrices the following are correct to rounding (1e-15):
- zheev and zhegv with UPLO='L',
- zheevd with UPLO='U',
- dsyev with UPLO='U' on real symmetric matrices of the same shape,
- zheev with UPLO='U' through the new LAPACK interface (-DACCELERATE_NEW_LAPACK).
Since only the eigenvectors are wrong, and only above N=129 (presumably where zhetrd switches to its blocked code), the fault seems to be in the blocked Householder reduction or in the back-transformation (zhetrd/zlatrd or zungtr/zungql) of the legacy library, for columns that are zero.
Asking LAPACK for the optimal LWORK (LWORK = -1 query) does not change the result.
We found this in a physics code (Quanty, quanty.org), where a matrix of this shape arises naturally: a block tridiagonal (block Lanczos) matrix whose last block is padded with zeros after deflation. The wrong eigenvectors silently degraded results by 1e-4 relative, instead of 1e-15.
STEPS TO REPRODUCE
- Save the attached zheev_upper_bug.c.
- clang -O2 zheev_upper_bug.c -o zheev_upper_bug -framework Accelerate
- ./zheev_upper_bug
- For comparison: clang -O2 -DACCELERATE_NEW_LAPACK zheev_upper_bug.c -o zheev_upper_bug_new -framework Accelerate && ./zheev_upper_bug_new
The program builds a random Hermitian matrix (fixed seed, entries in [-0.5,0.5), bandwidth 11 or dense), zeroes its last NZERO rows and columns, calls zheev/zhegv/zheevd and prints max|V^H V - 1| and max|A V - V diag(W)|.
EXPECTED RESULT Both numbers of order 1e-15 for every line, as with the new LAPACK interface.
ACTUAL RESULT (legacy interface) zheev UPLO=U N=129 zero trailing rows=2 INFO=0 max|V^H V - 1| = 4.0e-15 max|A V - V W| = 4.4e-15 ok zheev UPLO=U N=130 zero trailing rows=2 INFO=0 max|V^H V - 1| = 9.8e-01 max|A V - V W| = 1.1e+00 <-- WRONG zheev UPLO=U N=130 zero trailing rows=1 INFO=0 max|V^H V - 1| = 4.7e-15 max|A V - V W| = 3.3e-15 ok zheev UPLO=U N=138 zero trailing rows=5 INFO=0 max|V^H V - 1| = 1.1e+00 max|A V - V W| = 1.2e+00 <-- WRONG zheev UPLO=U N=138 zero trailing rows=5 INFO=0 max|V^H V - 1| = 2.1e-01 max|A V - V W| = 2.2e+00 <-- WRONG zheev UPLO=U N=300 zero trailing rows=2 INFO=0 max|V^H V - 1| = 7.0e-01 max|A V - V W| = 1.5e+00 <-- WRONG zhegv UPLO=U N=138 zero trailing rows=5 INFO=0 max|V^H V - 1| = 1.1e+00 max|A V - V W| = 1.2e+00 <-- WRONG
the same matrices with UPLO = 'L' or zheevd: zheev UPLO=L N=138 zero trailing rows=5 INFO=0 max|V^H V - 1| = 5.9e-15 max|A V - V W| = 3.6e-15 ok zhegv UPLO=L N=138 zero trailing rows=5 INFO=0 max|V^H V - 1| = 4.4e-15 max|A V - V W| = 3.7e-15 ok zheevd UPLO=U N=138 zero trailing rows=5 INFO=0 max|V^H V - 1| = 1.9e-15 max|A V - V W| = 1.4e-15 ok
With -DACCELERATE_NEW_LAPACK every line is "ok" (1e-15).
Code to reproduce the error:
