Generalized Eigenvalue Problem

[Eigen_solver]

DEar All,

I find it difficult to run a program, which calls the LAPACK library routines, I know which section of or which group of routines deals with generalized eigenvalue problems in LAPACK but doesn't succeed with setting up the proper compiling with GNU gfortran compiler. Some of the compiler errors are given below for the the eigenvalue example given at the bottom.

Your help will be appreciated!

1

Error: Syntax error in argument list at (1)
main_fortran.f90:39.7:
& 'V', N, A, LDA, B, LDB,
1
Error: Invalid character in name at (1)
main_fortran.f90:40.7:
& ALPHAR, ALPHAI, BETA, SCRATCH, 1, SCRATCH, 1,
1

Code: Select all

  PROGRAM TEST
       IMPLICIT NONE
 !C
       INTEGER           LDA, LDB, LDWORK, M, N
       PARAMETER        (M = 3)
       PARAMETER        (N = 3)
       PARAMETER        (LDA = M)
       PARAMETER        (LDB = M)
       PARAMETER        (LDWORK = 8 * N)
 !C
       DOUBLE PRECISION  A(LDA,N), ALPHAI(N), ALPHAR(N), B(LDB,N)
       DOUBLE PRECISION  BETA(N), SCRATCH, WORK(LDWORK)
       INTEGER           ICOL, INFO, IROW
 !C
       EXTERNAL          DGEGV
       INTRINSIC         ABS, DCMPLX
 !C
 !C     Initialize the arrays A and B to store the matrices A and B
 !C     shown below.
 !C
 !C         1  1  1        1  1  1
 !C     A = 2  1  1    B = 1  0  0
 !C         1  3  1        0  1  0
 !C

        DATA A / 1.0D0, 2.0D0, 1.0D0, 1.0D0, 1.0D0, 3.0D0, 1.0D0, 1.0D0, 1.0D0 /
       DATA B / 1.0D0, 1.0D0, 0.0D0, 1.0D0, 0.0D0, 1.0D0, 1.0D0, 0.0D0, 0.0D0 /
 !C
 !C     Print the initial value of the arrays.
 !C
       PRINT 1000
       PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, N)
       PRINT 1020
       PRINT 1010, ((B(IROW,ICOL), ICOL = 1, N), IROW = 1, N)
 !C
 !C     Compute and print the eigenvalues.
 !C
       CALL DGEGV ('N' ,
      &            'V', N, A, LDA, B, LDB,
      &            ALPHAR, ALPHAI, BETA, SCRATCH, 1, SCRATCH, 1,
      &            WORK, LDWORK, INFO)
       IF (INFO .NE. 0) THEN
         IF (INFO .LT. 0) THEN
           PRINT 1030, ABS(INFO)
           STOP 1
         ELSE
           PRINT 1040, INFO
           STOP 2
         END IF
       END IF
       PRINT 1050
       PRINT 1060, (ALPHAR(IROW), ALPHAI(IROW), BETA(IROW),
      &   IROW = 1, M)
       PRINT 1070
       PRINT 1080, (DCMPLX(ALPHAR(IROW), ALPHAI(IROW)) /
      &             BETA(IROW), IROW = 1, M)
 !C
  1000 FORMAT (1X, 'A:')
  1010 FORMAT (3(3X, F4.1))
  1020 FORMAT (/1X, 'B:')
  1030 FORMAT (1X, 'Illegal argument to DGEGV, argument #', I2)
  1040 FORMAT (1X, 'Internal failure in DGEGV, INFO = ', I4)
  1050 FORMAT (/1X, T6, 'ALPHA', T20, 'BETA')
  1060 FORMAT (1X, '(', F5.2, ',', F5.2, ')', T19, F5.2)
  1070 FORMAT (/1X, 'ALPHA / BETA')
  1080 FORMAT (1X, '(', F5.2, ',', F5.2, ')')
 ! C
       END


[sven]

In Fortran 90, when you want to continue a line, the ampersand needs to be on the line want to continue, not the continuation line (as one did in Fortran 77).

Best wishes,

Sven Hammarling.

[Eigen_solver]

Thank you very much indeed.

