## Non-eigenvectors from DGEEV

Open discussion regarding features, bugs, issues, vendors, etc.

### Non-eigenvectors from DGEEV

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...`
pav

Posts: 6
Joined: Tue Sep 16, 2008 2:58 pm

### Re: Non-eigenvectors from DGEEV

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
`clearA=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)); endrelgap(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-1510  -0.1887301587301589 1.525037e-18 0.000000e+0011  -0.1887301587301589 1.664842e-18 0.000000e+0012  -0.1887301587301586 9.324637e-19 1.617714e-1513  -0.0495989036825447 1.161271e-03 2.937903e-1514  -0.0495989036825445 1.165823e-03 1.678802e-1515  -0.0495989036825444 1.161071e-03 1.678802e-1516  -0.0495989036825441 1.176836e-03 6.155606e-1517  -0.0204389846397406 1.378855e-03 3.394931e-1518  -0.0204389846397405 1.380317e-03 2.376452e-1519  -0.0204389846397405 1.376691e-03 2.376452e-1520  -0.0204389846397403 1.377147e-03 7.978087e-1521  -0.0103199818359591 1.461548e-03 2.151599e-1422  -0.0103199818359589 1.459206e-03 1.008562e-1423  -0.0103199818359588 1.466116e-03 1.008562e-1424  -0.0103199818359586 1.462142e-03 1.512843e-1425  -0.0077181030285824 1.488029e-03 2.034081e-1426  -0.0077181030285823 1.480665e-03 9.327554e-1527  -0.0077181030285822 1.469812e-03 9.327554e-1528  -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.
Last edited by Julien Langou on Thu Apr 22, 2010 1:06 pm, edited 1 time in total.
Julien Langou

Posts: 792
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

### Re: Non-eigenvectors from DGEEV

In Matlab doing:
Code: Select all
`% method 1[Q,H]=hess(A);                                         [V,D]=eig(H);V = Q*V;`

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
Julien Langou

Posts: 792
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

### Re: Non-eigenvectors from DGEEV

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.
Julien Langou

Posts: 792
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

### Re: Non-eigenvectors from DGEEV

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.
dmday

Posts: 2
Joined: Mon Dec 03, 2007 6:18 pm

### Re: Non-eigenvectors from DGEEV

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 :)
pav

Posts: 6
Joined: Tue Sep 16, 2008 2:58 pm

### Re: Non-eigenvectors from DGEEV

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.
Julien Langou

Posts: 792
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA