DSYEVR returns non-orthogonal vectors?

[pav]

Hi,

(Forwarding http://projects.scipy.org/scipy/ticket/1159 )

For some matrices DSYEVR, when told to compute a small number of eigenvalues, returns non-orthogonal eigenvectors. Perhaps it possible for the algorithm to detect this and return with an error?

This applies at least when using the automatically chosen ABSTOL. However, at least for the example below, I wasn't able to find a suitable ABSTOL for which it would work.

Maybe, this might have to do with the fact that for the example problem below, the first eigenvalue is multiply degenerate, and the numerical second eigenvalue is eps-close to the first one.

The following test program illustrates this: http://projects.scipy.org/scipy/raw-att ... 9/test.f90
with the test data http://projects.scipy.org/scipy/raw-att ... 1159/L.txt

On Lapack 3.2.1 + reference BLAS, I get:

Code: Select all

$ gfortran -o test test.f90 -llapack -lblas
$ ./test < L.txt
 orth error (           2  eigs): -0.96599125529927221      HUGE ERROR
 orth error (           3  eigs): -0.96599125529927232      HUGE ERROR
 orth error (           4  eigs): -0.96599125529927221      HUGE ERROR
 orth error (           5  eigs): -2.52928802610462387E-017
 orth error (           6  eigs): -1.47893223140021721E-016
 orth error (           7  eigs): -1.18177613284446356E-016
 orth error (           8  eigs):   0.0000000000000000     
 orth error (           9  eigs): -2.84549495443531079E-016
 orth error (          10  eigs): -2.93456682158620488E-016


[admin]

Dear Pauli,
I just ran your test and it works fine on my machine.
I tried LAPACK 3.2.1, 3.1.1 and our current repository version with different BLAS (ref BLAS, veclib and ATLAS).
gfortran -o test test.f90 ~/Documents/Boulot/lapack-dev/lapack/tags/lapack-3.2.1/lapack_LINUX.a ~/Documents/Boulot/lapack-dev/lapack/tags/lapack-3.2.1/blas_LINUX.a
./test < L.txt
orth error ( 2 eigs): -2.67037960893108216E-032
orth error ( 3 eigs): 4.46452343408814271E-033
orth error ( 4 eigs): -2.67037960893108216E-032
orth error ( 5 eigs): 2.73435787900844218E-016
orth error ( 6 eigs): 1.45283091113057594E-017
orth error ( 7 eigs): -7.06450011644714449E-033
orth error ( 8 eigs): 0.0000000000000000
orth error ( 9 eigs): 3.53775168882020097E-016
orth error ( 10 eigs): -9.49761103097301884E-017

So I suspect a compilation issue, or something like this.
Could you tell me which configuration you are using? (machine, compiler, make.inc, etc..)

Regards,
Julie

[pav]

The issue occurred with Lapack+BLAS libraries shipped by Debian (libatlas3gf-base 3.6.0-22 + liblapack3gf 3.1.1-1), and also on those by Ubuntu (libblas3gf 1.2-2build1 + 3.2.1-2). Both are gfortran-built.

I also downloaded LAPACK 3.2.1 (http://netlib.org/lapack/lapack.tgz) and BLAS (http://netlib.org/blas/blas.tgz) built it using the default make.inc's as-is (except changing BLAS's to use gfortran). This means "-g" for LAPACK and "-O3" for BLAS.

On the Ubuntu machine (32-bit linux intel atom), I have

Code: Select all

$ gfortran --version
GNU Fortran (Ubuntu 4.4.3-4ubuntu5) 4.4.3
$ gfortran -o test test.f90 lapack-3.2.1/lapack_LINUX.a lapack-3.2.1/blas_LINUX.a
$ ./test < L.txt
 orth error (           2  eigs):  0.98018274937637251      HUGE ERROR
 orth error (           3  eigs):  0.98018274937637251      HUGE ERROR
 orth error (           4  eigs):  0.98018274937637251      HUGE ERROR
 orth error (           5  eigs):  1.25314924656247904E-016
 orth error (           6  eigs):  2.17735748322348682E-016


The same also with BLAS compiled with "-g" instead of "-O3". Moreover, the same occurs also with everything built with

Code: Select all

$ gfortran-4.3 --version
GNU Fortran (Ubuntu 4.3.4-10ubuntu1) 4.3.4


On the Debian machine (64-bit intel xeon), I get from the same procedure

Code: Select all

$ gfortran --version
GNU Fortran (Debian 4.3.2-1.1) 4.3.2
$ gfortran -o test test.f90 lapack-3.2.1/lapack_LINUX.a lapack-3.2.1/blas_LINUX.a
$ ./test < L.txt
 orth error (           2  eigs): -2.67037960893108216E-032
 orth error (           3  eigs):  4.46452343408814271E-033
 orth error (           4  eigs): -2.67037960893108216E-032
 orth error (           5  eigs):  2.73435787900844218E-016
 orth error (           6  eigs):  1.45283091113057594E-017
 orth error (           7  eigs): -7.06450011644714449E-033


So no error there, and adding "-O2" or "-O3" to LAPACK compilation doesn't seem to break anything.
Moreover,

Code: Select all

$ gfortran -o test test.f90 lapack-3.2.1/lapack_LINUX.a -lblas-3
22:59 pauli@phys:~/tmp$ ./test < L.txt
 orth error (           2  eigs): -1.80363506134932581E-032
 orth error (           3  eigs):  0.99999999999998068      HUGE ERROR
 orth error (           4  eigs):  0.99999999999998068      HUGE ERROR
 orth error (           5  eigs):   0.0000000000000000     
 orth error (           6  eigs):   0.0000000000000000     
 orth error (           7  eigs):   0.0000000000000000     
 orth error (           8  eigs): -3.13768108717304983E-016
 orth error (           9  eigs): -4.51028103753969845E-017
 orth error (          10  eigs):  1.40946282423115576E-016
$ gfortran -o test test.f90 lapack-3.2.1/lapack_LINUX.a BLAS/blas_LINUX.a
$ ./test < L.txt
 orth error (           2  eigs): -2.67037960893108216E-032
 orth error (           3  eigs):  4.46452343408814271E-033
 orth error (           4  eigs): -2.67037960893108216E-032
 orth error (           5  eigs):  2.73435787900844218E-016
$ gfortran -o test test.f90 -llapack-3 BLAS/blas_LINUX.a
$ ./test < L.txt
 orth error (           2  eigs):  9.82816896344506035E-019
 orth error (           3  eigs):  9.82816896344506035E-019
 orth error (           4  eigs):  9.82816896344506035E-019
 orth error (           5  eigs):  1.24132789084528224E-017


So it seems like the Debian-supplied BLAS is somehow faulty, at least in this case.

A third machine (Debian Lenny, Athlon 64-bit), shows similar behavior as the above 64-bit machine.

So, there seem to be problems with 32-bit self-compiled and all distributor-supplied packages. So it looks very much like a compiler-related problem, unless there's some bitness issues...

[admin]

Pav,
Yes I confirm, I managed to reproduce the problem on my 32 bits machine with 4.1.2 and the reference BLAS
I am going to take a closer look.
Julie

[miro_ilias]

Greetings, dear Lapack experts,

I am co-developer of the scientific code DIRAC, which can be improved by linking against proper DSYEVR routine.

We have discussed DSYEVR's non-orthogonal vectors at our web-manual,
http://diracprogram.org/doc/master/programmers/diag.html#independent-demonstration-program .

Following the post
http://www.netlib.org/lapack/Errata/index2.html#_strong_span_class_green_bug0056_span_strong_non_orthogonal_eigenvectors_returned_from_dsyevr

I (we) though that the non-orthogonal vectors isse got fixed (by Osni on Mon May 4th 2015) in the development version. So I cloned the fresh lapack (last commit from 17.6.2015) from https://github.com/live-clones/lapack and compiled it with GNU compilers.

However, with our demonstration program (http://diracprogram.org/doc/master/programmers/diag.html#independent-demonstration-program), I got non-orthogonal eigenvectors (gfortran) again:

Code: Select all

 **** LAPACK DSYEVR ****
 U^{+}*A*U - eps ?= 0> norm/diag:0.1616D-14  norm/offdiag:0.1419D-06
     U^{+}*U - I ?= 0> norm/diag:0.3454D-15  norm/offdiag:0.2838D-06
     U*U^{+} - I ?= 0> norm/diag:0.1047D-05  norm/offdiag:0.7950D-06

Any help, please ?

[osni]

This is an update about this potential bug:

1) If I run the test code provided by the user, in which DSYEVR is called with UPLO = 'U', I also get vectors that fail the orthogonality test.
2) If I run the test code provided by the user but changing it to call DSYEVR with UPLO = 'L', everything works fine (i.e. the vectors are orthogonal)
3) If I take the tridiagonal generated through 1) and run a separate tester everything works fine (i.e. the vectors are orthogonal).