1.) I've another question regarding the eigenvalue problems, please advice If need to start a new topic, or I can stick to the current one.
If I've the following (A-lambda*B)*X = C linear algebraic equation system which routine should I use to solve the lambda (eigenvalues) and corresponding X(eigenvectors) in LAPACK library. A and B are non symetric square matrices and C is non zero matrix.

2.) I also have another question regarding the results retrieved from DGGEV and DGEGV subroutines. Below is an outline of results, I feeded the A and B matrix where each is 9x9 and non-symetric, up to the 5th eigenvalue I do get the eigenvalue couples(which normally results in eigenvectors couples) than after 5th eigenvalue It somehow returns to singular eigenvalues and eigennvectors. While I was expecting 9 eigen values and 9 set of eigenvectors for each eigenvalue, Is it normal that DGGEV and DGEGV produces such results or am I missing some setting in the subroutine calls?

Code: Select all

 Eigenvalue( 1) = (-2.1475E-04, 3.8307E-04)

 Eigenvector( 1)
 ( 1.4954E-03,-5.3343E-03)
 ( 2.2766E-04,-5.7528E-04)
 ( 1.2741E-03,-4.0703E-03)
 (-7.4156E-07, 8.7622E-06)
 ( 1.0602E-06,-5.0425E-06)
 ( 8.6351E-05,-2.3437E-04)
 (-8.4175E-01,-1.5825E-01)
 ( 6.5278E-03,-1.6300E-02)
 ( 2.1464E-01,-6.2176E-01)

 Eigenvalue( 2) = (-2.1475E-04,-3.8307E-04)

 Eigenvector( 2)
 ( 1.4954E-03, 5.3343E-03)
 ( 2.2766E-04, 5.7528E-04)
 ( 1.2741E-03, 4.0703E-03)
 (-7.4156E-07,-8.7622E-06)
 ( 1.0602E-06, 5.0425E-06)
 ( 8.6351E-05, 2.3437E-04)
 (-8.4175E-01, 1.5825E-01)
 ( 6.5278E-03, 1.6300E-02)
 ( 2.1464E-01, 6.2176E-01)

.......................
CODE SNIPPET OMITTED HERE
 Eigenvalue( 5) =  1.2051E-07

 Eigenvector( 5)
  2.0277E-01
  6.8011E-03
  2.9270E-01
 -2.0820E-03
  6.8836E-04
  7.3883E-03
 -6.0946E-01
  1.2803E-01
 -1.0000E+00



Your advice will be appreciated.

Regards,

[Eigen_solver]

Does any body can help please !

[sven]

I am afraid that I am not sure how to solve your equation in 1). I have not come across that equation and need to give it some thought, but maybe someone else who knows how to solve the equation will answer.

On your 2), I am very unclear as to what you are asking. What do you mean by 'eigenvalue couples'? Are they not just complex eigenvalues and the others real? By the way, DGEGV is deprecated.

Sven Hammarling.

[CyLith]

Your equation (A-lambda B) X = C is strange. Let's consider a few cases.

1) Suppose (A,B) is a regular matrix pencil, then there exists at least one lambda that makes (A-lambda B) invertible, so X is just inv(A-lambda B)*C. So basically, pick any lambda that is not an eigenvalue of the pair (A,B), and you can always find an X to go with it that solves the equation. This assumes A and B don't have a common nullspace.
2) Suppose (A,B) is a singular matrix pencil, then (A-lambda B) is always singular for any lambda, and there can only be a nontrivial solution X if the columns of C do not lie in the shared nullspace. For instance, if C is not in the nullspace of A, then choosing lambda = 0 and A X = C can be solved.

In general, your equation seems to support an infinite number of solution pairs (lambda, X).

[Eigen_solver]

Dear Sven and CyLith thank you for taking time for helping me,

Explaining a bit further will help you to understand with what I'm dealing and what I'm assuming. First of all I'm not coming from comprehensive mathematical background but from engineering background, sure we studied the maths but not exactly at that level, mostly the PDEs and ODEs. Those type of problems are encountered usually in engineering when performing dynamic analysis, those A, B matrices are called stiffness and mass matrices respectively, the right hand side of equation is loading matrix. Generally simulteneous, subspace etc.. iteration methods are used to solve those kind of problems, I assumed that the eigenvalue subroutines in LAPACK can be used to solve
that. Probably I've overly simplified the problem definition and mis-leaded you.

