search

modelica / ModelicaStandardLibrary

Free (standard conforming) library to model mechanical (1D/3D), electrical (analog, digital, machines), magnetic, thermal, fluid, control systems and hierarchical state machines. Also numerical functions and functions for strings, files and streams are included.
https://doc.modelica.org
BSD 3-Clause "New" or "Revised" License
479 stars 169 forks source link

Bug in calculating eigenvectors with "eigenValues/dgeev" #918

Closed modelica-trac-importer closed 7 years ago

modelica-trac-importer commented 7 years ago

Reported by bastian on 4 Dec 2012 09:03 UTC dgeev from LAPACK is called as

    external "Fortran 77" dgeev("N", "V", n, Awork, n, eigenReal, eigenImag,
        eigenVectors, n, eigenVectors, n, work, size(work, 1), info)
        annotation (Library="Lapack");

i.e. eigenVectors is used twice (for the left as well as the right eigenvectors). By chance the first occurence (left eigenvectors) is used -- but these are not computed ("N"). So eigenVectors has the wrong entries (uninitialised or 0).

Additionally I found some inconsistencies in the documentation of some LAPACK functions.

The documentation for dgegv says Compute generalized eigenvalues and eigenvectors for a (A,B) system But this routine is called so that no eigenvectors are calculated.

However, dgegv is deprecated and dggev should be used. But here the documentation says Compute generalized eigenvalues for a (A,B) system Instead the left and right eigenvalues are calculated additionally. Moreover the purpose of the input parameter nA for dggev is not clear to me. I assume it is not needed and can be replaced by n in the call to the external dggev routine.

I attach a diff file with my changes.

Migrated-From: https://trac.modelica.org/Modelica/ticket/918

modelica-trac-importer commented 7 years ago

Comment by hansolsson on 5 Dec 2012 17:12 UTC No obvious error, but clearly the documentation is lacking:

The dgegv seems to work, since if VL="N" the left eigenvectors should not be referenced and thus not set to zero. But it seems odd to not have it documented.

The dggev could have nA to compute eigenvalues for a smaller sub-matrix; but the intent of the input is not documented.

modelica-trac-importer commented 7 years ago

Comment by otter on 8 Dec 2012 23:31 UTC Fixed in 6829d57388a2167ccf5a8b9eb56b597c7421faf3

Replying to [ticket:918 bastian@…]:

dgeev from LAPACK is called as

external "Fortran 77" dgeev("N", "V", n, Awork, n, eigenReal, eigenImag,
eigenVectors, n, eigenVectors, n, work, size(work, 1), info)
annotation (Library="Lapack");

i.e. eigenVectors is used twice (for the left as well as the right eigenvectors). By chance the first occurence (left eigenvectors) is used -- but these are not computed ("N"). So eigenVectors has the wrong entries (uninitialised or 0).

The LAPACK documentation states that the memory for the left eigenvectors (VL) is not referenced when option JOBVL="N". So, in principal the call is correct (since the provided memory for right eigenvectors is not utilized at the left eigenvector argument). However, I agree that this is not a good programming style and followed your advice to explicitly declare a left eigenvector array.

I also added to the documentation that this is not a complete interface to the LAPACK function.

Additionally I found some inconsistencies in the documentation of some LAPACK functions.

The documentation for dgegv says Compute generalized eigenvalues and eigenvectors for a (A,B) system But this routine is called so that no eigenvectors are calculated.

Corrected the documentation.

However, dgegv is deprecated and dggev should be used. But here the documentation says Compute generalized eigenvalues for a (A,B) system Instead the left and right eigenvalues are calculated additionally.

Corrected the documentation.

Moreover the purpose of the input parameter nA for dggev is not clear to me. I assume it is not needed and can be replaced by n in the call to the external dggev routine.

The actual computation is performed for A[1:nA,1:nA]. I checked and it seems that this function is used to compute the transmission zeros in Modelica_LinearSystems2. However, where it is called, nA is called with size(A,1). So, I also would prefer to remove this argument. However, for backwards compatibility this is not possible and it is also not possible to just ignore it, if there is still a place where this feature is used.

I just documented the behavior

I attach a diff file with my changes.

modelica-trac-importer commented 7 years ago

Comment by otter on 9 Dec 2012 14:48 UTC Additional fix in 2abb99ea8fd048c9a8588190638b087d72979523: Instead of introducing a dummy array for the left eigenvectors of dimension N, a dummy array of dimension 1 is introduced and reported to the LAPACK function (so, less memory consumption, and still clean implementation, since a leading dimension of 1 is reported).