SUBROUTINE DEIGV(n,A,WR,WI,ZR,ZI) C C****The following subroutines were taken from the NSWCMATH subroutines C****publicly available on the Queen's University mainframe C DOUBLE PRECISION A(n,n),WR(n),WI(n),ZR(n,n),ZI(n,n) ka=n LOW=1 IGH=n CALL DBAL(KA,N,A,LOW,IGH,ZI) CALL DORTH(KA,N,LOW,IGH,A,WR) CALL DORTRN(KA,N,LOW,IGH,A,WR,ZR) CALL DHQR2(KA,N,LOW,IGH,A,WR,WI,ZR,IERR) IF (IERR.NE.0) THEN PRINT *,'Error from the eigenvalue procedure is',IERR RETURN ENDIF CALL DBABK(KA,N,LOW,IGH,ZI,N,ZR) C DO 30 K=1,N IF (WI(K)) 30,10,20 10 DO 11 J=1,N 11 ZI(J,K)=0.D0 GO TO 30 20 KP1=K+1 DO 21 J=1,N ZI(J,K)=ZR(J,KP1) ZR(J,KP1)=ZR(J,K) 21 ZI(J,KP1)=-ZI(J,K) 30 CONTINUE RETURN END C C ------------------------------------------------------------------ C SUBROUTINE DBAL(NM,N,A,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC DOUBLE PRECISION A(NM,N),SCALE(N) DOUBLE PRECISION C,F,G,R,S,B2,RADIX C DOUBLE PRECISION DABS LOGICAL NOCONV C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). C C DBAL BALANCES A DOUBLE PRECISION REAL MATRIX AND ISOLATES C EIGENVALUES WHENEVER POSSIBLE. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C A CONTAINS THE INPUT MATRIX TO BE BALANCED. C C ON OUTPUT- C C A CONTAINS THE BALANCED MATRIX, C C LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) C IS EQUAL TO ZERO IF C (1) I IS GREATER THAN J AND C (2) J=1,...,LOW-1 OR I=IGH+1,...,N, C C SCALE CONTAINS INFORMATION DETERMINING THE C PERMUTATIONS AND SCALING FACTORS USED. C C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN C SCALE(J) = P(J), FOR J = 1,...,LOW-1 C = D(J,J), J = LOW,...,IGH C = P(J) J = IGH+1,...,N. C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, C THEN 1 TO LOW-1. C C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. C C THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN C DBAL IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS C K,L HAVE BEEN REVERSED.) C C ------------------------------------------------------------------ C C ********** RADIX IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE BASE OF THE MACHINE FLOATING POINT REPRESENTATION. C C ********** RADIX = 2 C B2 = RADIX * RADIX K = 1 L = N GO TO 100 C ********** IN-LINE PROCEDURE FOR ROW AND C COLUMN EXCHANGE ********** 20 SCALE(M) = J IF (J .EQ. M) GO TO 50 C DO 30 I = 1, L F = A(I,J) A(I,J) = A(I,M) A(I,M) = F 30 CONTINUE C DO 40 I = K, N F = A(J,I) A(J,I) = A(M,I) A(M,I) = F 40 CONTINUE C 50 GO TO (80,130), IEXC C ********** SEARCH FOR ROWS ISOLATING AN EIGENVALUE C AND PUSH THEM DOWN ********** 80 IF (L .EQ. 1) GO TO 280 L = L - 1 C ********** FOR J=L STEP -1 UNTIL 1 DO -- ********** 100 DO 120 JJ = 1, L J = L + 1 - JJ C DO 110 I = 1, L IF (I .EQ. J) GO TO 110 IF (A(J,I) .NE. 0.D0) GO TO 120 110 CONTINUE C M = L IEXC = 1 GO TO 20 120 CONTINUE C GO TO 140 C ********** SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE C AND PUSH THEM LEFT ********** 130 K = K + 1 C 140 DO 170 J = K, L C DO 150 I = K, L IF (I .EQ. J) GO TO 150 IF (A(I,J) .NE. 0.D0) GO TO 170 150 CONTINUE C M = K IEXC = 2 GO TO 20 170 CONTINUE C ********** NOW BALANCE THE SUBMATRIX IN ROWS K TO L ********** DO 180 I = K, L 180 SCALE(I) = 1.D0 C ********** ITERATIVE LOOP FOR NORM REDUCTION ********** 190 NOCONV = .FALSE. C DO 270 I = K, L C = 0.D0 R = 0.D0 C DO 200 J = K, L IF (J .EQ. I) GO TO 200 C = C + DABS(A(J,I)) R = R + DABS(A(I,J)) 200 CONTINUE C ********** GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ********** IF (C .EQ. 0.D0 .OR. R .EQ. 0.D0) GO TO 270 G = R / RADIX F = 1.D0 S = C + R 210 IF (C .GE. G) GO TO 220 F = F * RADIX C = C * B2 GO TO 210 220 G = R * RADIX 230 IF (C .LT. G) GO TO 240 F = F / RADIX C = C / B2 GO TO 230 C ********** NOW BALANCE ********** 240 IF ((C + R) / F .GE. 0.95D0 * S) GO TO 270 G = 1.D0 / F SCALE(I) = SCALE(I) * F NOCONV = .TRUE. C DO 250 J = K, N 250 A(I,J) = A(I,J) * G C DO 260 J = 1, L 260 A(J,I) = A(J,I) * F C 270 CONTINUE C IF (NOCONV) GO TO 190 C 280 LOW = K IGH = L RETURN C ********** LAST CARD OF DBAL ********** END C ------------------------------------------------------------------ C SUBROUTINE DORTH(NM,N,LOW,IGH,A,ORT) C INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW DOUBLE PRECISION A(NM,N),ORT(IGH) DOUBLE PRECISION F,G,H,SCALE C DOUBLE PRECISION DSQRT,DABS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE DBAL. IF DBAL HAS NOT BEEN USED THEN C SET LOW=1, IGH=N, C C A CONTAINS THE INPUT MATRIX. C C ON OUTPUT- C C A CONTAINS THE HESSENBERG MATRIX. INFORMATION ABOUT C THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION C IS STORED IN THE REMAINING TRIANGLE UNDER THE C HESSENBERG MATRIX, C C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C ------------------------------------------------------------------ C LA = IGH - 1 KP1 = LOW + 1 IF (LA .LT. KP1) GO TO 200 C DO 180 M = KP1, LA H = 0.D0 ORT(M) = 0.D0 SCALE = 0.D0 C ********** SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ********** DO 90 I = M, IGH 90 SCALE = SCALE + DABS(A(I,M-1)) C IF (SCALE .EQ. 0.D0) GO TO 180 MP = M + IGH C ********** FOR I=IGH STEP -1 UNTIL M DO -- ********** DO 100 II = M, IGH I = MP - II ORT(I) = A(I,M-1) / SCALE H = H + ORT(I) * ORT(I) 100 CONTINUE C G = DSQRT(H) IF (ORT(M) .GE. 0.D0) G = -G H = H - ORT(M) * G ORT(M) = ORT(M) - G C ********** FORM (I-(U*UT)/H) * A ********** DO 130 J = M, N F = 0.D0 C ********** FOR I=IGH STEP -1 UNTIL M DO -- ********** DO 110 II = M, IGH I = MP - II F = F + ORT(I) * A(I,J) 110 CONTINUE C F = F / H C DO 120 I = M, IGH 120 A(I,J) = A(I,J) - F * ORT(I) C 130 CONTINUE C ********** FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ********** DO 160 I = 1, IGH F = 0.D0 C ********** FOR J=IGH STEP -1 UNTIL M DO -- ********** DO 140 JJ = M, IGH J = MP - JJ F = F + ORT(J) * A(I,J) 140 CONTINUE C F = F / H C DO 150 J = M, IGH 150 A(I,J) = A(I,J) - F * ORT(J) C 160 CONTINUE C ORT(M) = SCALE * ORT(M) A(M,M-1) = SCALE * G 180 CONTINUE C 200 RETURN C ********** LAST CARD OF DORTH ********** END C C ------------------------------------------------------------------ C SUBROUTINE DORTRN(NM,N,LOW,IGH,A,ORT,Z) C INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1 DOUBLE PRECISION A(NM,IGH),ORT(IGH),Z(NM,N) DOUBLE PRECISION G C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTRANS, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C C THIS SUBROUTINE ACCUMULATES THE ORTHOGONAL SIMILARITY C TRANSFORMATIONS USED IN THE REDUCTION OF A REAL DOUBLE C PRECISION MATRIX TO UPPER HESSENBERG FORM BY DORTH. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE DBAL. IF DBAL HAS NOT BEEN USED THEN C SET LOW=1, IGH=N, C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION BY DORTH C IN ITS STRICT LOWER TRIANGLE, C C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS- C FORMATIONS USED IN THE REDUCTION BY DORTH. C ONLY ELEMENTS LOW THROUGH IGH ARE USED. C C ON OUTPUT- C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY DORTH, C C ORT HAS BEEN ALTERED. C C ------------------------------------------------------------------ C C ********** INITIALIZE Z TO IDENTITY MATRIX ********** DO 80 I = 1, N C DO 60 J = 1, N 60 Z(I,J) = 0.D0 C Z(I,I) = 1.D0 80 CONTINUE C KL = IGH - LOW - 1 IF (KL .LT. 1) GO TO 200 C ********** FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ********** DO 140 MM = 1, KL MP = IGH - MM IF (A(MP,MP-1) .EQ. 0.D0) GO TO 140 MP1 = MP + 1 C DO 100 I = MP1, IGH 100 ORT(I) = A(I,MP-1) C DO 130 J = MP, IGH G = 0.D0 C DO 110 I = MP, IGH 110 G = G + ORT(I) * Z(I,J) C ********** DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES. C DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ********** G = (G / ORT(MP)) / A(MP,MP-1) C DO 120 I = MP, IGH 120 Z(I,J) = Z(I,J) + G * ORT(I) C 130 CONTINUE C 140 CONTINUE C 200 RETURN C ********** LAST CARD OF DORTRN ********** END C C ------------------------------------------------------------------ C SUBROUTINE DHQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR) C INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN, X IGH,ITS,LOW,MP2,ENM2,IERR DOUBLE PRECISION H(NM,N),WR(N),WI(N),Z(NM,N) DOUBLE PRECISION A,B,D,P,Q,R,S,T,U,W,X,Y DOUBLE PRECISION RA,SA,VI,VR,ZZ,NORM,MACHEP C INTEGER MIN0 C DOUBLE PRECISION DSQRT,DABS C DOUBLE PRECISION DPMPAR LOGICAL NOTLAS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A REAL DOUBLE PRECISION UPPER HESSENBERG MATRIX BY THE C QR METHOD. THE EIGENVECTORS OF A REAL GENERAL MATRIX CAN C ALSO BE FOUND IF DORTH AND DORTRN HAVE BEEN USED TO REDUCE C THIS GENERAL MATRIX TO HESSENBERG FORM AND TO ACCUMULATE C THE SIMILARITY TRANSFORMATIONS. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE DBAL. IF DBAL HAS NOT BEEN USED THEN C SET LOW=1, IGH=N, C C H CONTAINS THE UPPER HESSENBERG MATRIX, C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY DORTRN C AFTER THE REDUCTION BY DORTH, IF PERFORMED. IF THE C EIGENVECTORS OF THE HESSENBERG MATRIX ARE DESIRED, C Z MUST CONTAIN THE IDENTITY MATRIX. C C ON OUTPUT- C C H HAS BEEN DESTROYED, C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N, C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. C IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z C CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX C WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH C COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS C EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN C ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 40 ITERATIONS. C C ------------------------------------------------------------------ C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER. ASSIGN C MACHEP THE VALUE U WHERE U IS THE SMALLEST POSITIVE C FLOATING POINT NUMBER SUCH THAT 1.0 + U .GT. 1.0 C IN THE DOUBLE PRECISION ARITHMETIC BEING USED. C C ********** MACHEP = 1.D-12 C IERR = 0 NORM = 0.D0 K = 1 C ********** STORE ROOTS ISOLATED BY DBAL C AND COMPUTE MATRIX NORM ********** DO 50 I = 1, N C DO 40 J = K, N 40 NORM = NORM + DABS(H(I,J)) C K = I IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 WR(I) = H(I,I) WI(I) = 0.D0 50 CONTINUE C EN = IGH T = 0.D0 C ********** SEARCH FOR NEXT EIGENVALUES ********** 60 IF (EN .LT. LOW) GO TO 340 ITS = 0 NA = EN - 1 ENM2 = NA - 1 C ********** LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT C FOR L=EN STEP -1 UNTIL LOW DO -- ********** 70 DO 80 LL = LOW, EN L = EN + LOW - LL IF (L .EQ. LOW) GO TO 100 S = DABS(H(L-1,L-1)) + DABS(H(L,L)) IF (S .EQ. 0.D0) S = NORM IF (DABS(H(L,L-1)) .LE. MACHEP * S) GO TO 100 80 CONTINUE C ********** FORM SHIFT ********** 100 X = H(EN,EN) IF (L .EQ. EN) GO TO 270 Y = H(NA,NA) W = H(EN,NA) * H(NA,EN) IF (L .EQ. NA) GO TO 280 IF (ITS .EQ. 40) GO TO 1000 IF (ITS .NE. 10 .AND. ITS .NE. 20 .AND. ITS .NE. 30) GO TO 130 C ********** FORM EXCEPTIONAL SHIFT ********** T = T + X C DO 120 I = LOW, EN 120 H(I,I) = H(I,I) - X C S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) X = .75D0 * S Y = X W = -.4375D0 * S * S 130 ITS = ITS + 1 C ********** LOOK FOR TWO CONSECUTIVE SMALL C SUB-DIAGONAL ELEMENTS. C FOR M=EN-2 STEP -1 UNTIL L DO -- ********** DO 140 MM = L, ENM2 M = ENM2 + L - MM ZZ = H(M,M) R = X - ZZ S = Y - ZZ P = (R * S - W) / H(M+1,M) + H(M,M+1) Q = H(M+1,M+1) - ZZ - R - S R = H(M+2,M+1) S = DABS(P) + DABS(Q) + DABS(R) P = P / S Q = Q / S R = R / S IF (M .EQ. L) GO TO 150 IF (DABS(H(M,M-1)) * (DABS(Q) + DABS(R)) .LE. MACHEP * DABS(P) X * (DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))) GO TO 150 140 CONTINUE C 150 MP2 = M + 2 C DO 160 I = MP2, EN H(I,I-2) = 0.D0 IF (I .EQ. MP2) GO TO 160 H(I,I-3) = 0.D0 160 CONTINUE C ********** DOUBLE QR STEP INVOLVING ROWS L TO EN AND C COLUMNS M TO EN ********** DO 260 K = M, NA NOTLAS = K .NE. NA IF (K .EQ. M) GO TO 170 P = H(K,K-1) Q = H(K+1,K-1) R = 0.D0 IF (NOTLAS) R = H(K+2,K-1) X = DABS(P) + DABS(Q) + DABS(R) IF (X .EQ. 0.D0) GO TO 260 P = P / X Q = Q / X R = R / X 170 S = DSQRT(P*P + Q*Q + R*R) IF (P .LT. 0.D0) S = -S IF (K .EQ. M) GO TO 180 H(K,K-1) = -S * X GO TO 190 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) 190 P = P + S X = P / S Y = Q / S ZZ = R / S Q = Q / P R = R / P C ********** ROW MODIFICATION ********** DO 210 J = K, N P = H(K,J) + Q * H(K+1,J) IF (.NOT. NOTLAS) GO TO 200 P = P + R * H(K+2,J) H(K+2,J) = H(K+2,J) - P * ZZ 200 H(K+1,J) = H(K+1,J) - P * Y H(K,J) = H(K,J) - P * X 210 CONTINUE C J = MIN0(EN,K+3) C ********** COLUMN MODIFICATION ********** DO 230 I = 1, J P = X * H(I,K) + Y * H(I,K+1) IF (.NOT. NOTLAS) GO TO 220 P = P + ZZ * H(I,K+2) H(I,K+2) = H(I,K+2) - P * R 220 H(I,K+1) = H(I,K+1) - P * Q H(I,K) = H(I,K) - P 230 CONTINUE C ********** ACCUMULATE TRANSFORMATIONS ********** DO 250 I = LOW, IGH P = X * Z(I,K) + Y * Z(I,K+1) IF (.NOT. NOTLAS) GO TO 240 P = P + ZZ * Z(I,K+2) Z(I,K+2) = Z(I,K+2) - P * R 240 Z(I,K+1) = Z(I,K+1) - P * Q Z(I,K) = Z(I,K) - P 250 CONTINUE C 260 CONTINUE C GO TO 70 C ********** ONE ROOT FOUND ********** 270 H(EN,EN) = X + T WR(EN) = H(EN,EN) WI(EN) = 0.D0 EN = NA GO TO 60 C ********** TWO ROOTS FOUND ********** 280 P = (Y - X) / 2.D0 Q = P * P + W ZZ = DSQRT(DABS(Q)) H(EN,EN) = X + T X = H(EN,EN) H(NA,NA) = Y + T IF (Q .LT. 0.D0) GO TO 320 C ********** REAL PAIR ********** IF (P .LT. 0.D0) ZZ = -ZZ ZZ = P + ZZ WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ .NE. 0.D0) WR(EN) = X - W / ZZ WI(NA) = 0.D0 WI(EN) = 0.D0 X = H(EN,NA) S = DABS(X) + DABS(ZZ) P = X / S Q = ZZ / S R = DSQRT(P*P+Q*Q) P = P / R Q = Q / R C ********** ROW MODIFICATION ********** DO 290 J = NA, N ZZ = H(NA,J) H(NA,J) = Q * ZZ + P * H(EN,J) H(EN,J) = Q * H(EN,J) - P * ZZ 290 CONTINUE C ********** COLUMN MODIFICATION ********** DO 300 I = 1, EN ZZ = H(I,NA) H(I,NA) = Q * ZZ + P * H(I,EN) H(I,EN) = Q * H(I,EN) - P * ZZ 300 CONTINUE C ********** ACCUMULATE TRANSFORMATIONS ********** DO 310 I = LOW, IGH ZZ = Z(I,NA) Z(I,NA) = Q * ZZ + P * Z(I,EN) Z(I,EN) = Q * Z(I,EN) - P * ZZ 310 CONTINUE C GO TO 330 C ********** COMPLEX PAIR ********** 320 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ 330 EN = ENM2 GO TO 60 C ********** ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND C VECTORS OF UPPER TRIANGULAR FORM ********** 340 IF (NORM .EQ. 0.D0) GO TO 1001 C ********** FOR EN=N STEP -1 UNTIL 1 DO -- ********** DO 800 NN = 1, N EN = N + 1 - NN P = WR(EN) Q = WI(EN) NA = EN - 1 IF (Q) 710, 600, 800 C ********** REAL VECTOR ********** 600 M = EN H(EN,EN) = 1.D0 IF (NA .EQ. 0) GO TO 800 C ********** FOR I=EN-1 STEP -1 UNTIL 1 DO -- ********** DO 700 II = 1, NA I = EN - II W = H(I,I) - P R = H(I,EN) IF (M .GT. NA) GO TO 620 C DO 610 J = M, NA 610 R = R + H(I,J) * H(J,EN) C 620 IF (WI(I) .GE. 0.D0) GO TO 630 ZZ = W S = R GO TO 700 630 M = I IF (WI(I) .NE. 0.D0) GO TO 640 T = W IF (W .EQ. 0.D0) T = MACHEP * NORM H(I,EN) = -R / T GO TO 700 C ********** SOLVE REAL EQUATIONS ********** 640 X = H(I,I+1) Y = H(I+1,I) Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) T = (X * S - ZZ * R) / Q H(I,EN) = T IF (DABS(X) .LE. DABS(ZZ)) GO TO 650 H(I+1,EN) = (-R - W * T) / X GO TO 700 650 H(I+1,EN) = (-S - Y * T) / ZZ 700 CONTINUE C ********** END REAL VECTOR ********** GO TO 800 C ********** COMPLEX VECTOR ********** 710 M = NA C ********** LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT C EIGENVECTOR MATRIX IS TRIANGULAR ********** IF (DABS(H(EN,NA)) .LE. DABS(H(NA,EN))) GO TO 720 H(NA,NA) = Q / H(EN,NA) H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA) GO TO 730 720 U = H(NA,NA) - P B = -H(NA,EN) / (U * U + Q * Q) H(NA,NA) = B * Q H(NA,EN) = B * U 730 H(EN,NA) = 0.D0 H(EN,EN) = 1.D0 ENM2 = NA - 1 IF (ENM2 .EQ. 0) GO TO 800 C ********** FOR I=EN-2 STEP -1 UNTIL 1 DO -- ********** DO 790 II = 1, ENM2 I = NA - II W = H(I,I) - P RA = 0.D0 SA = H(I,EN) C DO 760 J = M, NA RA = RA + H(I,J) * H(J,NA) SA = SA + H(I,J) * H(J,EN) 760 CONTINUE C IF (WI(I) .GE. 0.D0) GO TO 770 ZZ = W R = RA S = SA GO TO 790 770 M = I IF (WI(I) .NE. 0.D0) GO TO 780 D = W * W + Q * Q H(I,NA) = -(RA * W + SA * Q) / D H(I,EN) = (RA * Q - SA * W) / D GO TO 790 C ********** SOLVE COMPLEX EQUATIONS ********** 780 X = H(I,I+1) Y = H(I+1,I) VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q VI = (WR(I) - P) * 2.D0 * Q IF (VR .EQ. 0.D0 .AND. VI .EQ. 0.D0) VR = MACHEP * NORM X * (DABS(W) + DABS(Q) + DABS(X) + DABS(Y) + DABS(ZZ)) A = X * R - ZZ * RA + Q * SA B = X * S - ZZ * SA - Q * RA D = VR * VR + VI * VI H(I,NA) = (A * VR + B * VI) / D H(I,EN) = (B * VR - A * VI) / D IF (DABS(X) .LE. DABS(ZZ) + DABS(Q)) GO TO 785 H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X GO TO 790 785 A = -R - Y * H(I,NA) B = -S - Y * H(I,EN) D = ZZ * ZZ + Q * Q H(I+1,NA) = (A * ZZ + B * Q) / D H(I+1,EN) = (B * ZZ - A * Q) / D 790 CONTINUE C ********** END COMPLEX VECTOR ********** 800 CONTINUE C ********** END BACK SUBSTITUTION. C VECTORS OF ISOLATED ROOTS ********** DO 840 I = 1, N IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840 C DO 820 J = I, N 820 Z(I,J) = H(I,J) C 840 CONTINUE C ********** MULTIPLY BY TRANSFORMATION MATRIX TO GIVE C VECTORS OF ORIGINAL FULL MATRIX. C FOR J=N STEP -1 UNTIL LOW DO -- ********** DO 880 JJ = LOW, N J = N + LOW - JJ M = MIN0(J,IGH) C DO 880 I = LOW, IGH ZZ = 0.D0 C DO 860 K = LOW, M 860 ZZ = ZZ + Z(I,K) * H(K,J) C Z(I,J) = ZZ 880 CONTINUE C GO TO 1001 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 40 ITERATIONS ********** 1000 IERR = EN 1001 RETURN C ********** LAST CARD OF DHQR2 ********** END c c C C SUBROUTINE DBABK(NM,N,LOW,IGH,SCALE,M,Z) C INTEGER I,J,K,M,N,II,NM,IGH,LOW DOUBLE PRECISION SCALE(N),Z(NM,M) DOUBLE PRECISION S C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK, C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C BALANCED MATRIX DETERMINED BY DBAL. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C LOW AND IGH ARE INTEGERS DETERMINED BY DBAL, C C SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS C AND SCALING FACTORS USED BY DBAL, C C M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED, C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- C VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. C C ON OUTPUT- C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE C TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. C C ------------------------------------------------------------------ C IF (M .EQ. 0) GO TO 200 IF (IGH .EQ. LOW) GO TO 120 C DO 110 I = LOW, IGH S = SCALE(I) C ********** LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED C IF THE FOREGOING STATEMENT IS REPLACED BY C S=1.0/SCALE(I). ********** DO 100 J = 1, M 100 Z(I,J) = Z(I,J) * S C 110 CONTINUE C ********- FOR I=LOW-1 STEP -1 UNTIL 1, C IGH+1 STEP 1 UNTIL N DO -- ********** 120 DO 140 II = 1, N I = II IF (I .GE. LOW .AND. I .LE. IGH) GO TO 140 IF (I .LT. LOW) I = LOW - II K = SCALE(I) IF (K .EQ. I) GO TO 140 C DO 130 J = 1, M S = Z(I,J) Z(I,J) = Z(K,J) Z(K,J) = S 130 CONTINUE C 140 CONTINUE C 200 RETURN C ********** LAST CARD OF DBABK ********** END