'************************************************************************* '* 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.000000 2.000000 3.000000 -7.000000 12.000000 * '* 2.000000 4.000000 7.000000 3.000000 -1.000000 * '* 3.000000 7.000000 10.000000 8.000000 4.000000 * '* -7.000000 3.000000 8.000000 -0.750000 -9.000000 * '* 12.000000 -1.000000 4.000000 -9.000000 10.000000 * '* * '* Error code = 0 * '* * '* Eigebvalues: * '* .4633495279814548 0 * '* -7.774579729478302 0 * '* -10.48654516551558 0 * '* 18.29182060168545 0 * '* 23.75595476532698 0 * '* * '* Eigenvectors (real part in lines): * '* 0.208812 1.000000 -0.478027 -0.275000 -0.216915 * '* 1.000000 -0.326100 0.209620 -0.147689 -0.815422 * '* 0.467172 -0.007947 -0.507827 1.000000 0.264432 * '* 0.093053 0.594911 1.000000 0.455666 0.050740 * '* 0.705820 -0.012677 0.131486 -0.527500 1.000000 * '* * '* 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.000000 5.000000 3.000000 -7.000000 12.000000 * '* 2.000000 4.000000 7.000000 3.000000 -1.000000 * '* 3.000000 7.000000 10.000000 8.000000 4.000000 * '* -7.000000 3.000000 8.000000 -0.750000 -9.000000 * '* 12.000000 -1.000000 4.000000 -9.000000 10.000000 * '* Input Matrix (imag. part): * '* 1.000000 1.000000 1.000000 1.000000 1.000000 * '* 1.000000 1.000000 1.000000 1.000000 1.000000 * '* 1.000000 1.000000 1.000000 1.000000 1.000000 * '* 1.000000 1.000000 1.000000 1.000000 1.000000 * '* 1.000000 1.000000 1.000000 1.000000 1.000000 * '* * '* Error code = 0 * '* * '* Eigenvalues: * '* .8887133487236936 3.470700022098407D-02 * '* -8.284268423065674 2.498162361886612D-02 * '* -10.36367336069426 .8797082076605086 * '* 18.89479190264009 3.155582275674704 * '* 23.11443653239616 .9050208928249344 * '* * '* Eigenvectors (real part in lines): * '* 0.411502 1.000000 -0.532508 -0.217945 -0.417001 * '* 1.000000 -0.324189 0.269140 -0.277587 -0.872694 * '* 0.451233 -0.007290 -0.515266 1.000000 0.272732 * '* 0.179884 0.608600 1.000000 0.456406 0.069197 * '* 0.682657 -0.106663 -0.004456 -0.615170 1.000000 * '* Eigenvectors (imag. part in lines): * '* 0.000770 0.000000 -0.034750 0.009062 -0.005181 * '* 0.000000 0.001676 0.051622 -0.074120 -0.037648 * '* 0.043492 -0.070147 -0.050515 0.000000 -0.066367 * '* -0.198515 0.075502 0.000000 0.271030 -0.314928 * '* 0.044071 0.193806 0.292159 0.172780 0.000000 * '* * '* --------------------------------------------------------------------- * '* Reference: From Numath Library By Tuan Dang Trong in Fortran 77 * '* [BIBLI 18]. * '* * '* Basic Release 1.1 By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* --------------------------------------------------------------------- * '* Release 1.1: added subroutines 600/601 and Example #2 (12/06/2004). * '************************************************************************* ' LIST OF SUBROUTINES: ' 500: print matrix AR with caption ' 501: print matrix AI with caption ' 600: print transpose of matrix ZR with caption ' 601: print transpose of matrix ZI with caption '1000: CEIGEN (calculate eigenvalues/eigenvectors ' of a general complex matrix) '2000: CBALAN (balance a complex matrix) '3000: CUNITH (Unitary transformations to reduce a complex matrix ' to a Hessenberg form) '4000: COMQR2 (Apply QR algorithm to a complex matrix as output ' of 3000) '5000: CBALBK (Find the eigenvectors of a complex matrix, ' as output of subroutine COMQR2) '6000: NORMAL (Sort eigenvalues by absolute values in ascending ' order and normalize eigenvectors) '7000: EXCHG1 (Exchange elements of a complex matrix) '7500: ABSCD (Absolute value of a complex number AR+I*AI '7600: DCSQRT (Square root of a complex number A+I*B = SQR(X+I*Y)) '7700: DIVCD (Effect division (ZR+I*ZI)=(AR+I*AI)/(BR+I*BI)) ' ---- DO NOT USE IF BR=BI=0. 'PROGRAM TEST_CEIGEN DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 NM = 100 'Maximum size for N N = 5 'Size of complex square matrix DIM AR(N, N), AI(N, N) 'Real, Imag. parts of complex matrix DIM ZR(N, N), ZI(N, N) 'Real, Imag. parts of eigenvectors DIM WR(N), WI(N) 'Real, Imag. parts of eigenvalues DIM WORK(N), TMP1(N), TMP2(N) 'Real working vectors 'LOGICAL IAVEC (=0 NO EIGENVECTORS, =1 WITH EIGENVECTORS) 'INTEGER IERR (MUST BE ZERO) F\$ = "####.######" FOR I = 1 TO N FOR J = 1 TO N READ AR(I, J) AI(I, J) = 0# NEXT J NEXT I ' For #1: DATA 1#, 2#, 3#, -7#, 12# DATA 2#, 4#, 7#, 3#, -1# DATA 3#, 7#, 10#, 8#, 4# DATA -7#, 3#, 8#, -0.75, -9# DATA 12#, -1#, 4#, -9#, 10# ' For #2: ' Data 1#, 5#, 3#, -7#, 12# ' Data 2#, 4#, 7#, 3#, -1# ' Data 3#, 7#, 10#, 8#, 4# ' Data -7#, 3#, 8#, -0.75, -9# ' Data 12#, -1#, 4#, -9#, 10# ' For I = 1 To N ' For J = 1 To N ' AI(I, J) = 1# ' Next J ' Next I CLS A\$ = " Input Matrix (real part):" GOSUB 500 'print AR ' #2 only: ' A\$=" Input Matrix (imag. part):" ' GOSUB 501 'print AI ' INPUT " ", R\$ IAVEC = 1 'EIGENVECTORS ASKED FOR GOSUB 1000 'CALL CEIGEN(AR,AI,NM,N,TRUE,WR,WI,ZR,ZI,WORK,IERR) PRINT PRINT " Error code = "; IERR PRINT PRINT " Eigenvalues:" FOR I = 1 TO N PRINT " "; WR(I); " "; WI(I) NEXT I PRINT A\$ = " Eigenvectors (real part in lines):" GOSUB 600 'print ZR ' #2 only: ' A\$ = " Eigenvectors (imag. part in lines):" ' GOSUB 601 'print ZI END 'of main program 'print matrix AR with caption 500 PRINT A\$ FOR I = 1 TO N FOR J = 1 TO N PRINT USING F\$; AR(I, J); NEXT J PRINT NEXT I RETURN 'print matrix AI with caption 501 PRINT A\$ FOR I = 1 TO N FOR J = 1 TO N PRINT USING F\$; AI(I, J); NEXT J PRINT NEXT I RETURN 'print matrix ZR with caption 600 PRINT A\$ FOR I = 1 TO N FOR J = 1 TO N PRINT USING F\$; ZR(J, I); NEXT J PRINT NEXT I RETURN 'print matrix ZI with caption 601 PRINT A\$ FOR I = 1 TO N FOR J = 1 TO N PRINT USING F\$; ZI(J, I); NEXT J PRINT NEXT I RETURN 1000 '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. ' -------------------------------------------------------------------- ' INPUTS: ' AR = TABLE(N,N) OF REAL COEF. OF MATRIX A R*8 ' AI = TABLE(N,N) OF IMAGINARY COEF. OF MATRIX A R*8 ' N = DIMENSION OF MATRICES AR,AI,ZR & ZI I*4 ' IAVEC= CODE TO CALCULATE OR NOT EIGENVECTORS L*8 ' = 0 NO EIGENVECTORS ' = 1 WITH EIGENVECTORS ' OUTPUTS: ' WR,WI= TABLEA(N) OF EIGENVALUES R*8 ' ZR,ZI= TABLE(N,N) OF EIGENVECTORS R*8 ' IERR = ERROR CODE I*4 ' = 0 NO ERROR ' < OU = N NO CONVERGENCE ' = 9999 DIMENSION ERROR IN TABLES AR, AI, ZR, ZI ' WORKING SPACE: ' WORK = TABLE(N) R*8 ' REFERENCE: ' J.H.WILKINSON & C.REINSCH ' LINEAR ALGEBRA SPRINGER-VERLAG 1971 ' -------------------------------------------------------------------- IF N <= NM THEN N1 = N + 1 N2 = N + N1 GOSUB 2000 'CALL CBALAN(NM,N,AR,AI,LOW,IUPP,WORK) GOSUB 3000 'CALL CUNITH(NM,N,LOW,IUPP,AR,AI,TMP1,TMP2) GOSUB 4000 'CALL COMQR2(NM,N,LOW,IUPP,TMP,TMP2,AR,AI,WR,WI, 'IAVEC,ZR,ZI,IERR) IF IERR = 0 THEN GOSUB 5000 'CALL CBALBK(NM,N,LOW,IUPP,N,WORK,ZR,ZI) GOSUB 6000 'CALL NORMAL(NM,N,WR,WI,IAVEC,ZR,ZI,TMP1,TMP2) END IF ELSE IERR = 9999 END IF RETURN 2000 'SUBROUTINE CBALAN(NM,N,AR,AI,LOW,IUPP,WORK) ' BALANCE A COMPLEX MATRIX (VERSION DERIVED FROM BALANC) 'Double B2,C,F,G,R,S 'LOGICAL NOCONV 'INTEGER IB IB = 16 B2 = IB * IB L = 1 K = N 2010 FOR J = K TO 1 STEP -1 J1 = J K1 = K L1 = L FOR I = 1 TO K IF I <> J THEN IF AR(J, I) <> 0# OR AI(J, I) <> 0# THEN GOTO 2011 END IF NEXT I M1 = K GOSUB 7000 'CALL EXCHG1(K,NM,N,AR,AI,WORK,J1,K1,L1) K = K - 1 GOTO 2010 2011 NEXT J 2020 FOR J = L TO K J1 = J K1 = K L1 = L FOR I = L TO K IF I <> J THEN IF AR(I, J) <> 0# OR AI(I, J) <> 0# THEN GOTO 2022 END IF NEXT I M1 = L GOSUB 7000 'CALL EXCHG1(L,NM,N,AR,AI,WORK,J1,K1,L1) L = L + 1 GOTO 2020 2022 NEXT J LOW = L IUPP = K FOR I = L TO K WORK(I) = 1# NEXT I 2030 NOCONV = 0 FOR I = L TO K C = 0# R = 0# FOR J = L TO K IF J <> I THEN C = C + ABS(AR(J, I)) + ABS(AI(J, I)) R = R + ABS(AR(I, J)) + ABS(AI(I, J)) END IF NEXT J G = R / IB F = 1# S = C + R 2040 IF C < G THEN F = F * IB C = C * B2 GOTO 2040 END IF G = R * IB 2050 IF C >= G THEN F = F / IB C = C / B2 GOTO 2050 END IF IF (C + R) / F < .95 * S THEN G = 1# / F WORK(I) = WORK(I) * F NOCONV = 1 FOR J = L TO N AR(I, J) = AR(I, J) * G AI(I, J) = AI(I, J) * G NEXT J FOR J = 1 TO K AI(J, I) = AI(J, I) * F AR(J, I) = AR(J, I) * F NEXT J END IF NEXT I IF NOCONV <> 0 THEN GOTO 2030 RETURN 3000 'SUBROUTINE CUNITH(NM,N,LOW,IUPP,AR,AI,TMP1,TMP2) ' UNITARY TRANSFORMATIONS TO REDUCE A COMPLEX MATRIX ' TO A HESSENBERG FORM ' Double XR,XI,G,F,H,S IF IUPP - 1 < LOW + 1 THEN RETURN FOR M = LOW + 1 TO IUPP - 1 TMP1(M) = 0# TMP2(M) = 0# S = 0# H = 0# FOR I = M TO IUPP S = S + ABS(AR(I, M - 1)) + ABS(AI(I, M - 1)) NEXT I IF S <> 0# THEN FOR I = IUPP TO M STEP -1 TMP1(I) = AR(I, M - 1) / S TMP2(I) = AI(I, M - 1) / S H = H + TMP1(I) * TMP1(I) + TMP2(I) * TMP2(I) NEXT I G = SQR(H) AR = TMP1(M): AI = TMP2(M): GOSUB 7500 F = ABSCD IF F = 0# THEN TMP1(M) = G AR(M, M - 1) = S ELSE H = H + F * G G = G / F TMP1(M) = (1# + G) * TMP1(M) TMP2(M) = (1# + G) * TMP2(M) END IF FOR J = M TO N XR = 0# XI = 0# FOR I = IUPP TO M STEP -1 XR = XR + TMP1(I) * AR(I, J) + TMP2(I) * AI(I, J) XI = XI + TMP1(I) * AI(I, J) - TMP2(I) * AR(I, J) NEXT I XR = XR / H XI = XI / H FOR I = M TO IUPP AR(I, J) = AR(I, J) - XR * TMP1(I) + XI * TMP2(I) AI(I, J) = AI(I, J) - XR * TMP2(I) - XI * TMP1(I) NEXT I NEXT J FOR I = 1 TO IUPP XR = 0# XI = 0# FOR J = IUPP TO M STEP -1 XR = XR + TMP1(J) * AR(I, J) - TMP2(J) * AI(I, J) XI = XI + TMP1(J) * AI(I, J) + TMP2(J) * AR(I, J) NEXT J XR = XR / H XI = XI / H FOR J = M TO IUPP AR(I, J) = AR(I, J) - XR * TMP1(J) - XI * TMP2(J) AI(I, J) = AI(I, J) + XR * TMP2(J) - XI * TMP1(J) NEXT J NEXT I TMP1(M) = TMP1(M) * S TMP2(M) = TMP2(M) * S AR(M, M - 1) = -G * AR(M, M - 1) AI(M, M - 1) = -G * AI(M, M - 1) END IF NEXT M RETURN 4000 'SUBROUTINE COMQR2(NM,N,LOW,IUPP,TMP1,TMP2,AR,AI,WR,WI,IAVEC,ZR,ZI,IERR) ' --------------------------------------------------------------------- ' QR ALGORITH APPLIED TO A COMPLEX HESSENBERG MATRIX ' --------------------------------------------------------------------- 'Double SR,SI,TR,TI,YR,YI,XR,XI,VR,VI,XMACHEP,XNORM,EPS,EPS1 'INTEGER I,J,K,L,M,N,IN0,NA,NM,ITS,LOW,IUPP,IERR (Note: IN is a reserved word). 'LOGICAL IAVEC ' SEEK FLOATING POINT RELATIVE PRECISION FOR TYPE DOUBLE EPS = 1# 4005 EPS = EPS * .5# EPS1 = 1# + EPS IF EPS1 > 1# THEN GOTO 4005 XMACHEP = 2# * EPS IERR = 0 IF IAVEC <> 0 THEN FOR I = 1 TO N FOR J = 1 TO N ZR(I, J) = 0# ZI(I, J) = 0# NEXT J ZR(I, I) = 1# NEXT I ' ACCUMULATE UNITARY TRANSFORMATIONS FOR I = IUPP - 1 TO LOW + 1 STEP -1 IF (AR(I, I - 1) = 0# AND AI(I, I - 1) = 0#) OR (TMP1(I) = 0# AND TMP2(I) = 0#) THEN GOTO 4007 XNORM = AR(I, I - 1) * TMP1(I) + AI(I, I - 1) * TMP2(I) FOR K = I + 1 TO IUPP TMP1(K) = AR(K, I - 1) TMP2(K) = AI(K, I - 1) NEXT K FOR J = I TO IUPP SR = 0# SI = 0# FOR K = I TO IUPP SR = SR + TMP1(K) * ZR(K, J) + TMP2(K) * ZI(K, J) SI = SI + TMP1(K) * ZI(K, J) - TMP2(K) * ZR(K, J) NEXT K SR = SR / XNORM SI = SI / XNORM FOR K = I TO IUPP ZR(K, J) = ZR(K, J) + SR * TMP1(K) - SI * TMP2(K) ZI(K, J) = ZI(K, J) + SR * TMP2(K) + SI * TMP1(K) NEXT K NEXT J 4007 NEXT I END IF ' CREATE REAL SUBDIAGONAL L = LOW + 1 FOR I = L TO IUPP IF I + 1 < IUPP THEN M = I + 1 ELSE M = IUPP END IF IF AI(I, I - 1) = 0# THEN GOTO 4009 AR = AR(I, I - 1): AI = AI(I, I - 1): GOSUB 7500 XNORM = ABSCD YR = AR(I, I - 1) / XNORM YI = AI(I, I - 1) / XNORM AR(I, I - 1) = XNORM AI(I, I - 1) = 0# FOR J = I TO N SI = YR * AI(I, J) - YI * AR(I, J) AR(I, J) = YR * AR(I, J) + YI * AI(I, J) AI(I, J) = SI NEXT J FOR J = 1 TO M SI = YR * AI(J, I) + YI * AR(J, I) AR(J, I) = YR * AR(J, I) - YI * AI(J, I) AI(J, I) = SI NEXT J IF IAVEC <> 0 THEN FOR J = LOW TO IUPP SI = YR * ZI(J, I) + YI * ZR(J, I) ZR(J, I) = YR * ZR(J, I) - YI * ZI(J, I) ZI(J, I) = SI NEXT J END IF 4009 NEXT I ' STORE ROOTS ALREADY ISOLATED BY SUBROUTINE BALANC FOR I = 1 TO N IF I < LOW OR I > IUPP THEN WR(I) = AR(I, I) WI(I) = AI(I, I) END IF NEXT I IN0 = IUPP TR = 0# TI = 0# ' SEEK EIGENVALUES 4010 IF IN0 < LOW THEN GOTO 4050 ITS = 0 NA = IN0 - 1 ' WE LOOK AT REAL PART OF SUBDIAGONAL ELEMENT 4015 FOR L = IN0 TO LOW STEP -1 SR = XMACHEP * (ABS(AR(L - 1, L - 1)) + ABS(AI(L - 1, L - 1)) + ABS(AR(L, L)) + ABS(AI(L, L))) IF ABS(AR(L, L - 1)) <= SR OR L = LOW THEN GOTO 4020 NEXT L ' SHIFTING 4020 IF L = IN0 THEN GOTO 4045 IF ITS = 30 THEN GOTO 4060 IF ITS = 10 OR ITS = 20 THEN ' EXTRA SHIFTING SR = ABS(AR(IN0, NA)) + ABS(AR(NA, IN0 - 2)) SI = 0# ELSE SR = AR(IN0, IN0) SI = AI(IN0, IN0) XR = AR(NA, IN0) * AR(IN0, NA) XI = AI(NA, IN0) * AR(IN0, NA) IF XR <> 0# OR XI <> 0# THEN YR = (AR(NA, NA) - SR) * .5 YI = (AI(NA, NA) - SI) * .5 'CALL DCSQRT(YR*YR-YI*YI+XR,2.D0*YR*YI+XI,VR,VI) X = YR * YR - YI * YI + XR: Y = 2# * YR * YI + XI: GOSUB 7600 VR = A: VI = B IF YR * VR + YI * VI < 0# THEN VR = -VR VI = -VI END IF 'CALL DIVCD(XR,XI,YR+VR,YI+VI,VR,VI) AR = XR: AI = XI: BR = YR + VR: BI = YI + VI: GOSUB 7700 VR = ZR: VI = ZI SR = SR - VR SI = SI - VI END IF END IF FOR I = LOW TO IN0 AR(I, I) = AR(I, I) - SR AI(I, I) = AI(I, I) - SI NEXT I TR = TR + SR TI = TI + SI ITS = ITS + 1 ' TRIANGULARY DECOMPOSITION H = Q * R FOR I = L + 1 TO IN0 SR = AR(I, I - 1) VR = AR(I - 1, I - 1) VI = AI(I - 1, I - 1) XNORM = SQR(VR * VR + VI * VI + SR * SR) XR = VR / XNORM XI = VI / XNORM WR(I - 1) = XR WI(I - 1) = XI AR(I - 1, I - 1) = XNORM AI(I - 1, I - 1) = 0# AR(I, I - 1) = 0# AI(I, I - 1) = SR / XNORM FOR J = I TO N YR = AR(I - 1, J) YI = AI(I - 1, J) VR = AR(I, J) VI = AI(I, J) AR(I - 1, J) = XR * YR + XI * YI + AI(I, I - 1) * VR AI(I - 1, J) = XR * YI - XI * YR + AI(I, I - 1) * VI AR(I, J) = XR * VR - XI * VI - AI(I, I - 1) * YR AI(I, J) = XR * VI + XI * VR - AI(I, I - 1) * YI NEXT J NEXT I SI = AI(IN0, IN0) IF SI = 0# THEN GOTO 4035 'XNORM = ABSCD(AR(IN0, IN0), SI) AR = AR(IN0, IN0): AI = SI: GOSUB 7500 XNORM = ABSCD SR = AR(IN0, IN0) / XNORM SI = SI / XNORM AI(IN0, IN0) = 0# AR(IN0, IN0) = XNORM IF IN0 <> N THEN FOR J = IN0 + 1 TO N YR = AR(IN0, J) YI = AI(IN0, J) AR(IN0, J) = SR * YR + SI * YI AI(IN0, J) = SR * YI - SI * YR NEXT J END IF ' INVERSE OPERATION 4035 FOR J = L + 1 TO IN0 XR = WR(J - 1) XI = WI(J - 1) FOR I = 1 TO J YR = AR(I, J - 1) YI = 0# VR = AR(I, J) VI = AI(I, J) IF I <> J THEN YI = AI(I, J - 1) AI(I, J - 1) = XR * YI + XI * YR + AI(J, J - 1) * VI END IF AR(I, J - 1) = XR * YR - XI * YI + AI(J, J - 1) * VR AR(I, J) = XR * VR + XI * VI - AI(J, J - 1) * YR AI(I, J) = XR * VI - XI * VR - AI(J, J - 1) * YI NEXT I IF IAVEC <> 0 THEN FOR I = LOW TO IUPP 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 + AI(J, J - 1) * VR ZI(I, J - 1) = XR * YI + XI * YR + AI(J, J - 1) * VI ZR(I, J) = XR * VR + XI * VI - AI(J, J - 1) * YR ZI(I, J) = XR * VI - XI * VR - AI(J, J - 1) * YI NEXT I END IF NEXT J IF SI = 0# THEN GOTO 4015 FOR I = 1 TO IN0 YR = AR(I, IN0) YI = AI(I, IN0) AR(I, IN0) = SR * YR - SI * YI AI(I, IN0) = SR * YI + SI * YR NEXT I IF IAVEC <> 0 THEN FOR I = LOW TO IUPP YR = ZR(I, IN0) YI = ZI(I, IN0) ZR(I, IN0) = SR * YR - SI * YI ZI(I, IN0) = SR * YI + SI * YR NEXT I END IF GOTO 4015 ' ONE ROOT HAS BEEN FOUND 4045 AR(IN0, IN0) = AR(IN0, IN0) + TR AI(IN0, IN0) = AI(IN0, IN0) + TI WR(IN0) = AR(IN0, IN0) WI(IN0) = AI(IN0, IN0) IN0 = NA GOTO 4010 ' ALL ROOTS HAVE BEEN FOUND ' START CALCULATION OF EIGENVECTORS BY SUBSTITUTION 4050 IF IAVEC <> 0 THEN XNORM = 0# FOR I = 1 TO N FOR J = I TO N XNORM = XNORM + ABS(AR(I, J)) + ABS(AI(I, J)) NEXT J NEXT I IF XNORM = 0# THEN GOTO 4060 FOR IN0 = N TO 2 STEP -1 XR = WR(IN0) XI = WI(IN0) AR(IN0, IN0) = 1# AI(IN0, IN0) = 0# FOR I = IN0 - 1 TO 1 STEP -1 VR = AR(I, IN0) VI = AI(I, IN0) IF I <> IN0 - 1 THEN FOR J = I + 1 TO IN0 - 1 VR = VR + AR(I, J) * AR(J, IN0) - AI(I, J) * AI(J, IN0) VI = VI + AR(I, J) * AI(J, IN0) + AI(I, J) * AR(J, IN0) NEXT J END IF YR = XR - WR(I) YI = XI - WI(I) IF YR = 0# AND YI = 0# THEN YR = XNORM * XMACHEP 'CALL DIVCD(VR,VI,YR,YI,TR,TI) AR = VR: AI = VI: BR = YR: BI = YI: GOSUB 7700 TR = ZR: TI = ZI AR(I, IN0) = TR AI(I, IN0) = TI NEXT I NEXT IN0 FOR I = 1 TO N - 1 IF I < LOW OR I > IUPP THEN FOR J = I + 1 TO N ZR(I, J) = AR(I, J) ZI(I, J) = AI(I, J) NEXT J END IF NEXT I ' END CALCULATION OH HESSENBERG EIGENVECTORS FOR J = N TO LOW + 1 STEP -1 'M=MIN(J-1,IUPP) IF J - 1 < IUPP THEN M = J - 1 ELSE M = IUPP END IF FOR I = LOW TO IUPP VR = ZR(I, J) VI = ZI(I, J) FOR K = LOW TO M VR = VR + ZR(I, K) * AR(K, J) - ZI(I, K) * AI(K, J) VI = VI + ZR(I, K) * AI(K, J) + ZI(I, K) * AR(K, J) NEXT K ZR(I, J) = VR ZI(I, J) = VI NEXT I NEXT J END IF 'IAVEC<>0 ' END CALCULATION OF EIGENVECTORS OF INPUT MATRIX RETURN ' NO CONVERGENCE AFTER 30 ITERATIONS! 4060 IERR = IN0 RETURN END 5000 'SUBROUTINE CBALBK(NM,N,LOW,IUPP,M,WORK ZR ZI) ' This subroutine finds the eigenvectors of a complex matrix, ' as output of subroutine COMQR2. ' Double S M = N FOR I = LOW TO IUPP S = WORK(I) FOR J = 1 TO M ZR(I, J) = ZR(I, J) * S ZI(I, J) = ZI(I, J) * S NEXT J NEXT I FOR L = 1 TO N I = L IF I >= LOW AND I <= IUPP THEN GOTO 5010 IF I < LOW THEN I = LOW - L K = WORK(I) IF K <> I THEN FOR J = 1 TO M S = ZR(I, J) ZR(I, J) = ZR(K, J) ZR(K, J) = S S = ZI(I, J) ZI(I, J) = ZI(K, J) ZI(K, J) = S NEXT J END IF 5010 NEXT L RETURN 6000 'SUBROUTINE NORMAL(NM,N,WR,WI,IAVEC,ZR,ZI,TMP1,TMP2) ' Sort eigenvalues by absolute values in ascending order and ' normalize eigenvectors. ' SORT SOLUTIONS IN ASCENDING ORDER FOR J = 2 TO N VR = WR(J) VI = WI(J) IF IAVEC <> 0 THEN FOR K = 1 TO N TMP1(K) = ZR(K, J) TMP2(K) = ZI(K, J) NEXT K END IF FOR I = J - 1 TO 1 STEP -1 'IF ABSCD(WR(I),WI(I)) <= ABSCD(VR,VI) THEN GOTO 6005 AR = WR(I): AI = WI(I): GOSUB 7500: ABSCD1 = ABSCD AR = VR: AI = VI: GOSUB 7500 IF ABSCD1 <= ABSCD THEN GOTO 6005 WR(I + 1) = WR(I) WI(I + 1) = WI(I) IF IAVEC <> 0 THEN FOR K = 1 TO N ZR(K, I + 1) = ZR(K, I) ZI(K, I + 1) = ZI(K, I) NEXT K END IF NEXT I I = 0 6005 WR(I + 1) = VR WI(I + 1) = VI IF IAVEC <> 0 THEN FOR K = 1 TO N ZR(K, I + 1) = TMP1(K) ZI(K, I + 1) = TMP2(K) NEXT K END IF NEXT J ' NORMALIZE WITH RESPECT TO BIGGEST ELEMENT IF IAVEC <> 0 THEN FOR J = N TO 1 STEP -1 ZM = 0# FOR I = 1 TO N VR = ZR(I, J) VI = ZI(I, J) 'Z1=ABSCD (VR,VI) AR = VR: AI = VI: GOSUB 7500: Z1 = ABSCD IF Z1 >= ZM THEN IM = I ZM = Z1 END IF NEXT I Z1 = ZR(IM, J) Z2 = ZI(IM, J) FOR I = 1 TO N 'CALL DIVCD(ZR(I,J),ZI(I,J),Z1,Z2,ZR,ZI) AR = ZR(I, J): AI = ZI(I, J): BR = Z1: BI = Z2: GOSUB 7700 ZR(I, J) = ZR ZI(I, J) = ZI NEXT I NEXT J END IF RETURN ' Utility subroutines ' Exchange elements of a complex matrix 7000 'SUBROUTINE EXCHG1(M1,NM,N,AR,AI,WORK,J1,K1,L1) 'double F WORK(M1) = J1 IF J1 <> M1 THEN FOR I = 1 TO K1 F = AR(I, J) AR(I, J) = AR(I, M1) AR(I, M1) = F F = AI(I, J) AI(I, J) = AI(I, M1) AI(I, M1) = F NEXT I FOR I = L1 TO N F = AR(J, I) AR(J, I) = AR(M1, I) AR(M1, I) = F F = AI(J, I) AI(J, I) = AI(M1, I) AI(M1, I) = F NEXT I END IF RETURN 7500 'FUNCTION ABSCD(AR,AI) ' ABSOLUTE VALUE OF A COMPLEX NUMBER AR+I*AI ' ABSCD=SQRT(AR**2+AI**2) ' Double XR,XI,W XR = ABS(AR) XI = ABS(AI) IF XR <= XI THEN W = XR XR = XI XI = W END IF IF XI = 0# THEN ABSCD = XR ELSE ABSCD = XR * SQR(1# + (XI / XR) ^ 2) END IF RETURN 7600 'SUBROUTINE DCSQRT(X,Y,A,B) ' SQUARE ROOT OF A COMPLEX NUMBER A+I*B = SQR(X+I*Y) IF X = 0# AND Y = 0# THEN A = 0# B = 0# ELSE AR=X: AI=Y: GOSUB 7500 A = SQR((ABS(X) + ABSCD) * .5) IF X >= 0# THEN B = Y / (A + A) ELSE IF Y < 0# THEN B = -A ELSE B = A A = Y / (B + B) END IF END IF END IF RETURN 7700 '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. 'Double AR,AI,BR,BI,YR,YI,ZR,ZI,W YR = BR YI = BI IF ABS(YR) > 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 END IF RETURN 'end of file tceigen.bas