n your 2), I am very unclear as to what you are asking. What do you mean by 'eigenvalue couples'? Are they not just complex eigenvalues and the others real? By the way, DGEGV is deprecated.

I've used both DGEGV and DGGEV they produced the same result. What I was intending to ask was: Why do I get the two value for the 1st eigenvalue value below, is it normal mathematically? I was expecting -2.1475E-04 or 3.8307E-04 but not both of them.

Eigenvalue( 1) = (-2.1475E-04, 3.8307E-04)

Eigenvector( 1)
( 1.4954E-03,-5.3343E-03)
( 2.2766E-04,-5.7528E-04)
( 1.2741E-03,-4.0703E-03)
.........
(-7.4156E-07, 8.7622E-06)
[/code]

[CyLith]

You have a 9x9 real matrix generalized eigenvalue problem. You expect 9 eigenvalues with corresponding eigenvectors. It is not difficult to show that the eigenvalues will either be real, or come in complex conjugate pairs. Now, DGGEV returns the eigenvalues and eigenvectors in a somewhat unusual way, due to the need to fit all the data into real arrays. For real eigenvalues, the returned ALPHAI is zero, and so the value in ALPHAR is the eigenvalue, with the corresponding column in VL and VR being the left and right eigenvectors. For complex conjugate pairs, ALPHAI is nonzero for two successive values in the array, and the two successive columns in the eigenvector matrices are the real and imaginary parts, where you have to negate the imaginary part (second column) for the second conjugate eigenvalue. Real eigenvalues have real eigenvectors, so they only take up one column. You may already know this, but I'm writing this out just in case, because it wasn't clear from what you said. Basically, the eigenvector matrices do not simply contain real and imaginary parts of each eigenvector, since then they would have to have twice as many columns as they have rows (which they do not; they are square). The storage is compressed so that it all fits in a single square matrix.

That said, since your problem is 9x9, and 9 is an odd number, you must have at least one purely real eigenvalue. Your first eigenvalue is -2.1475E-04 + 3.8307E-04 * i, and your second eigenvalue is -2.1475E-04 - 3.8307E-04 * i, and so on, and you have real eigenvalue 1.2051E-07. If you said that you only have 5 eigenvalues, then you are misinterpreting the output; you should have 9, always. There is however no guarantee about their ordering, so you could have a conjugate pair followed by a singlet real, followed by more conjugate pairs, etc.

[Eigen_solver]

You have a 9x9 real matrix generalized eigenvalue problem. You expect 9 eigenvalues with corresponding eigenvectors. It is not difficult to show that the eigenvalues will either be real, or come in complex conjugate pairs. Now, DGGEV returns the eigenvalues and eigenvectors in a somewhat unusual way, due to the need to fit all the data into real arrays. For real eigenvalues, the returned ALPHAI is zero, and so the value in ALPHAR is the eigenvalue, with the corresponding column in VL and VR being the left and right eigenvectors. For complex conjugate pairs, ALPHAI is nonzero for two successive values in the array, and the two successive columns in the eigenvector matrices are the real and imaginary parts, where you have to negate the imaginary part (second column) for the second conjugate eigenvalue. Real eigenvalues have real eigenvectors, so they only take up one column. You may already know this, but I'm writing this out just in case, because it wasn't clear from what you said.

Thanks for pointing that out. Actually the code omission for sake of clarity created a confusion here, the actual output is below. I do have 9 eigenvalues. If understand your comment correctly my first 4 eigenvalues are complex conjugate pairs(this is also a hint that I didn't accurately populated the inner products of A or B matrix). Since I'm not focusing on complex results, my point of interest is with real eigenvalues which starts from 5th eigenvalue. I think that I do understand the logic of outputs here to some extend.

I've got couple of questions regarding the solution process of the Lapack for eigen problems, I'll be very grateful If you clarify them.
1.) Having the Left and Right eigenvector set probably is correct from mathematical point of view, but engineeringly I'm not supposed to get the left and right vectors, instead I should get simply a set of eigenvectors for each eigenvalue. I got feeling that I'm misinterpreting something here.
2.) Do we have iterative solution methods in Lapack e.g. simultaneous, subspace etc.. ?
3.) It may sound a bit silly question, but in order overcome that C matrix in the (A-lamda*B=C) equation, may I take that in the form of (D=A-C) ; (D-lamda*B)=0

