LAPACK/ScaLAPACK Development

Hi,

Any ideas why for the following matrix: http://projects.scipy.org/numpy/raw-att ... 54/mat.txt

DGEEV returns vectors for which the relative error |A v - lambda v|/||A||_F is of the order of 1e-2? Test program here: http://projects.scipy.org/numpy/raw-att ... 4/test.f90

And so, with LAPACK 3.2.1,

Code: Select all

$ gfortran -o test test.f90 -llapack -lblas
$ ./test < mat.txt
  4.09736095556244961E-026
  4.43089563288366399E-023
   3.4612399980026800      LARGE ERROR
   8.9078382202313087      LARGE ERROR
   8.8674244567312162      LARGE ERROR
   8.9066982533988046      LARGE ERROR
  3.83826475736065412E-032
...


That's a really cool nasty little matrix that you got there! Congratulations ...
Numpy does not work: http://projects.scipy.org/numpy/ticket/1454
Octave does not work.
Matlab does not work.
LAPACK does not work.
Thanks for reporting this

Code: Select all

clear
A=load('~/Downloads/mat.txt', 'ascii');               
n = 28;
[V,D]=eig(A);                                         
[oo,I] = sort( diag(D) ); V = V(:,I); D = D(I,I);
relgap = zeros(n,1);
relgap(1) = (D(2,2)-D(1,1))/abs(D(1,1));
for i=2:27, relgap(i) = min( [ ((D(i,i)-D(i-1,i-1)))   (D(i+1,i+1)-D(i,i)) ])/abs(D(i,i)); end
relgap(n) = (D(n,n)-D(n-1,n-1))/abs(D(n,n));
fprintf('    eig                 relerr       relgap\n');
for i=1:28, fprintf('%2d %20.16f %e %e\n',i,D(i,i),norm(A*V(:,i)-V(:,i)*D(i,i))/norm(A),relgap(i));, end;


Code: Select all

    eig                 relerr       relgap
 1  -3.1751079149941352 2.033247e-02 4.195976e-16
 2  -3.1751079149941339 2.445700e-02 1.398659e-16
 3  -3.1751079149941335 2.161849e-02 1.398659e-16
 4  -3.1751079149941308 3.274260e-13 8.391952e-16
 5  -2.0498377961923029 1.524756e-02 1.299876e-15
 6  -2.0498377961923002 1.524127e-02 2.166460e-16
 7  -2.0498377961922998 1.468891e-11 2.166460e-16
 8  -2.0498377961922971 1.499134e-02 1.299876e-15
 9  -0.1887301587301591 3.198971e-18 1.176519e-15
10  -0.1887301587301589 1.525037e-18 0.000000e+00
11  -0.1887301587301589 1.664842e-18 0.000000e+00
12  -0.1887301587301586 9.324637e-19 1.617714e-15
13  -0.0495989036825447 1.161271e-03 2.937903e-15
14  -0.0495989036825445 1.165823e-03 1.678802e-15
15  -0.0495989036825444 1.161071e-03 1.678802e-15
16  -0.0495989036825441 1.176836e-03 6.155606e-15
17  -0.0204389846397406 1.378855e-03 3.394931e-15
18  -0.0204389846397405 1.380317e-03 2.376452e-15
19  -0.0204389846397405 1.376691e-03 2.376452e-15
20  -0.0204389846397403 1.377147e-03 7.978087e-15
21  -0.0103199818359591 1.461548e-03 2.151599e-14
22  -0.0103199818359589 1.459206e-03 1.008562e-14
23  -0.0103199818359588 1.466116e-03 1.008562e-14
24  -0.0103199818359586 1.462142e-03 1.512843e-14
25  -0.0077181030285824 1.488029e-03 2.034081e-14
26  -0.0077181030285823 1.480665e-03 9.327554e-15
27  -0.0077181030285822 1.469812e-03 9.327554e-15
28  -0.0077181030285821 1.485513e-03 1.562084e-14


So a lots of cluster of eigenvalues. But, still DGEEV should do its job.
I also looked at "spy(A)" and observed a great structure in the matrix.

Yes, there is a bug in DGEEV. Thanks for reporting it. We will investigate.

Best wishes,
Julien.

In Matlab doing:

Code: Select all

% method 1
[Q,H]=hess(A);                                         
[V,D]=eig(H);
V = Q*V;


instead of

Code: Select all

% method 2
[V,D]=eig(A);


removes the problem.

(If Matlab does what I think it does ...)

Method 2 (eig) calls DGEEV directly. So that means: (1) balancing on A, (2) reduce to Hessenberg form, (3) DGEHRD, etc.
Method 1 does: (1) reduce to Hessenberg form, (2) balancing on H, (3) DGEHRD, etc.

Method 1 works, method 2 does not.

Julien

Following on the previous clue that it might be the balancing the guilty part, it appears that balancing can be turned off in Matlab.

[V,D]=eig(A,'nobalance');

works great.
Julien.

To whom it may concern,
Matlab worked for me. 7.9.0.529 (R2009b)
The matrix is, modulo absolute perturbations the size of the machine precision, a Kronecker product.
A = numpy54;
amax = max( max( abs(A))), % approximately 122
A11 = A(1:7,1:7);
I = eye(4,4);
norm( kron(I,A11) - A )/amax, % approximately eps
[V11,D11] = eig(A11);
V = kron(I,V11);
D = kron(I,D11); % [V,D] = eig(A);
Each eigenvalue has multiplicity 4. The condition number of V
is sensitive to the algorithm used. Balance reduces A11 to
Schur (lower triangular) form.

Julien Langou wrote:That's a really cool nasty little matrix that you got there! Congratulations ... Yes, there is a bug in DGEHRD. Thanks for reporting it. We will investigate.

Very good, thanks! (Though the honor of finding the matrix goes to the person who reported this bug for Numpy :)

This problem is referenced as BUG0057 in our bug list: http://www.netlib.org/lapack/bug_list.html.

The same bug was recently (March 28th, 2013) reported on the forum, see topic 4270: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=4270 and has been reported as BUG0106.

Both bugs (BUG0057 and BUG0106) have the same origin and we fixed the problem in our SVN repository. (See SVN revision 1413 on May 26th, 2013.) We will release the fix in LAPACK 3.5.0.

Cheers,
Julien.