ZGGEVX bug - wrong eigenvectors

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

ZGGEVX bug - wrong eigenvectors

Postby Hopito » Sun May 14, 2017 6:10 am

Hello,

I am using the subroutine ZGGEVX
(it computes for a pair of n by n complex nonsymmetric matrices (A,B) the generalized eigenvalues and, optionally, the left and/or right generalized eigenvectors using the QZ algorithm). Even though the eigenvalues are correct, all the eigenvectors (left and right) but one are wrong. I used the subroutine on 3 different probles (2x2 complex matrices (A,B), 3x3 complex matrices (A,B),, 4x4 complex matrices (A,B) and for all three problems the eigenvalues were correct and the last eigenvector (i.e. for a problem with 2x2 complex matrices the second eigenvector and for 4x4 matrices the fourth eigenvector) that is given by the subroutine is correct.
I have checked the following 1. Whether the workspace is incorrect, I gave on purpose a very large value (nmax). Same result.
2. Whether the matrices are being read incorrectly from the txt file. They are read correct.
3. Whether the type of the variables specified is wrong. It is not. (at least i do not think so)
4. Whether something goes wrong on the storage of the eigenvectors and the components of them are being messed up. It's
not that either.
5. Whether the mistake is that i define complex(kind=8) i changed it to complex*16. Same result. (i am using microsoft visual studio 2010 - fortran90)
I cannot find the mistake. Please take a look to see if you can find something which might be obvious to you but not for me.

program MyProgramI
IMPLICIT NONE
! Variables
INTEGER nin ! Name of INPUT.TXT file and OUTPUT.TXT FILES
PARAMETER (nin=20)
INTEGER nb, nmax ! Number of common blocks, number of rows
PARAMETER (nb=64,nmax=10)
INTEGER lda, ldb,ldvr, lwork
PARAMETER (lda=nmax,ldvr=nmax, ldb=nmax,lwork =nb*nmax+2*nmax*nmax)
!*************** LOCAL SCALARS **********************************************************
DOUBLE PRECISION abnorm, abnrm, bbnrm, eps, erbnd,rcnd, small, tol
INTEGER i, Ihi,Ilo, info, j, lwkopt, n
!**************** LOCAL ARRAYS
COMPLEX(KIND=8) a(lda,nmax), alpha(nmax), b(ldb,nmax),beta(nmax), dummy(1,1) ,vr(ldvr, nmax), work(lwork)
DOUBLE PRECISION lscale(nmax), rconde(nmax), rcondv(nmax), rscale(nmax), rwork(6*nmax)
INTEGER iwork(nmax+2)
LOGICAL bwork(nmax)
!******** EXTERNAL FUNCTIONS
DOUBLE PRECISION dlamch
EXTERNAL dlamch
EXTERNAL ZGGEVX
!******** INTRINSIC FUNCTIONS
INTRINSIC abs,idint,sqrt
!****** EXECUTABLE STATEMENTS
WRITE(*,*) 'ZGGEVX EXAMPLE PROGRAM RESULTS'
! Skip heading in data file
OPEN(nin,file='C:\Users\Elpida\Desktop\data3.txt',STATUS = 'OLD',ACTION = 'READ')
READ(nin,*)
READ(nin,*) n
IF (n<=nmax) THEN
READ(nin,*) ((a(i,j),j=1,n),i=1,n)
READ(nin,*) ((b(i,j),j=1,n),i=1,n)


! SOLVE THE GENERALIZED EIGENVALUE PROBLEM
CALL zggevx('B','N','V','B',n,a,lda,b,ldb,alpha,beta,dummy,1 ,vr,ldvr,Ilo,Ihi,lscale,rscale,abnrm,bbnrm,rconde,rcondv,work,lwork,rwork,iwork,bwork,info)
DO i=1,n
DO j=1,n
Write(*,*) a(i,j)
end do
end do
IF (info>0) THEN
WRITE(*,*)
WRITE(*,100) 'Failure in zggevx. INFO = ', info
ELSE
!Compute the machine precision, the safe range parameter
eps = dlamch('Eps')
small = dlamch('Sfmin')
abnorm = (sqrt(abnrm**2+ bbnrm**2))
tol = eps*abnorm

! Print out the eigenvalues and vectors and associated condition number and bounds

DO j=1,n
WRITE(*,*)
IF ((abs(alpha(j)))*small>abs(beta(j))) THEN
WRITE(*,110) 'Eigenvalue(',j,')',' is numerically infinite or undetermined','ALPHA(',j,')= ',alpha(j),', BETA(',j,')=', beta(j)
ELSE
WRITE(*,120) 'Eigenvalue(',j,') = ', alpha(j)/beta(j)
END IF
rcnd = rconde(j)
WRITE(*,*)
WRITE(*,130) 'Reciprocal condition number = ', rcnd
IF (rcnd>0.0D0) THEN
erbnd=tol/rcnd
WRITE(*,130) 'Error bound = ',erbnd
ELSE
WRITE(*,*) 'Error bound is infinite'
END IF

! PRINT OUT INFORMATION ON THE JTH EIGENVECTOR

WRITE(*,*)
WRITE(*,140) 'Eigenvector(',j,')',(vr(i,j),i=1,n)
rcnd=rcondv(j)
WRITE(*,*)
WRITE(*,130) 'Reciprocal condition number = ',rcnd
IF (rcnd>0.0D0) THEN
erbnd = tol/rcnd
WRITE(*,130) 'Error bound is = ',erbnd
ELSE
WRITE(*,*) 'Error bound is infinite'
END IF
END DO
write(*,*) work(1),lwork,lwkopt
lwkopt=(work(1))
IF (lwork<=lwkopt) THEN
WRITE(*,*)
WRITE(*,150) 'Optimum workspace required = ',lwkopt, 'Workspace provided = ',lwork
END IF
END IF
ELSE
WRITE(*,*)
WRITE(*,*) 'Nmax too small'
END IF
CLOSE(nin)
STOP

100 FORMAT(1X, A, I4)
110 FORMAT(1X, A, I2, 2A, /1X, 2(A,I2,A,'(',1P,E11.4,',',1P,E11.4,')'))
120 FORMAT(1X, A, I2, A, /1X, '(',1P,E11.4,',',1P,E11.4,')')
130 FORMAT(1X, A, 1P, E8.1)
140 FORMAT(1X, A, I2, A, /3(1X, '(',1P,E11.4,',',1P,E11.4,')',:))
150 FORMAT(1X, A, I5, /1X, A, I5)


end program MyProgramI

TXT FILE

3
( -1.00, 0.00) ( 23.00, -5.00) ( -6.15, -3.00)
(-15.00, 3.00) ( 2.00, -7.50) ( -9.20, 43.00)
( 19.00,-66.00) (-89.00, -9.00) ( 7.15, 62.00) : End of A
( 31.00, 33.00) ( 83.00,-90.00) ( 16.15, 13.00)
(-15.00, 3.00) ( 2.00, -7.50) (-39.40, 29.00)
( 13.00,-53.00) ( -3.00, -9.45) ( -7.15, 2.00) : End of B
Hopito
 
Posts: 1
Joined: Sun May 14, 2017 5:48 am

Return to Algorithm / Data

Who is online

Users browsing this forum: No registered users and 2 guests