Your help will be appreciated,

Regards,

Code: Select all

Eigenvalue( 1) = (-2.1474E-04, 3.8307E-04)

 Eigenvector( 1)
 ( 1.4954E-03,-5.3343E-03)
 ( 2.2765E-04,-5.7528E-04)
 ( 1.2741E-03,-4.0702E-03)
 (-7.4154E-07, 8.7621E-06)
 ( 1.0602E-06,-5.0425E-06)
 ( 8.6351E-05,-2.3437E-04)
 (-8.4176E-01,-1.5824E-01)
 ( 6.5278E-03,-1.6300E-02)
 ( 2.1464E-01,-6.2176E-01)

 Eigenvalue( 2) = (-2.1474E-04,-3.8307E-04)

 Eigenvector( 2)
 ( 1.4954E-03, 5.3343E-03)
 ( 2.2765E-04, 5.7528E-04)
 ( 1.2741E-03, 4.0702E-03)
 (-7.4154E-07,-8.7621E-06)
 ( 1.0602E-06, 5.0425E-06)
 ( 8.6351E-05, 2.3437E-04)
 (-8.4176E-01, 1.5824E-01)
 ( 6.5278E-03, 1.6300E-02)
 ( 2.1464E-01, 6.2176E-01)

 Eigenvalue( 3) = (-1.8696E-07, 1.3561E-07)

 Eigenvector( 3)
 ( 2.9570E-01, 1.0354E-01)
 (-3.6520E-03, 1.5151E-03)
 ( 4.5721E-02, 2.0330E-02)
 (-5.2474E-04,-1.0180E-04)
 ( 1.9830E-04, 5.2686E-05)
 (-1.0456E-04, 8.0890E-04)
 (-3.8798E-02,-1.2676E-01)
 (-1.1402E-01, 4.2922E-02)
 ( 6.9592E-01,-3.0408E-01)

 Eigenvalue( 4) = (-1.8696E-07,-1.3561E-07)

 Eigenvector( 4)
 ( 2.9570E-01,-1.0354E-01)
 (-3.6520E-03,-1.5151E-03)
 ( 4.5721E-02,-2.0330E-02)
 (-5.2474E-04, 1.0180E-04)
 ( 1.9830E-04,-5.2686E-05)
 (-1.0456E-04,-8.0890E-04)
 (-3.8798E-02, 1.2676E-01)
 (-1.1402E-01,-4.2922E-02)
 ( 6.9592E-01, 3.0408E-01)

 Eigenvalue( 5) =  1.2051E-07

 Eigenvector( 5)
  2.0277E-01
  6.8011E-03
  2.9270E-01
 -2.0820E-03
  6.8836E-04
  7.3883E-03
 -6.0946E-01
  1.2803E-01
 -1.0000E+00

 Eigenvalue( 6) =  3.0218E-10

 Eigenvector( 6)
  8.0521E-01
  6.1644E-02
  1.7393E-01
 -8.2557E-04
  4.7075E-04
  4.7113E-03
 -8.4556E-01
  1.2357E-01
 -1.0000E+00

 Eigenvalue( 7) = -1.1468E-14

 Eigenvector( 7)
  1.6418E-01
  4.0329E-02
  6.6482E-02
  1.0555E-02
 -4.1019E-02
  1.6015E-03
  6.1565E-02
  1.3209E-01
 -1.0000E+00

 Eigenvalue( 8) =  2.4515E-16

 Eigenvector( 8)
  3.4255E-01
 -4.3142E-01
  2.0170E-01
 -1.6028E-01
  9.4246E-02
  8.5530E-01
 -5.2758E-04
 -4.7140E-01
 -1.0000E+00

 Eigenvalue( 9) = -3.8840E-16

 Eigenvector( 9)
 -2.9768E-01
 -1.5401E-02
 -3.9154E-02
  4.9394E-01
  7.4078E-02
 -1.0000E+00
  1.7489E-01
  5.9632E-01
 -9.3858E-01


