!************************************************************************* !* Eigenvalues and Eigenvectors of a general complex square matrix * !* using the QR method for complex matrices * !* --------------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Example #1: * !* * !* Input matrix (real symmetric for comparison with Jacobi): * !* 1, 2, 3, -7, 12 * !* 2, 4, 7, 3, -1 * !* 3, 7, 10, 8, 4 * !* -7, 3, 8, -0.75,-9 * !* 12, -1, 4, -9, 10 * !* * !* Input matrix (real part): * !*1 1.00000000 2.00000000 3.00000000 -7.00000000 12.00000000 * !*2 2.00000000 4.00000000 7.00000000 3.00000000 -1.00000000 * !*3 3.00000000 7.00000000 10.00000000 8.00000000 4.00000000 * !*4 -7.00000000 3.00000000 8.00000000 -0.75000000 -9.00000000 * !*5 12.00000000 -1.00000000 4.00000000 -9.00000000 10.00000000 * !* * !* Error code = 0 * !* * !* Eigebvalues: * !* 0.463349527981464 0.000000000000000E+000 * !* -7.77457972947830 0.000000000000000E+000 * !* -10.4865451655156 0.000000000000000E+000 * !* 18.2918206016854 0.000000000000000E+000 * !* 23.7559547653270 0.000000000000000E+000 * !* * !* Eigenvectors (real part in lines): * !*1 0.20881245 1.00000000 -0.47802717 -0.27500001 -0.21691482 * !*2 1.00000000 -0.32610016 0.20961973 -0.14768922 -0.81542193 * !*3 0.46717183 -0.00794702 -0.50782668 1.00000000 0.26443199 * !*4 0.09305260 0.59491100 1.00000000 0.45566550 0.05074045 * !*5 0.70581969 -0.01267740 0.13148643 -0.52749949 1.00000000 * !* * !* Example #2: * !* * !* Input matrix (real part): * !* 1, 5, 3, -7, 12 * !* 2, 4, 7, 3, -1 * !* 3, 7, 10, 8, 4 * !* -7, 3, 8, -0.75,-9 * !* 12, -1, 4, -9, 10 * !* Input matrix (imaginary part): * !* 1, 1, 1, 1, 1 * !* 1, 1, 1, 1, 1 * !* 1, 1, 1, 1, 1 * !* 1, 1, 1, 1 1 * !* 1, 1, 1, 1, 1 * !* * !* Input Matrix (real part): * !*1 1.00000000 5.00000000 3.00000000 -7.00000000 12.00000000 * !*2 2.00000000 4.00000000 7.00000000 3.00000000 -1.00000000 * !*3 3.00000000 7.00000000 10.00000000 8.00000000 4.00000000 * !*4 -7.00000000 3.00000000 8.00000000 -0.75000000 -9.00000000 * !*5 12.00000000 -1.00000000 4.00000000 -9.00000000 10.00000000 * !* Input Matrix (imag. part): * !*1 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 * !*2 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 * !*3 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 * !*4 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 * !*5 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 * !* * !* Error code = 0 * !* * !* Eigenvalues: * !* 0.888713348723694 3.470700022098407E-002 * !* -8.28426842306568 2.498162361886443E-002 * !* -10.3636733606943 0.879708207660510 * !* 18.8947919026401 3.15558227567471 * !* 23.1144365323962 0.905020892824929 * !* * !* Eigenvectors (real part in lines): * !*1 0.41150176 1.00000000 -0.53250763 -0.21794537 -0.41700119 * !*2 1.00000000 -0.32418914 0.26914039 -0.27758700 -0.87269351 * !*3 0.45123301 -0.00728949 -0.51526552 1.00000000 0.27273187 * !*4 0.17988445 0.60859892 1.00000000 0.45640602 0.06919679 * !*5 0.68265703 -0.10666293 -0.00445564 -0.61516974 1.00000000 * !* Eigenvectors (imag. part in lines): * !*1 0.00076996 0.00000000 -0.03474965 0.00906216 -0.00518060 * !*2 0.00000000 0.00167602 0.05162180 -0.07411949 -0.03764759 * !*3 0.04349160 -0.07014688 -0.05051509 0.00000000 -0.06636677 * !*4 -0.19851484 0.07550215 0.00000000 0.27102952 -0.31492765 * !*5 0.04407051 0.19380611 0.29215881 0.17278034 0.00000000 * !* * !* --------------------------------------------------------------------- * !* Reference: From Numath Library By Tuan Dang Trong in Fortran 77 * !* [BIBLI 18]. * !* * !* F90 Release 1.1 By J-P Moreau, Paris * !* (www.jpmoreau.fr) * !* --------------------------------------------------------------------- * !* Release 1.1: added subroutine MTPRINT and Example #2 (12/04/2004). * !************************************************************************* PROGRAM TEST_CEIGEN parameter(NM=5,N=5,LMA=5) REAL*8 AR(NM,N), AI(NM,N),ZR(LMA,N),ZI(LMA,N) REAL*8 WR(N),WI(N),WORK(3*N) LOGICAL AVEC INTEGER IERR ! For #1: DATA AR /1.d0,2.d0,3.d0,-7.d0,12.d0, & 2.d0,4.d0,7.d0,3.d0,-1.d0, & 3.d0,7.d0,10.d0,8.d0,4.d0, & -7.d0,3.d0,8.d0,-0.75d0,-9.d0, & 12.d0,-1.d0,4.d0,-9.d0,10.d0/ AI = 0.d0 ! For #2: ! DATA AR /1.d0,2.d0,3.d0,-7.d0,12.d0, & ! 5.d0,4.d0,7.d0,3.d0,-1.d0, & ! 3.d0,7.d0,10.d0,8.d0,4.d0, & ! -7.d0,3.d0,8.d0,-0.75d0,-9.d0, & ! 12.d0,-1.d0,4.d0,-9.d0,10.d0/ ! AI = 1.d0 CALL MPRINT(N,NM,N,' Input Matrix (real part):',AR) ! #2 only: ! CALL MPRINT(N,NM,N,' Input Matrix (imag. part):',AI) ! pause CALL CEIGEN(AR,AI,NM,N,.TRUE.,WR,WI,ZR,ZI,WORK,IERR) print *,' ' print *,' Error code =', IERR print *,' ' print *,' Eigenvalues:' do i=1,N print *, WR(i), WI(i) end do print *,' ' CALL MTPRINT(N,NM,N,' Eigenvectors (real part in lines):',ZR) ! #2 only: ! CALL MTPRINT(N,NM,N,' Eigenvectors (imag. part in lines):',ZI) print *,' ' END !of main program !print a matrix with caption SUBROUTINE MPRINT(N,MMAX,M,LABEL,A) REAL *8 A(MMAX,1) CHARACTER*(*) LABEL PRINT 2,LABEL DO 1 I=1,N 1 PRINT 3,I,(A(I,J),J=1,M) RETURN 2 FORMAT(A) 3 FORMAT(I1,5(F13.8)) END !print transpose matrix with caption SUBROUTINE MTPRINT(N,MMAX,M,LABEL,A) REAL *8 A(MMAX,1) CHARACTER*(*) LABEL PRINT 2,LABEL DO 1 I=1,N 1 PRINT 3,I,(A(J,I),J=1,M) RETURN 2 FORMAT(A) 3 FORMAT(I1,5(F13.8)) END SUBROUTINE CEIGEN(AR,AI,NM,N,AVEC,WR,WI,ZR,ZI,WORK,IERR) ! -------------------------------------------------------------------- ! This subroutine calculates the eigenvalues and eigenvectors ! of a general complex matrix. The used method begins with a ! reduction by unitary transformations to a Hessenberg form. ! The QR algorithm is then applied to get the eigenvalues. ! Finally, the corresponding eigenvectors are found by a simple ! substitution. ! -------------------------------------------------------------------- ! CALLING MODE: ! CALL CEIGEN (AR,AI,LMA,N,AVEC,WR,WI,ZR,ZI,WORK,IERR) ! INPUTS: ! AR = TABLE(LMA,N) OF REAL COEF. OF MATRIX A R*8 ! AI = TABLE(LMA,N) OF IMAGINARY COEF. OF MATRIX A R*8 ! NM = 1ST DIMENSION OF MATRICES AR,AI,ZR & ZI ! AS DECLARED IN CALLING PROGRAM I*4 ! N = 2ND DIMENSION OF MATRICES AR,AI,ZR & ZI I*4 ! AVEC = CODE TO CALCULATE OR NOT EIGENVECTORS L*8 ! = .FALSE. NO EIGENVECTORS ! = .TRUE. WITH EIGENVECTORS ! OUTPUTS: ! WR,WI= TABLEA(N) OF EIGENVALUES R*8 ! ZR,ZI= TABLE(LMA,N) OF EIGENVECTORS R*8 ! IERR = ERROR CODE I*4 ! = 0 NO ERROR ! < OU = N NO CONVERGENCE ! = 9999 DIMENSION ERROR IN TABLES A OR ZR,ZI ! WORKING SPACE: ! WORK = TABLE(3*N) R*8 ! REFERENCE: ! J.H.WILKINSON & C.REINSCH ! LINEAR ALGEBRA SPRINGER-VERLAG 1971 ! -------------------------------------------------------------------- REAL*8 AR(NM,*),AI(NM,*),ZR(NM,*),ZI(NM,*),WR(*),WI(*),WORK(*) INTEGER UPP LOGICAL AVEC IF(N.LE.NM) THEN N1=N+1 N2=N+N1 CALL CBALAN(NM,N,AR,AI,LOW,UPP,WORK(1)) CALL CUNITH(NM,N,LOW,UPP,AR,AI,WORK(N1),WORK(N2)) CALL COMQR2(NM,N,LOW,UPP,WORK(N1),WORK(N2),AR,AI,WR,WI, & AVEC,ZR,ZI,IERR) IF(IERR.EQ.0) THEN CALL CBALBK(NM,N,LOW,UPP,N,WORK(1),ZR,ZI) CALL NORMAL(NM,N,WR,WI,AVEC,ZR,ZI,WORK(N1),WORK(N2)) ENDIF ELSE IERR= 9999 ENDIF RETURN END ! SUBROUTINE CUNITH(NM,N,LOW,UPP,AR,AI,DR,DI) ! UNITARY TRANSFORMATIONS TO REDUCE A COMPLEX MATRIX ! TO A HESSENBERG FORM REAL*8 AR(NM,*),AI(NM,*),DR(*),DI(*),XR,XI,G,F,H,S INTEGER UPP IF(UPP-1.LT.LOW+1) RETURN DO M=LOW+1,UPP-1 DR(M)=0.D0 DI(M)=0.D0 S=0.D0 H=0.D0 DO I=M,UPP S=S+ABS(AR(I,M-1))+ABS(AI(I,M-1)) ENDDO ! WRITE(*,'(1X,A,I4,E15.6)') ' M,S ',M,S IF(S.NE.0.D0) THEN DO I=UPP,M,-1 DR(I)=AR(I,M-1)/S DI(I)=AI(I,M-1)/S H=H+DR(I)*DR(I)+DI(I)*DI(I) ENDDO G=SQRT(H) F=ABSCD(DR(M),DI(M)) IF(F.EQ.0.D0) THEN DR(M)=G AR(M,M-1)=S ELSE H=H+F*G G=G/F DR(M)=(1.D0+G)*DR(M) DI(M)=(1.D0+G)*DI(M) ENDIF DO J=M,N XR=0.D0 XI=0.D0 DO I=UPP,M,-1 XR=XR+DR(I)*AR(I,J)+DI(I)*AI(I,J) XI=XI+DR(I)*AI(I,J)-DI(I)*AR(I,J) ENDDO XR=XR/H XI=XI/H DO I=M,UPP AR(I,J)=AR(I,J)-XR*DR(I)+XI*DI(I) AI(I,J)=AI(I,J)-XR*DI(I)-XI*DR(I) ENDDO ENDDO DO I=1,UPP XR=0.D0 XI=0.D0 DO J=UPP,M,-1 XR=XR+DR(J)*AR(I,J)-DI(J)*AI(I,J) XI=XI+DR(J)*AI(I,J)+DI(J)*AR(I,J) ENDDO XR=XR/H XI=XI/H DO J=M,UPP AR(I,J)=AR(I,J)-XR*DR(J)-XI*DI(J) AI(I,J)=AI(I,J)+XR*DI(J)-XI*DR(J) ENDDO ENDDO DR(M)=DR(M)*S DI(M)=DI(M)*S AR(M,M-1)=-G*AR(M,M-1) AI(M,M-1)=-G*AI(M,M-1) ENDIF ENDDO RETURN END SUBROUTINE COMQR2(NM,N,LOW,UPP,GR,GI,HR,HI,WR,WI,AVEC,ZR,ZI,IERR) ! --------------------------------------------------------------------- ! QR ALGORITH APPLIED TO A COMPLEX HESSENBERG MATRIX ! --------------------------------------------------------------------- REAL*8 HR(NM,*),HI(NM,*),ZR(NM,*),ZI(NM,*),WR(*),WI(*),GR(*), & GI(*),SR,SI,TR,TI,YR,YI,XR,XI,VR,VI,MACHEP,NORM,EPS,EPS1 INTEGER I,J,K,L,M,N,IN,NA,NM,ITS,LOW,UPP,IERR LOGICAL AVEC ! SEEK FLOATING POINT RELATIVE PRECISION FOR TYPE REAL*8 EPS=1.D0 5 EPS=EPS*0.5D0 EPS1=1.D0+EPS IF(EPS1.GT.1.D0) GO TO 5 MACHEP=2.D0*EPS IERR=0 IF(AVEC) THEN DO I=1,N DO J=1,N ZR(I,J)=0.D0 ZI(I,J)=0.D0 ENDDO ZR(I,I)=1.D0 ENDDO ! ACCUMULATE UNITARY TRANSFORMATIONS DO I=UPP-1,LOW+1,-1 IF((HR(I,I-1).EQ.0.D0.AND.HI(I,I-1).EQ.0.D0) & .OR.(GR(I).EQ.0.D0.AND.GI(I).EQ.0.D0)) GO TO 7 NORM=HR(I,I-1)*GR(I)+HI(I,I-1)*GI(I) DO K=I+1,UPP GR(K)=HR(K,I-1) GI(K)=HI(K,I-1) ENDDO DO J=I,UPP SR=0.D0 SI=0.D0 DO K=I,UPP SR=SR+GR(K)*ZR(K,J)+GI(K)*ZI(K,J) SI=SI+GR(K)*ZI(K,J)-GI(K)*ZR(K,J) ENDDO SR=SR/NORM SI=SI/NORM DO K=I,UPP ZR(K,J)=ZR(K,J)+SR*GR(K)-SI*GI(K) ZI(K,J)=ZI(K,J)+SR*GI(K)+SI*GR(K) ENDDO ENDDO 7 CONTINUE ENDDO ENDIF ! CREATE REAL SUBDIAGONAL L=LOW+1 DO I=L,UPP M=MIN(I+1,UPP) IF(HI(I,I-1).EQ.0.D0) GO TO 9 NORM=ABSCD(HR(I,I-1),HI(I,I-1)) YR=HR(I,I-1)/NORM YI=HI(I,I-1)/NORM HR(I,I-1)=NORM HI(I,I-1)=0.D0 DO J=I,N SI= YR*HI(I,J)-YI*HR(I,J) HR(I,J)=YR*HR(I,J)+YI*HI(I,J) HI(I,J)=SI ENDDO DO J=1,M SI= YR*HI(J,I)+YI*HR(J,I) HR(J,I)=YR*HR(J,I)-YI*HI(J,I) HI(J,I)=SI ENDDO IF(AVEC) THEN DO J=LOW,UPP SI= YR*ZI(J,I)+YI*ZR(J,I) ZR(J,I)=YR*ZR(J,I)-YI*ZI(J,I) ZI(J,I)=SI ENDDO ENDIF 9 CONTINUE ENDDO ! STORE ROOTS ALREADY ISOLATED BY SUBROUTINE BALANC DO I=1,N IF(I.LT.LOW.OR.I.GT.UPP) THEN WR(I)=HR(I,I) WI(I)=HI(I,I) ENDIF ENDDO IN=UPP TR=0.D0 TI=0.D0 ! SEEK EIGENVALUES 10 IF(IN.LT.LOW) GO TO 50 ITS=0 NA=IN-1 ! WE LOOK AT REAL PART OF SUBDIAGONAL ELEMENT 15 DO L=IN,LOW,-1 SR=MACHEP*(ABS(HR(L-1,L-1))+ABS(HI(L-1,L-1)) & +ABS(HR(L,L))+ABS(HI(L,L))) IF (ABS(HR(L,L-1)).LE.SR.OR.L.EQ.LOW) GO TO 20 ENDDO ! SHIFTING 20 IF(L.EQ.IN) GO TO 45 IF(ITS.EQ.30) GO TO 55 IF(ITS.EQ.10.OR.ITS.EQ.20) THEN ! EXTRA SHIFTING SR=ABS(HR(IN,NA))+ABS(HR(NA,IN-2)) SI=0.D0 ELSE SR=HR(IN,IN) SI=HI(IN,IN) XR=HR(NA,IN)*HR(IN,NA) XI=HI(NA,IN)*HR(IN,NA) IF(XR.NE.0.D0.OR.XI.NE.0.D0) THEN YR=(HR(NA,NA)-SR)*.5D0 YI=(HI(NA,NA)-SI)*.5D0 CALL DCSQRT(YR*YR-YI*YI+XR,2.D0*YR*YI+XI,VR,VI) IF(YR*VR+YI*VI.LT.0.D0) THEN VR=-VR VI=-VI ENDIF CALL DIVCD(XR,XI,YR+VR,YI+VI,VR,VI) SR=SR-VR SI=SI-VI ENDIF ENDIF DO I=LOW,IN HR(I,I)=HR(I,I)-SR HI(I,I)=HI(I,I)-SI ENDDO TR=TR+SR TI=TI+SI ITS=ITS+1 ! TRIANGULARY DECOMPOSITION H = Q * R DO I=L+1,IN SR=HR(I,I-1) VR=HR(I-1,I-1) VI=HI(I-1,I-1) NORM=SQRT(VR*VR+VI*VI+SR*SR) XR=VR/NORM XI=VI/NORM WR(I-1)=XR WI(I-1)=XI HR(I-1,I-1)=NORM HI(I-1,I-1)=0.D0 HR(I,I-1)=0.D0 HI(I,I-1)=SR/NORM DO J=I,N YR=HR(I-1,J) YI=HI(I-1,J) VR=HR(I,J) VI=HI(I,J) HR(I-1,J)=XR*YR+XI*YI+HI(I,I-1)*VR HI(I-1,J)=XR*YI-XI*YR+HI(I,I-1)*VI HR(I,J) =XR*VR-XI*VI-HI(I,I-1)*YR HI(I,J) =XR*VI+XI*VR-HI(I,I-1)*YI ENDDO ENDDO SI=HI(IN,IN) IF(SI.EQ.0.D0) GO TO 35 NORM=ABSCD(HR(IN,IN),SI) SR=HR(IN,IN)/NORM SI=SI/NORM HI(IN,IN)=0.D0 HR(IN,IN)=NORM IF(IN.NE.N) THEN DO J=IN+1,N YR=HR(IN,J) YI=HI(IN,J) HR(IN,J)=SR*YR+SI*YI HI(IN,J)=SR*YI-SI*YR ENDDO ENDIF ! INVERSE OPERATION 35 DO J=L+1,IN XR=WR(J-1) XI=WI(J-1) DO I=1,J YR=HR(I,J-1) YI=0.D0 VR=HR(I,J) VI=HI(I,J) IF(I.NE.J) THEN YI=HI(I,J-1) HI(I,J-1)=XR*YI+XI*YR+HI(J,J-1)*VI ENDIF HR(I,J-1)=XR*YR-XI*YI+HI(J,J-1)*VR HR(I,J) =XR*VR+XI*VI-HI(J,J-1)*YR HI(I,J) =XR*VI-XI*VR-HI(J,J-1)*YI ENDDO IF(AVEC) THEN DO I=LOW,UPP YR=ZR(I,J-1) YI=ZI(I,J-1) VR=ZR(I,J) VI=ZI(I,J) ZR(I,J-1)=XR*YR-XI*YI+HI(J,J-1)*VR ZI(I,J-1)=XR*YI+XI*YR+HI(J,J-1)*VI ZR(I,J) =XR*VR+XI*VI-HI(J,J-1)*YR ZI(I,J) =XR*VI-XI*VR-HI(J,J-1)*YI ENDDO ENDIF ENDDO IF(SI.EQ.0.D0) GO TO 15 DO I=1,IN YR=HR(I,IN) YI=HI(I,IN) HR(I,IN)=SR*YR-SI*YI HI(I,IN)=SR*YI+SI*YR ENDDO IF(AVEC) THEN DO I=LOW,UPP YR=ZR(I,IN) YI=ZI(I,IN) ZR(I,IN)=SR*YR-SI*YI ZI(I,IN)=SR*YI+SI*YR ENDDO ENDIF GO TO 15 ! ONE ROOT HAS BEEN FOUND 45 HR(IN,IN)=HR(IN,IN)+TR HI(IN,IN)=HI(IN,IN)+TI WR(IN)=HR(IN,IN) WI(IN)=HI(IN,IN) IN=NA GO TO 10 ! ALL ROOTS HAVE BEEN FOUND ! START CALCULATION OF EIGENVECTORS BY SUBSTITUTION 50 IF(AVEC) THEN NORM=0.D0 DO I=1,N DO J=I,N NORM=NORM+ABS(HR(I,J))+ABS(HI(I,J)) ENDDO ENDDO IF(NORM.EQ.0.D0) GO TO 60 DO IN=N,2,-1 XR=WR(IN) XI=WI(IN) HR(IN,IN)=1.D0 HI(IN,IN)=0.D0 DO I=IN-1,1,-1 VR=HR(I,IN) VI=HI(I,IN) IF(I.NE.IN-1) THEN DO J=I+1,IN-1 VR=VR+HR(I,J)*HR(J,IN)-HI(I,J)*HI(J,IN) VI=VI+HR(I,J)*HI(J,IN)+HI(I,J)*HR(J,IN) ENDDO ENDIF YR=XR-WR(I) YI=XI-WI(I) IF(YR.EQ.0.D0.AND.YI.EQ.0.D0) YR=NORM*MACHEP CALL DIVCD(VR,VI,YR,YI,TR,TI) HR(I,IN)=TR HI(I,IN)=TI ENDDO ENDDO DO I=1,N-1 IF(I.LT.LOW.OR.I.GT.UPP) THEN DO J=I+1,N ZR(I,J)=HR(I,J) ZI(I,J)=HI(I,J) ENDDO ENDIF ENDDO ! END CALCULATION OH HESSENBERG EIGENVECTORS DO J=N,LOW+1,-1 M=MIN(J-1,UPP) DO I=LOW,UPP VR=ZR(I,J) VI=ZI(I,J) DO K=LOW,M VR=VR+ZR(I,K)*HR(K,J)-ZI(I,K)*HI(K,J) VI=VI+ZR(I,K)*HI(K,J)+ZI(I,K)*HR(K,J) ENDDO ZR(I,J)=VR ZI(I,J)=VI ENDDO ENDDO ENDIF ! END CALCULATION OF EIGENVECTORS OF INPUT MATRIX GO TO 60 ! NO CONVERGENCE AFTER 30 ITERATIONS! 55 IERR = IN 60 RETURN END SUBROUTINE CBALAN(NM,N,AR,AI,LOW,UPP,D) ! BALANCE A COMPLEX MATRIX (VERSION DERIVED FROM BALANC) REAL*8 AR(NM,N),AI(NM,N),D(N),B2,C,F,G,R,S LOGICAL NOCONV INTEGER UPP,B DATA B/16/ B2=B*B L=1 K=N 10 DO J=K,1,-1 J1=J K1=K L1=L DO I=1,K IF(I.NE.J) THEN IF(AR(J,I).NE.0.D0 .OR. AI(J,I).NE.0.D0) GO TO 11 ENDIF ENDDO CALL EXCHG1(K,NM,N,AR,AI,D,J1,K1,L1) K=K-1 GO TO 10 11 CONTINUE ENDDO 20 DO J=L,K J1=J K1=K L1=L DO I=L,K IF(I.NE.J) THEN IF(AR(I,J).NE.0.D0 .OR. AI(I,J).NE.0.D0) GO TO 22 ENDIF ENDDO CALL EXCHG1(L,NM,N,AR,AI,D,J1,K1,L1) L=L+1 GO TO 20 22 CONTINUE ENDDO LOW=L UPP=K DO I=L,K D(I)=1.D0 ENDDO 30 NOCONV=.FALSE. DO I=L,K C=0.D0 R=0.D0 DO J=L,K IF(J.NE.I) THEN C=C+ABS(AR(J,I))+ABS(AI(J,I)) R=R+ABS(AR(I,J))+ABS(AI(I,J)) ENDIF ENDDO G=R/B F=1.D0 S=C+R 40 IF(C.LT.G) THEN F=F*B C=C*B2 GO TO 40 ENDIF G=R*B 50 IF(C.GE.G) THEN F=F/B C=C/B2 GO TO 50 ENDIF IF((C+R)/F.LT.0.95D0*S) THEN G=1.D0/F D(I)=D(I)*F NOCONV=.TRUE. DO J=L,N AR(I,J)=AR(I,J)*G AI(I,J)=AI(I,J)*G ENDDO DO J=1,K AI(J,I)=AI(J,I)*F AR(J,I)=AR(J,I)*F ENDDO ENDIF ENDDO IF(NOCONV) GO TO 30 RETURN END SUBROUTINE CBALBK(NM,N,LOW,UPP,M,D,Z,T) ! This subroutine finds the eigenvectors of a complex matrix, ! as output of subroutine CBALAN. REAL*8 Z(NM,*),T(NM,*),D(*),S INTEGER UPP DO I=LOW,UPP S=D(I) DO J=1,M Z(I,J)=Z(I,J)*S T(I,J)=T(I,J)*S ENDDO ENDDO DO L=1,N I=L IF(I.GE.LOW.AND.I.LE.UPP) GO TO 10 IF(I.LT.LOW) I=LOW-L K=D(I) IF(K.NE.I) THEN DO J=1,M S=Z(I,J) Z(I,J)=Z(K,J) Z(K,J)=S S=T(I,J) T(I,J)=T(K,J) T(K,J)=S ENDDO ENDIF 10 CONTINUE ENDDO RETURN END SUBROUTINE NORMAL(NM,N,WR,WI,AVEC,ZR,ZI,TR,TI) ! Sort eigenvalues by absolute values in ascending order and ! normalize. REAL*8 ZR(NM,*),ZI(NM,*),WR(*),WI(*),TR(*),TI(*),ZM,Z1,Z2,VR,VI, & ABSCD LOGICAL AVEC ! SORT SOLUTIONS IN ASCENDING ORDER DO J=2,N VR=WR(J) VI=WI(J) IF(AVEC) THEN DO K=1,N TR(K)=ZR(K,J) TI(K)=ZI(K,J) ENDDO ENDIF DO I=J-1,1,-1 IF( ABSCD(WR(I),WI(I)).LE. ABSCD(VR,VI) ) GO TO 5 WR(I+1)=WR(I) WI(I+1)=WI(I) IF(AVEC) THEN DO K=1,N ZR(K,I+1)=ZR(K,I) ZI(K,I+1)=ZI(K,I) ENDDO ENDIF ENDDO I=0 5 WR(I+1)=VR WI(I+1)=VI IF(AVEC) THEN DO K=1,N ZR(K,I+1)=TR(K) ZI(K,I+1)=TI(K) ENDDO ENDIF ENDDO ! NORMALIZE WITH RESPECT TO BIGGEST ELEMENT IF(AVEC) THEN DO J=N,1,-1 ZM= 0.D0 DO I=1,N VR=ZR(I,J) VI=ZI(I,J) Z1=ABSCD (VR,VI) IF(Z1.GE.ZM) THEN IM = I ZM = Z1 ENDIF ENDDO Z1 = ZR(IM,J) Z2 = ZI(IM,J) DO I = 1,N CALL DIVCD(ZR(I,J),ZI(I,J),Z1,Z2,VR,VI) ZR(I,J) = VR ZI(I,J) = VI ENDDO ENDDO ENDIF END ! Utility routines FUNCTION ABSCD(AR,AI) ! ABSOLUTE VALUE OF A COMPLEX NUMBER C=AR+I*AI ! ABSCD,AR,AI REAL*8 ! ABSCD=SQRT(AR**2+AI**2) REAL*8 ABSCD,AR,AI,XR,XI,W XR=ABS(AR) XI=ABS(AI) IF(XR.LE.XI) THEN W=XR XR=XI XI=W ENDIF IF(XI.EQ.0.D0) THEN ABSCD=XR ELSE ABSCD=XR*SQRT(1.D0+(XI/XR)**2) ENDIF RETURN END SUBROUTINE DCSQRT(X,Y,A,B) ! SQUARE ROOT OF A COMPLEX NUMBER A+I*B = RAC(X+I*Y) REAL*8 X,Y,A,B,DCABS IF(X.EQ.0.D0.AND.Y.EQ.0.D0) THEN A=0.D0 B=0.D0 ELSE A=SQRT((ABS(X)+DCABS(X,Y))*0.5D0) IF(X.GE.0.D0) THEN B=Y/(A+A) ELSE IF(Y.LT.0.D0) THEN B=-A ELSE B=A A=Y/(B+B) ENDIF ENDIF ENDIF END ! ABSOLUTE VALUE OF A COMPLEX NUMBER X+I*Y FUNCTION DCABS(XX,YY) REAL*8 DCABS,XX,YY,X,Y,W X=ABS(XX) Y=ABS(YY) IF(X.EQ.0.D0) THEN W=Y ELSE IF(Y.EQ.0) THEN W=X ELSE IF(X.GT.Y) THEN W=X*SQRT(1.D0+(Y/X)**2) ELSE W=Y*SQRT(1.D0+(X/Y)**2) ENDIF ENDIF ENDIF DCABS=W RETURN END SUBROUTINE DIVCD(AR,AI,BR,BI,ZR,ZI) ! EFFECT DIVISION Z=(ZR+I*ZI)=(AR+I*AI)/(BR+I*BI) !-----DO NOT USE IF BR=BI=0. REAL*8 AR,AI,BR,BI,YR,YI,ZR,ZI,W YR=BR YI=BI IF(ABS(YR).GT.ABS(YI)) THEN W=YI/YR YR=W*YI+YR ZR=(AR+W*AI)/YR ZI=(AI-W*AR)/YR ELSE W=YR/YI YI=W*YR+YI ZR=(W*AR+AI)/YI ZI=(W*AI-AR)/YI ENDIF RETURN END ! Exchange elements of a complex matrix SUBROUTINE EXCHG1(M,NM,N,A,B,D,J,K,L) REAL*8 A(NM,N),B(NM,N),D(N),F D(M)=J IF(J.NE.M) THEN DO I=1,K F=A(I,J) A(I,J)=A(I,M) A(I,M)=F F=B(I,J) B(I,J)=B(I,M) B(I,M)=F ENDDO DO I=L,N F=A(J,I) A(J,I)=A(M,I) A(M,I)=F F=B(J,I) B(J,I)=B(M,I) B(M,I)=F ENDDO ENDIF RETURN END !end of file tceigen.f90