Complex EigenVectors

[mjerdem]

Hi,

I am new to using fortran. In my work I need to calculate the eigenvalues and eigenvectors of various matrices. Sometimes they have complex eigenvectors sometimes real. You can find my simple code that calculates a simple example of 2*2 matrix's eigenvectors and eigenvalues. They are supposed to be complex, however, I do not get complex eigenvectors. Can you please let me know what I am doing wrong.

Thank you

(please do not worry about the extra inputs. This is actually a part of a bigger code)

program sec
implicit none


real(8),dimension,,allocatable :,EigVecint
real(8),dimension),allocatable ::WORK, EigValintr,EigValinti
integer :: INFO, Nz,i,j
CHARACTER (LEN=100):: D2intmat, InvD2intmat, EigValintmat, D2eigintmat,cornerinvmat, cornermat
CHARACTER (LEN=100) :: D2intcolmat, D2introwmat, EigVecintmat

allocate (D(2,2), EigValintr(2), EigValinti(2), EigVecint(2,2), WORK(4*(2)))
Nz=3

D(1,1)=3
D(2,1)=4
D(1,2)=-2
D(2,2)=-1


cornermat='cornermat'
OPEN(16,FILE=cornermat)
DO i=1,2
WRITE(16,*) (D(i,j),j=1,2)
ENDDO
CLOSE(16)


CALL DGEEV('N','V',2,D,2,EigValintr,EigValinti,EigVecint,2,EigVecint,Nz-1,WORK,4*(2),INFO)

EigValintmat='EigValintmat'
OPEN(21,FILE=EigValintmat)
DO i=1,Nz-1
WRITE(21,*)i,EigValintr(i),EigValinti(i)
ENDDO
CLOSE(21)

EigVecintmat='EigVecintmat'
OPEN(22,FILE=EigVecintmat)
DO i=1,Nz-1
WRITE(22,*) (EigVecint(i,j),j=1,Nz-1)
ENDDO
CLOSE(22)


end program sec

 

[salgerman]

without looking too much into the program...I don't see any variable in your program having been declared COMPLEX? How are you supposed to get complex number if you don't have one? I presume that when your variable gets assigned a complex number it gets typecast down to real and the imaginary part gets lost.

 

[gummibaer]

mjerdem,

this does not seem to be a problem of your code, for it does just read and write input and output respectively.
The question is your call to DGEEV and what all those parameters indicate what DGEEV should do with the data.

Best is, you compare these parameters to what your documentation on DGEEV says and find out if everything is as it should be.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[gummibaer]

p. s.

I just looked up DGEEV in lapack documentation. To distinguish between left and right eigenvectors is new to me. What happens if you exchange 'V' and 'N' in the parameters? What happens if you set both to 'V'?

What does INFO show after DGEEV executed?

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[mjerdem]

Thank you gummibaer and salgerman

I forget to asign Eigenvector as complex for sure; however, that did not solved my problem. I made that change and corrected the input matrix.

I am interested in right eigenvector in my work that is why I am not calculating the left eigenvector.

The right eigenvector v(j) of A satisfies the following formula:

A*v(j) = lambda(j)*v(j)

where

lambda(j) is its eigenvalue


(left eigenvector u(j)H*A = lambda(j)*u(j)H where u(j)H denotes the conjugate transpose of u(j). )

This is a really simple 2by2 matrix I do not know why I can not get the imaginary eigenvector while getting the right eigenvalues.

I really appreciate your help.


program sec
implicit none

real(8),dimension,,allocatable ::A
real(8),dimension),allocatable ::WORK, EigValinti, EigValintr
complex (8),dimension,,allocatable :: EigVecint
integer :: INFO, Nz,i,j
CHARACTER (LEN=100):: D2intmat, InvD2intmat, EigValintmat, D2eigintmat,cornerinvmat, cornermat
CHARACTER (LEN=100) :: D2intcolmat, D2introwmat, EigVecintmat

Nz=2
allocate (A(Nz,Nz), EigValintr(Nz), EigValinti(Nz), EigVecint(Nz,Nz), WORK(4*(Nz)))

A(1,1)=3
A(2,1)=4
A(1,2)=-2
A(2,2)=-1

cornermat='cornermat'
OPEN(16,FILE=cornermat)
DO i=1,Nz
WRITE(16,*) (A(i,j),j=1,Nz)
ENDDO
CLOSE(16)

CALL DGEEV('N','V',Nz,A,Nz,EigValintr,EigValinti,EigVecint,Nz,EigVecint,Nz,WORK,4*(Nz),INFO)

EigValintmat='EigValintmat'
OPEN(21,FILE=EigValintmat)
DO i=1,Nz
WRITE(21,*)i,EigValintr(i),EigValinti(i)
ENDDO
CLOSE(21)

