DGGEV does not converge

[zibi14]

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

[Julien Langou]

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

[zibi14]

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