[CyLith]

1. What people mean by THE eigenvectors is usually the right eigenvectors. The left eigenvectors are THE eigenvectors of the transposed problem.
2. All eigenvalue solution methods are iterative. But I think what you mean is does Lapack contain Krylov subspace type methods like Arnoldi, Lanczos, or Jacobi-Davidson where you only need to apply matrix-vector products. The answer is no. Lapack contains algorithms for dense problems which requires you to have access to the full matrix. For DGGEV, it first performs a QR factorization of B, and applies the Q to A, then it gets the new A into upper Hessenberg form. At this point, it does QZ iteration to get both matrices into upper triangular form (generalized Schur form). The eigenvalues can then be read off the diagonals, and you can perform backward substitution to get the eigenvectors.
3. So... If your problem comes from discretizing a mechanical system, where B is a mass matrix and A is a stiffness matrix, then usually the eigenvalues have the interpretation of the (squared) resonant frequencies of the system. For all of this the right hand side C = 0. If you want to then solve for the response of the system subject to external forcing, then C is nonzero, and you have to choose a (squared) frequency lambda at which to compute the response. I think it would help immensely if you gave the actual background information for the problem you're trying to solve. And perhaps this discussion is best taken over to scicomp.stackexchange.com, where actual math notation can be used.

[Eigen_solver]

Thanks CyLith for the valuable information,

2.)I recently learned that general eigenvalue problems are iterative by nature due to Abbel's proff. Then can I use the Arpack library for the subspace iterations, AFAIU Arpack uses the Blas and Lapack core files extensively and in the web site of Arpack it shows that it's capable to implement the iterative Kyrlov subsapce, Arnoldi methods.
Is that forum is also valid for Arpack questions?

3.) Problem comes from solid mechanics and structural dynamics. The origin is from motion equation M*X''(t)+C*X'(t)+K*X(t)=F(t), left side is mass, damping and stiffness matrices respectively. In my case I don't have a damping matrix for sake of simplicity (I'll take it to next step If can eliminate the current) and the solution I'm trying to obtain is not time dependent (also this is the next step), that is my X(t) is constant. So this morphs into eigenvalue problem somehow ( Galerkin approach, Duhammel integral soution, Green theorm etc..) and I end up into (A-lamda*B=C). So I'm not sure that left side can be eliminated by simply subtracting from LHS. I've tried to summarize it, if its still not verbose and elaborated enough I'll try to explain it as you suggest at stackexchange.


Can you shed some light also on Newton-Raphson method availability in Lapack for linear algebraic equation systems?

Regards,

[CyLith]

2) Yes, the Abel-Ruffini theorem says you can't generally solve eigenvalue problems of dimension greater than 4 algebraically. Arpack is for solving large sparse eigenvalue problems. Generally this site is not for asking about Arpack (but really I don't think there is a great place to ask, so you might as well ask here). Arpack basically projects the huge problem onto smaller subspaces and uses Lapack to solve the small problems. It is not very easy to use, and if you want to find the smallest eigenvalues, then it's really only good and finding a handful. It gets very slow the more eigenpairs you want.

3) Let me preface this with the fact that I'm an electrical engineer, so I may be completely wrong. Your equation is not an eigenvalue problem. If you are seeking steady state solutions, you assume X(t) = x(w) e^{i w t}, where w (omega) is the angular frequency (essentially Fourier transforming the equation). You then get an equation -M*Z*w^2 + K*Z = f(w), where I'm using lowercase letters to represent the Fourier transformed quantities. Now in the absence of forcing (F = 0), you have an eigenvalue problem in frequency squared (you are trying to find the resonant modes of the system). With forcing, this is not an eigenvalue problem, it is just a linear system, and you must find the spectral response of the system for all frequencies w. This is not equivalent to [ (A-C) - lambda*B ] x = 0, since C (the forcing function) should not multiply x. So, if you really are after the spectral response, then you must solve the linear system A(w)x=f for all the frequencies w you are interested in, where A(w) = K-M*w^2. Then the solution x(w) at each frequency is the steady state excitation of whatever degrees of freedom it represents.

Lapack is focused on dense linear algebra, and has no Newton-Raphson method. You could write a Newton method using Lapack to perform the basic operations.