EigVecintmat='EigVecintmat'
OPEN(22,FILE=EigVecintmat)
DO i=1,Nz
WRITE(22,*) (EigVecint(i,j),j=1,Nz)
ENDDO
CLOSE(22)

end program sec
INCLUDE 'libhar.f95'
INCLUDE 'sgedi.f95'

 

[gummibaer]

I guess you should look up the docs, e.g. [link ]

http://www.physics.orst.edu/~rubin/nacphy/lapack/routines/dgeev.html

[/url]

I take it that:

Real and imaginary parts of the eigenvalues are saved as double precision arrays (=real) entities.

Then there is this part on the left eigenvector:

VL	  (output) DOUBLE PRECISION array, dimension (LDVL,N)
	  If JOBVL = 'V', the left eigenvectors	u(j) are stored	one after
	  another in the columns of VL,	in the same order as their eigen-
	  values.  If JOBVL = 'N', VL is not referenced.  If the j-th
	  eigenvalue is	real, then u(j)	= VL(:,j), the j-th column of VL.  If
	  the j-th and (j+1)-st	eigenvalues form a complex conjugate pair,
	  then u(j) = VL(:,j) +	i*VL(:,j+1) and
	  u(j+1) = VL(:,j) = i*VL(:,j+1).

I understand this, that of conjugate complex eigenvectors you have both real and imaginary part one after the other in the array holding your eigenvectors,

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[mikrom]

I found this example

http://software.intel.com/sites/pro.../mkl_sa/11/mkl_lapack_examples/dgeev_ex.f.htm

and adapted it for your matrix

http://www.sosmath.com/matrix/eigen3/eigen3.html

I got these results:

 Eigenvalues: 
 (  1.00,  2.00) (  1.00, -2.00)

 Left eigenvectors:
 (  0.82,  0.00) (  0.82, -0.00)
 ( -0.41, -0.41) ( -0.41,  0.41)

 Right eigenvectors:
 (  0.41,  0.41) (  0.41, -0.41)
 (  0.82,  0.00) (  0.82, -0.00)

 

[mikrom]

Btw, I used this code

eigenvectors.f95

