DGEES gives different results in versions 3.9 and 3.10

I encountered a problem with Lapack library. As a maintainer of scipy library in Fedora I noticed that test suite of scipy fails on Schur decomposition.

I found out that from version 3.10 of Lapack, the DGEES function gives different results in comparison to previous versions. This behavior was caused by re-implementation [1] of DLARTG function.

The DLARTG function still generates a plane rotation so that

    [  C  S  ]  .  [ F ]  =  [ R ]
    [ -S  C  ]     [ G ]     [ 0 ]

    where C**2 + S**2 = 1

but when F is negative, signs on C, S and R are swapped (in comparison to previous implementation). This causes different results of DGEES function.

With my limited mathematical knowledge I cannot judge if the behavior of DGEES is correct, I can only state that its output in versions 3.9 and 3.10 is different (example below).

Was this change of DGEES desired? Does it still produce valid results?

Example (pseudocode): (Full program can be found here, please consider I'm not a Fortran expert, this is my first Fortran program.)

FUNCTION LHS(R,I) result (RES)
    DOUBLE PRECISION R, I
    LOGICAL RES

    RES = R < 0.0
END FUNCTION

A = [  4     3     1   -1   ]
    [ -4.5  -3.5  -1    1   ]
    [  9     6    -4    4.5 ]
    [  6     4    -3    3.5 ]

CALL DGEES('V', 'S', LHS, 4, A, 4, SDIM, WR, WI, VS, 4, WORK, 12, BWORK, INFO)

Lapack 3.9 implementation returns:

A =  [ -.141421356E+01   0.145572146E+00   -.115815859E+02   -.771743518E+01 ]
     [ 0.000000000E+00   -.500000000E+00   0.944719662E+01   -.718440573E+00 ]
     [ 0.000000000E+00   0.000000000E+00   0.141421356E+01   -.145572146E+00 ]
     [ 0.000000000E+00   0.000000000E+00   0.000000000E+00   0.500000000E+00 ]

VS = [ 0.113395338E+00   0.543632173E+00   0.831628257E+00   -.104555566E-13 ]
     [ -.113395338E+00   -.824476349E+00   0.554418838E+00   -.799688216E-14 ]
     [ -.821281688E+00   0.130774409E+00   0.264978233E-01   -.554700196E+00 ]
     [ -.547521126E+00   0.871829392E-01   0.176652155E-01   0.832050294E+00 ]

Lapack 3.10 implementation returns:

A =  [ -.141421356E+01   0.145572146E+00   0.115815859E+02   0.771743518E+01 ]
     [ 0.000000000E+00   -.500000000E+00   -.944719662E+01   0.718440573E+00 ]
     [ 0.000000000E+00   0.000000000E+00   0.141421356E+01   -.145572146E+00 ]
     [ 0.000000000E+00   0.000000000E+00   0.000000000E+00   0.500000000E+00 ]

VS = [ -.113395338E+00   -.543632173E+00   0.831628257E+00   -.104555566E-13 ]
     [ 0.113395338E+00   0.824476349E+00   0.554418838E+00   -.799688216E-14 ]
     [ 0.821281688E+00   -.130774409E+00   0.264978233E-01   -.554700196E+00 ]
     [ 0.547521126E+00   -.871829392E-01   0.176652155E-01   0.832050294E+00 ]

Checklist

[x] I've included a minimal example to reproduce the issue
[ ] I'd be willing to make a PR to solve this issue

Hi @mmuzila,

"Full program can be found here, please consider I'm not a Fortran expert, this is my first Fortran program." => Congratulations for your Fortran program!

With an eye check, the outputs of DGESS 3.9 and 3.10 look good to me, the only difference that I can spot is the signs of the columns of V, so first and second columns have opposite signs (between 3.9 and 3.10), and then third and fourth have the same signs (between 3.9 and 3.10). And then accordingly S(1,3) and S(1,4), and S(2,3) and S(2,4) have opposite signs. This is all good then.

So in short the Schur representation is unique up to a few things: (1) order of the eigenvalues on the diagonal of S, the Schur form (2) signs (or multiplication by a complex number of modulus 1, e^(i.theta), in complex) (3) if we have a repeated eigenvalues, then this is a little more complicated to explain

In the case above, (1) same order on the diagonal, and (3) no repeated eigenvalues, (2) and real case, so this is just up to a possible change of signs.

For the real case, what you want to check after Schur is that, (1) S is upper quasi-triangular with 1-by-1 and 2-by-2 blocks, and the 2-by-2 blocks will be standardized in the form [[ a, b ],[ c a ]] where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc). (2) || A - V * S * V^T || / || A || is small (3) || I - V * V^T || is small.

For the complex case, what you want to check after Schur is that, (1) S is upper triangular (2) || A - V * S * V^H || / || A || is small (3) || I - V * V^H || is small.

For the change of signs, takes any diagonal matrix D with +/- 1 on the diagonal, and then you will note that if the pair V and S is a solution then the pair ( V * D ) and ( D * S * D ) is also a solution since, D^2 = I, and V * S * V^T = ( V * D ) * ( D * S * D ) * ( V * D )^T and V * V^T = ( V * D ) * ( V * D )^T. ( V * D ) means that you can change the signs of any columns of V. ( D * S * D ) means that, e.g, if you changed signs of V1, and did not change signs of V3, then sign of S(1,3) needs to changed.

"I found out that from version 3.10 of Lapack, the DGEES function gives different results in comparison to previous versions. This behavior was caused by re-implementation [1] of DLARTG function." => I am not sure about this. There were major updates in DGEES in itself in 3.10. See: #421. However you are correct that we also changed DLARTG.

@langou, thank you for explanation!

"I found out that from version 3.10 of Lapack, the DGEES function gives different results in comparison to previous versions. This behavior was caused by re-implementation [1] of DLARTG function." => I am not sure about this. There were major updates in DGEES in itself in 3.10. See: #421. However you are correct that we also changed DLARTG.

I tried to use current version of Lapack (3.10) where I changed definition of DLARTG to the old one (3.9). It gave results equal to the old version (3.9), therefore I came to opinion, that the change of DGEES output was caused by DLARTG re-implementation.

I tried to use current version of Lapack (3.10) where I changed definition of DLARTG to the old one (3.9). It gave results equal to the old version (3.9), therefore I came to opinion, that the change of DGEES output was caused by DLARTG re-implementation.

I see. Thanks for explaining. This makes sense. So then, indeed, the changes in the output of DGEES that you observed where due to the DLARTG changes. Good to know. In fairness, and thinking more about it, I am not sure how much the DGEES algorithm has changed from 3.9 to 3.10 for 4x4 matrices. For larger matrices, the DGEES algorithm is very different. I am not sure for 4x4 matrices.

Reference-LAPACK / lapack

DGEES gives different results in versions 3.9 and 3.10 #628