There is something weird going on but I need to dig a little more before declaring it a bug. By the way, if there is a bug, I don't think it is related to bug 56: in bug 56, dsyevr is called with RANGE = 'I' (the IL-th through IU-th eigenvalues will be found), which triggers a different set of calculations, while in the test code provide by the user dsyevr is called with RANGE = 'A'.

Osni

[admin]

Dear Miro

Osni managed to reproduce your issue.

Here is his update

I managed to reproduce the bug with our standalone tester, after saving the tridiagonal matrix generated by DSYTRD in DSYEVR with a lot of digits... Looking at the eigenvectors computed by DSTEMR (which is called by DSYEVR), it seems that the problem is related to some splitting of the tridiagonal. (I compared the eigenvectors computed by DSTEMR and DSTEVX to check that hypothesis. Note, however, that this bug is not related to bug 0056.)

We opened a new bug for it - bug 126 - see http://www.netlib.org/lapack/Errata/index2.html

Thank you again for your hard work
We are going to further investigate the issue and hopefully correct it very soon
We will keep you updated
Julie

[miro_ilias]

Hi,

I put the dsyevr test here,

https://github.com/miroi/lapack-dsyevr-test

Lapack 3.6.0 still contains this bug, labeled as bug 126.

[TriMike]

Hi,

I put the dsyevr test here,

https://github.com/miroi/lapack-dsyevr-test

Lapack 3.6.0 still contains this bug, labeled as bug 126.

Hi Miro, is this bug still in Lapack or was it fixed?

[jakers87]

Hi,

I put the dsyevr test here,

https://github.com/miroi/lapack-dsyevr-test

Lapack 3.6.0 still contains this bug, labeled as bug 126.

Is there any fix for that bug yet?

[miro_ilias]

Hi,
Is there any fix for that bug yet?

No, https://github.com/miroi/lapack-dsyevr-test , this bug remains with lapack-3.7.0, see also https://github.com/Reference-LAPACK/lapack/issues/151