DGGEV does not converge

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

DGGEV does not converge

Postby zibi14 » Tue Jan 30, 2007 7:26 pm

I noticed that the LAPACK function DGGEV does not converge for this pair of
matrices:

A = [0 0 1; 0 1 0; 1 0 0];
B = [2 0 0; 0 0 1; 0 1 0];

This is a C program I used (linking with CLAPACK ver. 3.0).

int main(void)
{
int i, j;
char jobvl = 'N', jobvr = 'N';
integer nn = 3, one = 1, lwork, info;
doublereal tmp, *work;
doublereal VR, VL, alphar[3], alphai[3], beta[3];
doublereal pA[9] = {0., 0., 1., 0., 1., 0., 1., 0., 0.};
doublereal pB[9] = {2., 0., 0., 0., 0., 1., 0., 1., 0.};

lwork = -1; /* query */
dggev_(&jobvl, &jobvr, &nn, pA, &nn, pB, &nn,
alphar, alphai, beta, &VL, &one, &VR, &one,
&tmp, &lwork, &info);
lwork = (integer) tmp;
printf("\ndggev: lwork = %ld\n", lwork);

work = (doublereal *) malloc(sizeof(doublereal)*lwork);

dggev_(&jobvl, &jobvr, &nn, pA, &nn, pB, &nn,
alphar, alphai, beta, &VL, &one, &VR, &one,
work, &lwork, &info);

printf("dggev: info = %ld\n", info);
return 0;
}

The program returns "info = 3". This matrix pair seems to be simple, and it
is easy to compute eigenvalues of the generalized eigenvalue problem

A.x = s*B.x

I get three eigenvalues (one real and one conjugate pair):
0.793701
-0.39685 - 0.687365*i
-0.39685 + 0.687365*i

Therefore it is interesting to know why DGGEV does not converge.

Zbigniew
zibi14
 
Posts: 9
Joined: Tue Jan 30, 2007 6:57 pm
Location: College Station, Texas

Postby Julien Langou » Wed Jan 31, 2007 11:49 am

Hello Zbigniew,

you are right, LAPACK current DGGEV can fail on surprisingly easy matrix pencils (A,B).
For example, the one you provide:
Code: Select all
     [ 0 0 1 ]        [ 2 0 0 ]
A =  [ 0 1 0 ]    B = [ 0 0 1 ]
     [ 1 0 0 ]        [ 0 1 0 ]

LAPACK DGGEV is unable to find the roots (1/2)^(4/3)* ( -1 +/- sqrt(3).i ) and (1/2)^(1/3).

You can even take an easier one (provided to me by Daniel Kressner who referenced
Vasile Sima)
Code: Select all
     [ 0 0 0 1 ]        [ 1 0 0 0 ]
A =  [ 1 0 0 0 ]    B = [ 0 1 0 0 ]
     [ 0 1 0 0 ]        [ 0 0 1 0 ]
     [ 0 0 1 0 ]        [ 0 0 0 1 ]

Same thing, LAPACK DGGEV is unable to find the roots 1, -1, +i, -i.

This is a known problem. The solution is to apply the so-called 'Exceptional shift
strategy' in the QZ algorithm and LAPACK is currently not doing this.

MATLAB is using a modified version of LAPACK DGEHQZ (the working horse routine of
DGEEV) that implements the exceptional shift strategy. So no suprise that you can get
the correct eigenpairs using MATLAB. Mathworks has kindly provided us with their
implementation but we have encountered side problems when merging both Mathworks
DGEHQZ in LAPACK. So the current LAPACK-3.1.0 still does not have QZ with
Exceptional shift strategy'.

On the midterm, Daniel Kressner and Bo Kagstrom are currently working for a
completely new version of DGEHQZ. This is why few efforts have been spent on the
current routine.

On the short term, a quick patch might do it ...

Julien
Julien Langou
 
Posts: 734
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Postby zibi14 » Wed Jan 31, 2007 1:03 pm

Julien,

Thank you for your prompt reply. This is what I wanted to know - that this is a bug in dggev, and not in my code. And yes, probably a patch will be necessary.

Zbigniew
zibi14
 
Posts: 9
Joined: Tue Jan 30, 2007 6:57 pm
Location: College Station, Texas


Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot], Yahoo [Bot] and 3 guests