[COLOR=#a020f0]Program[/color] EigenVectors
  [COLOR=#0000ff]! Finding Eigenvalues and Eigenvectors of a matrix A using LAPACK[/color]
  [COLOR=#2e8b57][b]Implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]

  [COLOR=#0000ff]! declarations[/color]
  [COLOR=#2e8b57][b]integer[/b][/color] :: n, lda, ldvl, ldvr, lwork, lwmax, info
  [COLOR=#2e8b57][b]parameter[/b][/color](n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color], lwmax [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1000[/color])
  [COLOR=#2e8b57][b]parameter[/b][/color](lda [COLOR=#804040][b]=[/b][/color] n, ldvl [COLOR=#804040][b]=[/b][/color] n, ldvr[COLOR=#804040][b]=[/b][/color] n)
  [COLOR=#2e8b57][b]double precision[/b][/color] :: A(lda,n), VL(ldvl,n), VR(ldvl,n), [highlight #ffff00][COLOR=#0000ff]&[/color][/highlight]
                      WR(n), WI(n), work(lwmax)
  
  [COLOR=#0000ff]! matrix A[/color]
 
  A([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]3.00[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2.00[/color][COLOR=#804040][b]/[/b][/color])
  A([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]4.00[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1.00[/color][COLOR=#804040][b]/[/b][/color])  

  [COLOR=#0000ff]! workspace query: calculates the optimal size of the WORK array[/color]
  lwork [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
  [COLOR=#a020f0]call[/color] DGEEV([COLOR=#ff00ff]'Vectors'[/color], [COLOR=#ff00ff]'Vectors'[/color], n, A, lda, WR, WI, VL, LDVL, [highlight #ffff00][COLOR=#0000ff]&[/color][/highlight]
              VR, LDVR, work, lwork, info)
  [COLOR=#0000ff]! compute workspace size[/color]
  lwork [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]min[/color](lwmax, [COLOR=#008080]int[/color](work([COLOR=#ff00ff]1[/color])))
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'lwork = '[/color], lwork

  [COLOR=#0000ff]! solve[/color]
  [COLOR=#a020f0]call[/color] DGEEV([COLOR=#ff00ff]'Vectors'[/color], [COLOR=#ff00ff]'Vectors'[/color], n, A, lda, WR, WI, VL, LDVL, [highlight #ffff00][COLOR=#0000ff]&[/color][/highlight]
             VR, LDVR, work, lwork, info)

  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'info = '[/color], info
  [COLOR=#804040][b]if[/b][/color] (info [COLOR=#804040][b].eq.[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Solution was succesfully computed:'[/color]
    [COLOR=#0000ff]! print the solution[/color]
    [COLOR=#a020f0]call[/color] print_eigenvalues([COLOR=#ff00ff]'Eigenvalues: '[/color], n, WR, WI )
    [COLOR=#a020f0]call[/color] print_eigenvectors([COLOR=#ff00ff]'Left eigenvectors:'[/color], n, WI, VL, ldvl )
    [COLOR=#a020f0]call[/color] print_eigenvectors([COLOR=#ff00ff]'Right eigenvectors:'[/color], n, WI, VR, ldvr )
  [COLOR=#804040][b]else[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'* Error computing solution!'[/color]
  [COLOR=#804040][b]end if[/b][/color]
[COLOR=#a020f0]end program[/color] EigenVectors

[COLOR=#0000ff]!-------------------------- Subroutines --------------------------[/color]
[COLOR=#a020f0]subroutine[/color] print_eigenvalues (desc, n, WR, WI)
  [COLOR=#2e8b57][b]character[/b][/color][COLOR=#804040][b]*[/b][/color]([COLOR=#804040][b]*[/b][/color]) :: desc
  [COLOR=#2e8b57][b]integer[/b][/color] :: n
  [COLOR=#2e8b57][b]double precision[/b][/color] :: WR([COLOR=#804040][b]*[/b][/color]), WI([COLOR=#804040][b]*[/b][/color])
  [COLOR=#2e8b57][b]integer[/b][/color] :: j
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) 
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) desc

  [COLOR=#804040][b]do[/b][/color] j[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
     [COLOR=#804040][b]if[/b][/color] (WI(j) [COLOR=#804040][b]==[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]9998[/color],[COLOR=#804040][b]advance[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]'no'[/color]) WR(j)
     [COLOR=#804040][b]else[/b][/color]
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]9999[/color],[COLOR=#804040][b]advance[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]'no'[/color]) WR(j), WI(j)
     [COLOR=#804040][b]end if[/b][/color]   
  [COLOR=#804040][b]end do[/b][/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) 

  [COLOR=#ff00ff]9998[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#ff00ff]11[/color](:,[COLOR=#008080]1X[/color],[COLOR=#008080]F6.2[/color]) )
  [COLOR=#ff00ff]9999[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#ff00ff]11[/color](:,[COLOR=#008080]1X[/color],[COLOR=#ff00ff]'('[/color],[COLOR=#008080]F6.2[/color],[COLOR=#ff00ff]','[/color],[COLOR=#008080]F6.2[/color],[COLOR=#ff00ff]')'[/color]))
[COLOR=#a020f0]end subroutine[/color] print_eigenvalues

[COLOR=#a020f0]subroutine[/color] print_eigenvectors (desc, n, WI, V, ldv)
  [COLOR=#2e8b57][b]character[/b][/color][COLOR=#804040][b]*[/b][/color]([COLOR=#804040][b]*[/b][/color]) :: desc
  [COLOR=#2e8b57][b]integer[/b][/color] :: n, ldv
  [COLOR=#2e8b57][b]double precision[/b][/color] :: WI([COLOR=#804040][b]*[/b][/color]), V(ldv, [COLOR=#804040][b]*[/b][/color])
  [COLOR=#2e8b57][b]integer[/b][/color] :: i, j
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) 
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) desc

  [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
     j [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color]
     [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] (j [COLOR=#804040][b]<=[/b][/color] n)
        [COLOR=#804040][b]if[/b][/color] (WI(j) [COLOR=#804040][b]==[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
           [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]9998[/color],[COLOR=#804040][b]advance[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]'no'[/color]) V(i,j)
           j [COLOR=#804040][b]=[/b][/color] j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
        [COLOR=#804040][b]else[/b][/color]
           [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]9999[/color],[COLOR=#804040][b]advance[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]'no'[/color]) V(i,j), V(i,j[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color])
           [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]9999[/color],[COLOR=#804040][b]advance[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]'no'[/color]) V(i,j),[COLOR=#804040][b]-[/b][/color]V(i,j[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color])
           j [COLOR=#804040][b]=[/b][/color] j [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]2[/color]
        [COLOR=#804040][b]end if[/b][/color]   
     [COLOR=#804040][b]end do[/b][/color]
     [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) 
  [COLOR=#804040][b]end do[/b][/color]

  [COLOR=#ff00ff]9998[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#ff00ff]11[/color](:,[COLOR=#008080]1X[/color],[COLOR=#008080]F6.2[/color]) )
  [COLOR=#ff00ff]9999[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#ff00ff]11[/color](:,[COLOR=#008080]1X[/color],[COLOR=#ff00ff]'('[/color],[COLOR=#008080]F6.2[/color],[COLOR=#ff00ff]','[/color],[COLOR=#008080]F6.2[/color],[COLOR=#ff00ff]')'[/color]))
[COLOR=#a020f0]end subroutine[/color] print_eigenvectors

and compiled it and linked with lapack DLL with this command

$ gfortran -o eigenvectors eigenvectors.f95 -L. -lliblapack

 

[mjerdem]

Thank you mikrom. I will try to incorporate the subroutine into my main code.

I appreciate all the help.