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:
zheev_upper_bug.c.txt
