C 535_double.for modified rmatin, can generate random A and B SUBROUTINE CQZHES(NM,N,AR,AI,BR,BI,MATZ,ZR,ZI) IMPLICIT NONE INTEGER I,J,K,L,N,K1,LB,L1,NM,NK1,NM1 DOUBLE PRECISION AR(NM,N),AI(NM,N),BR(NM,N),BI(NM,N), $ ZR(NM,N),ZI(NM,N) DOUBLE PRECISION R,S,T,TI,U1,U2,XI,XR,YI,YR,RHO,U1I DOUBLE PRECISION DSQRT,DABS LOGICAL MATZ C THIS SUBROUTINE IS A COMPLEX ANALOGUE OF THE FIRST STEP OF THE C QZ ALGORITHM FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF COMPLEX GENERAL MATRICES AND C REDUCES ONE OF THEM TO UPPER HESSENBERG FORM WITH REAL (AND NON- C NEGATIVE) SUBDIAGONAL ELEMENTS AND THE OTHER TO UPPER TRIANGULAR C FORM USING UNITARY TRANSFORMATIONS. IT IS USUALLY FOLLOWED BY C CQZVAL AND POSSIBLY CQZVEC. 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 MATRICES, C C A=(AR,AI) CONTAINS A COMPLEX GENERAL MATRIX, C C B=(BR,BI) CONTAINS A COMPLEX GENERAL MATRIX, C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE. C C ON OUTPUT- C C A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS C BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO, AND THE C SUBDIAGONAL ELEMENTS HAVE BEEN MADE REAL (AND NON-NEGATIVE), C C B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS C BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO, C C Z=(ZR,ZI) CONTAINS THE PRODUCT OF THE RIGHT HAND C TRANSFORMATIONS IF MATZ HAS BEEN SET TO .TRUE. C OTHERWISE, Z IS NOT REFERENCED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C INITIALIZE Z IF (MATZ) then DO 3 I = 1, N DO 2 J = 1, N ZR(I,J) = 0.0 ZI(I,J) = 0.0 2 CONTINUE ZR(I,I) = 1.0 3 CONTINUE endif C REDUCE B TO UPPER TRIANGULAR FORM WITH C TEMPORARILY REAL DIAGONAL ELEMENTS IF (N .LE. 1) return NM1 = N - 1 DO 100 L = 1, NM1 L1 = L + 1 S = 0.0 DO 20 I = L, N S = S + DABS(BR(I,L)) + DABS(BI(I,L)) 20 CONTINUE IF (S .NE. 0.0) then RHO = 0.0 DO 25 I = L, N BR(I,L) = BR(I,L) / S BI(I,L) = BI(I,L) / S RHO = RHO + BR(I,L)**2 + BI(I,L)**2 25 CONTINUE R = DSQRT(RHO) XR = dsqrt(BR(L,L)**2+BI(L,L)**2) IF (XR .NE. 0.0) then RHO = RHO + XR * R U1 = -BR(L,L) / XR U1I = -BI(L,L) / XR YR = R / XR + 1.0 BR(L,L) = YR * BR(L,L) BI(L,L) = YR * BI(L,L) else BR(L,L) = R U1 = -1.0 U1I = 0.0 endif DO 50 J = L1, N T = 0.0 TI = 0.0 DO 30 I = L, N T = T + BR(I,L) * BR(I,J) + BI(I,L) * BI(I,J) TI = TI + BR(I,L) * BI(I,J) - BI(I,L) * BR(I,J) 30 CONTINUE T = T / RHO TI = TI / RHO DO 40 I = L, N BR(I,J) = BR(I,J) - T * BR(I,L) + TI * BI(I,L) BI(I,J) = BI(I,J) - T * BI(I,L) - TI * BR(I,L) 40 CONTINUE XI = U1 * BI(L,J) - U1I * BR(L,J) BR(L,J) = U1 * BR(L,J) + U1I * BI(L,J) BI(L,J) = XI 50 CONTINUE DO 80 J = 1, N T = 0.0 TI = 0.0 DO 60 I = L, N T = T + BR(I,L) * AR(I,J) + BI(I,L) * AI(I,J) TI = TI + BR(I,L) * AI(I,J) - BI(I,L) * AR(I,J) 60 CONTINUE T = T / RHO TI = TI / RHO DO 70 I = L, N AR(I,J) = AR(I,J) - T * BR(I,L) + TI * BI(I,L) AI(I,J) = AI(I,J) - T * BI(I,L) - TI * BR(I,L) 70 CONTINUE XI = U1 * AI(L,J) - U1I * AR(L,J) AR(L,J) = U1 * AR(L,J) + U1I * AI(L,J) AI(L,J) = XI 80 CONTINUE BR(L,L) = R * S BI(L,L) = 0.0 DO 90 I = L1, N BR(I,L) = 0.0 BI(I,L) = 0.0 90 CONTINUE write(6,*) '100 CONTINUE L=', L endif 100 CONTINUE C REDUCE A TO UPPER HESSENBERG FORM WITH REAL SUBDIAGONAL C ELEMENTS, WHILE KEEPING B TRIANGULAR DO 160 K = 1, NM1 K1 = K + 1 C SET BOTTOM ELEMENT IN K-TH COLUMN OF A REAL IF (AI(N,K) .NE. 0.0) then R = DSQRT(AR(N,K)**2+AI(N,K)**2) U1 = AR(N,K) / R U1I = AI(N,K) / R AR(N,K) = R AI(N,K) = 0.0 DO 103 J = K1, N XI = U1 * AI(N,J) - U1I * AR(N,J) AR(N,J) = U1 * AR(N,J) + U1I * AI(N,J) AI(N,J) = XI 103 CONTINUE XI = U1 * BI(N,N) - U1I * BR(N,N) BR(N,N) = U1 * BR(N,N) + U1I * BI(N,N) BI(N,N) = XI endif IF (K .EQ. NM1) return NK1 = NM1 - K C FOR L=N-1 STEP -1 UNTIL K+1 DO -- DO 150 LB = 1, NK1 L = N - LB L1 = L + 1 C ZERO A(L+1,K) S = DABS(AR(L,K)) + DABS(AI(L,K)) + AR(L1,K) IF (S .NE. 0.0) then U1 = AR(L,K) / S U1I = AI(L,K) / S U2 = AR(L1,K) / S R = DSQRT(U1*U1+U1I*U1I+U2*U2) U1 = U1 / R U1I = U1I / R U2 = U2 / R AR(L,K) = R * S AI(L,K) = 0.0 AR(L1,K) = 0.0 DO 110 J = K1, N XR = AR(L,J) XI = AI(L,J) YR = AR(L1,J) YI = AI(L1,J) AR(L,J) = U1 * XR + U1I * XI + U2 * YR AI(L,J) = U1 * XI - U1I * XR + U2 * YI AR(L1,J) = U1 * YR - U1I * YI - U2 * XR AI(L1,J) = U1 * YI + U1I * YR - U2 * XI 110 CONTINUE XR = BR(L,L) BR(L,L) = U1 * XR BI(L,L) = -U1I * XR BR(L1,L) = -U2 * XR DO 120 J = L1, N XR = BR(L,J) XI = BI(L,J) YR = BR(L1,J) YI = BI(L1,J) BR(L,J) = U1 * XR + U1I * XI + U2 * YR BI(L,J) = U1 * XI - U1I * XR + U2 * YI BR(L1,J) = U1 * YR - U1I * YI - U2 * XR BI(L1,J) = U1 * YI + U1I * YR - U2 * XI 120 CONTINUE C ZERO B(L+1,L) S = DABS(BR(L1,L1)) + DABS(BI(L1,L1)) + DABS(BR(L1,L)) IF (S .NE. 0.0) then U1 = BR(L1,L1) / S U1I = BI(L1,L1) / S U2 = BR(L1,L) / S R = DSQRT(U1*U1+U1I*U1I+U2*U2) U1 = U1 / R U1I = U1I / R U2 = U2 / R BR(L1,L1) = R * S BI(L1,L1) = 0.0 BR(L1,L) = 0.0 DO 130 I = 1, L XR = BR(I,L1) XI = BI(I,L1) YR = BR(I,L) YI = BI(I,L) BR(I,L1) = U1 * XR + U1I * XI + U2 * YR BI(I,L1) = U1 * XI - U1I * XR + U2 * YI BR(I,L) = U1 * YR - U1I * YI - U2 * XR BI(I,L) = U1 * YI + U1I * YR - U2 * XI 130 CONTINUE DO 140 I = 1, N XR = AR(I,L1) XI = AI(I,L1) YR = AR(I,L) YI = AI(I,L) AR(I,L1) = U1 * XR + U1I * XI + U2 * YR AI(I,L1) = U1 * XI - U1I * XR + U2 * YI AR(I,L) = U1 * YR - U1I * YI - U2 * XR AI(I,L) = U1 * YI + U1I * YR - U2 * XI 140 CONTINUE IF (MATZ) then DO 145 I = 1, N XR = ZR(I,L1) XI = ZI(I,L1) YR = ZR(I,L) YI = ZI(I,L) ZR(I,L1) = U1 * XR + U1I * XI + U2 * YR ZI(I,L1) = U1 * XI - U1I * XR + U2 * YI ZR(I,L) = U1 * YR - U1I * YI - U2 * XR ZI(I,L) = U1 * YI + U1I * YR - U2 * XI 145 CONTINUE endif endif endif write(6,*) '150 CONTINUE, LB=', LB 150 continue write(6,*) '160 CONTINUE, K=', K 160 CONTINUE return END SUBROUTINE CQZVAL(NM,N,AR,AI,BR,BI,EPS1,ALFR,ALFI,BETA, X MATZ,ZR,ZI,IERR) IMPLICIT NONE INTEGER I,J,K,L,N,EN,K1,K2,LL,L1,NA,NM,ITS,KM1,LM1, X ENM2,IERR,LOR1,ENORN DOUBLE PRECISION AR(NM,N),AI(NM,N),BR(NM,N),BI(NM,N), 1 ALFR(N),ALFI(N), X BETA(N),ZR(NM,N),ZI(NM,N) DOUBLE PRECISION R,S,A1,A2,EP,SH,U1,U2,XI,XR,YI,YR,ANI,A1I, 1 A33,A34,A43,A44, X BNI,B11,B33,B44,SHI,U1I,A33I,A34I,A43I,A44I,B33I,B44I, X EPSA,EPSB,EPS1,ANORM,BNORM,B3344,B3344I DOUBLE PRECISION DSQRT,CDABS,DABS INTEGER MAX0 LOGICAL MATZ DOUBLE COMPLEX Z3 DOUBLE COMPLEX CDSQRT,DCMPLX DOUBLE PRECISION DREAL,DIMAG C THIS SUBROUTINE IS A COMPLEX ANALOGUE OF STEPS 2 AND 3 OF THE C QZ ALGORITHM FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, C AS MODIFIED IN TECHNICAL NOTE NASA TN E-7305(1973) BY WARD. C C THIS SUBROUTINE ACCEPTS A PAIR OF COMPLEX MATRICES, ONE OF THEM C IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM, C THE HESSENBERG MATRIX MUST FURTHER HAVE REAL SUBDIAGONAL ELEMENTS. C IT REDUCES THE HESSENBERG MATRIX TO TRIANGULAR FORM USING C UNITARY TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM C OF THE OTHER MATRIX AND FURTHER MAKING ITS DIAGONAL ELEMENTS C REAL AND NON-NEGATIVE. IT THEN RETURNS QUANTITIES WHOSE RATIOS C GIVE THE GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY C CQZHES AND POSSIBLY FOLLOWED BY CQZVEC. 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 MATRICES, C C A=(AR,AI) CONTAINS A COMPLEX UPPER HESSENBERG MATRIX C WITH REAL SUBDIAGONAL ELEMENTS, C C B=(BR,BI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX, C C EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. C EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN C ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF C ERROR TIMES THE NORM OF ITS MATRIX. IF THE INPUT EPS1 IS C POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE C IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX. A C POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, C BUT LESS ACCURATE RESULTS, C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE, C C Z=(ZR,ZI) CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE C TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION C BY CQZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. C IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. C C ON OUTPUT- C C A HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS C BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO, C C B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS C HAVE BEEN ALTERED. IN PARTICULAR, ITS DIAGONAL HAS BEEN SET C REAL AND NON-NEGATIVE. THE LOCATION BR(N,1) IS USED TO C STORE EPS1 TIMES THE NORM OF B FOR LATER USE BY CQZVEC, C C ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE C DIAGONAL ELEMENTS OF THE TRIANGULARIZED A MATRIX, C C BETA CONTAINS THE REAL NON-NEGATIVE DIAGONAL ELEMENTS OF THE C CORRESPONDING B. THE GENERALIZED EIGENVALUES ARE THEN C THE RATIOS ((ALFR+I*ALFI)/BETA), C C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS C (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE., C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF AR(J,J-1) HAS NOT BECOME C ZERO AFTER 50 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY IERR = 0 C COMPUTE EPSA,EPSB ANORM = 0.0 BNORM = 0.0 C DO 30 I = 1, N ANI = 0.0 IF (I .NE. 1) ANI = DABS(AR(I,I-1)) BNI = 0.0 C DO 20 J = I, N ANI = ANI + DABS(AR(I,J)) + DABS(AI(I,J)) BNI = BNI + DABS(BR(I,J)) + DABS(BI(I,J)) 20 CONTINUE C IF (ANI .GT. ANORM) ANORM = ANI IF (BNI .GT. BNORM) BNORM = BNI 30 CONTINUE C IF (ANORM .EQ. 0.0) ANORM = 1.0 IF (BNORM .EQ. 0.0) BNORM = 1.0 EP = EPS1 IF (EP .GT. 0.0) GO TO 50 C COMPUTE ROUNDOFF LEVEL IF EPS1 IS ZERO EP = 1.0 40 EP = EP / 2.0 IF (1.0 + EP .GT. 1.0) GO TO 40 50 EPSA = EP * ANORM EPSB = EP * BNORM C REDUCE A TO TRIANGULAR FORM, WHILE C KEEPING B TRIANGULAR LOR1 = 1 ENORN = N EN = N C BEGIN QZ STEP 60 IF (EN .EQ. 0) GO TO 1001 IF (.NOT. MATZ) ENORN = EN ITS = 0 NA = EN - 1 ENM2 = NA - 1 C CHECK FOR CONVERGENCE OR REDUCIBILITY. C FOR L=EN STEP -1 UNTIL 1 DO -- 70 DO 80 LL = 1, EN LM1 = EN - LL L = LM1 + 1 IF (L .EQ. 1) GO TO 95 IF (DABS(AR(L,LM1)) .LE. EPSA) GO TO 90 80 CONTINUE C 90 AR(L,LM1) = 0.0 C SET DIAGONAL ELEMENT AT TOP OF B REAL 95 B11 = DSQRT(BR(L,L)**2+BI(L,L)**2) write(6,*) 'L= ', L IF (B11 .EQ. 0.0) GO TO 98 U1 = BR(L,L) / B11 U1I = BI(L,L) / B11 C DO 97 J = L, ENORN XI = U1 * AI(L,J) - U1I * AR(L,J) AR(L,J) = U1 * AR(L,J) + U1I * AI(L,J) AI(L,J) = XI XI = U1 * BI(L,J) - U1I * BR(L,J) BR(L,J) = U1 * BR(L,J) + U1I * BI(L,J) BI(L,J) = XI 97 CONTINUE C BI(L,L) = 0.0 98 IF (L .NE. EN) GO TO 100 C 1-BY-1 BLOCK ISOLATED ALFR(EN) = AR(EN,EN) ALFI(EN) = AI(EN,EN) BETA(EN) = B11 EN = NA write(6,*)'GO TO 60, EN =',EN GO TO 60 C CHECK FOR SMALL TOP OF B 100 L1 = L + 1 IF (B11 .GT. EPSB) GO TO 120 BR(L,L) = 0.0 S = DABS(AR(L,L)) + DABS(AI(L,L)) + DABS(AR(L1,L)) U1 = AR(L,L) / S U1I = AI(L,L) / S U2 = AR(L1,L) / S R = DSQRT(U1*U1+U1I*U1I+U2*U2) U1 = U1 / R U1I = U1I / R U2 = U2 / R AR(L,L) = R * S AI(L,L) = 0.0 C DO 110 J = L1, ENORN XR = AR(L,J) XI = AI(L,J) YR = AR(L1,J) YI = AI(L1,J) AR(L,J) = U1 * XR + U1I * XI + U2 * YR AI(L,J) = U1 * XI - U1I * XR + U2 * YI AR(L1,J) = U1 * YR - U1I * YI - U2 * XR AI(L1,J) = U1 * YI + U1I * YR - U2 * XI XR = BR(L,J) XI = BI(L,J) YR = BR(L1,J) YI = BI(L1,J) BR(L1,J) = U1 * YR - U1I * YI - U2 * XR BR(L,J) = U1 * XR + U1I * XI + U2 * YR BI(L,J) = U1 * XI - U1I * XR + U2 * YI BI(L1,J) = U1 * YI + U1I * YR - U2 * XI 110 CONTINUE C LM1 = L L = L1 write(6,*) 'GO TO 90, L = ', L GO TO 90 C ITERATION STRATEGY 120 CONTINUE write(6,*) 'ITS', ITS IF (ITS .EQ. 50) GO TO 1000 IF (ITS .EQ. 10) GO TO 135 C DETERMINE SHIFT B33 = BR(NA,NA) B33I = BI(NA,NA) IF (DSQRT(B33**2+B33I**2) .GE. EPSB) GO TO 122 B33 = EPSB B33I = 0.0 122 B44 = BR(EN,EN) B44I = BI(EN,EN) IF (DSQRT(B44**2+B44I**2) .GE. EPSB) GO TO 124 B44 = EPSB B44I = 0.0 124 B3344 = B33 * B44 - B33I * B44I B3344I = B33 * B44I + B33I * B44 A33 = AR(NA,NA) * B44 - AI(NA,NA) * B44I A33I = AR(NA,NA) * B44I + AI(NA,NA) * B44 A34 = AR(NA,EN) * B33 - AI(NA,EN) * B33I X - AR(NA,NA) * BR(NA,EN) + AI(NA,NA) * BI(NA,EN) A34I = AR(NA,EN) * B33I + AI(NA,EN) * B33 X - AR(NA,NA) * BI(NA,EN) - AI(NA,NA) * BR(NA,EN) A43 = AR(EN,NA) * B44 A43I = AR(EN,NA) * B44I A44 = AR(EN,EN) * B33 - AI(EN,EN) * B33I - AR(EN,NA) * BR(NA,EN) A44I = AR(EN,EN) * B33I + AI(EN,EN) * B33 - AR(EN,NA) * BI(NA,EN) SH = A44 SHI = A44I XR = A34 * A43 - A34I * A43I XI = A34 * A43I + A34I * A43 IF (XR .EQ. 0.0 .AND. XI .EQ. 0.0) GO TO 140 YR = (A33 - SH) / 2.0 YI = (A33I - SHI) / 2.0 Z3 = CDSQRT(DCMPLX(YR**2-YI**2+XR,2.0*YR*YI+XI)) U1 = DREAL(Z3) U1I = DIMAG(Z3) IF (YR * U1 + YI * U1I .GE. 0.0) GO TO 125 U1 = -U1 U1I = -U1I 125 Z3 = (DCMPLX(SH,SHI) - DCMPLX(XR,XI) / DCMPLX(YR+U1,YI+U1I)) X / DCMPLX(B3344,B3344I) SH = DREAL(Z3) SHI = DIMAG(Z3) GO TO 140 C AD HOC SHIFT 135 SH = AR(EN,NA) + AR(NA,ENM2) SHI = 0.0 C DETERMINE ZEROTH COLUMN OF A 140 A1 = AR(L,L) / B11 - SH A1I = AI(L,L) / B11 - SHI A2 = AR(L1,L) / B11 ITS = ITS + 1 IF (.NOT. MATZ) LOR1 = L C MAIN LOOP DO 260 K = L, NA K1 = K + 1 K2 = K + 2 KM1 = MAX0(K-1,L) C ZERO A(K+1,K-1) IF (K .EQ. L) GO TO 170 A1 = AR(K,KM1) A1I = AI(K,KM1) A2 = AR(K1,KM1) 170 S = DABS(A1) + DABS(A1I) + DABS(A2) U1 = A1 / S U1I = A1I / S U2 = A2 / S R = DSQRT(U1*U1+U1I*U1I+U2*U2) U1 = U1 / R U1I = U1I / R U2 = U2 / R C DO 180 J = KM1, ENORN XR = AR(K,J) XI = AI(K,J) YR = AR(K1,J) YI = AI(K1,J) AR(K,J) = U1 * XR + U1I * XI + U2 * YR AI(K,J) = U1 * XI - U1I * XR + U2 * YI AR(K1,J) = U1 * YR - U1I * YI - U2 * XR AI(K1,J) = U1 * YI + U1I * YR - U2 * XI XR = BR(K,J) XI = BI(K,J) YR = BR(K1,J) YI = BI(K1,J) BR(K,J) = U1 * XR + U1I * XI + U2 * YR BI(K,J) = U1 * XI - U1I * XR + U2 * YI BR(K1,J) = U1 * YR - U1I * YI - U2 * XR BI(K1,J) = U1 * YI + U1I * YR - U2 * XI 180 CONTINUE C IF (K .EQ. L) GO TO 240 AI(K,KM1) = 0.0 AR(K1,KM1) = 0.0 AI(K1,KM1) = 0.0 C ZERO B(K+1,K) 240 S = DABS(BR(K1,K1)) + DABS(BI(K1,K1)) + DABS(BR(K1,K)) U1 = BR(K1,K1) / S U1I = BI(K1,K1) / S U2 = BR(K1,K) / S R = DSQRT(U1*U1+U1I*U1I+U2*U2) U1 = U1 / R U1I = U1I / R U2 = U2 / R IF (K .EQ. NA) GO TO 245 XR = AR(K2,K1) AR(K2,K1) = U1 * XR AI(K2,K1) = -U1I * XR AR(K2,K) = -U2 * XR C 245 DO 250 I = LOR1, K1 XR = AR(I,K1) XI = AI(I,K1) YR = AR(I,K) YI = AI(I,K) AR(I,K1) = U1 * XR + U1I * XI + U2 * YR AI(I,K1) = U1 * XI - U1I * XR + U2 * YI AR(I,K) = U1 * YR - U1I * YI - U2 * XR AI(I,K) = U1 * YI + U1I * YR - U2 * XI XR = BR(I,K1) XI = BI(I,K1) YR = BR(I,K) YI = BI(I,K) BR(I,K1) = U1 * XR + U1I * XI + U2 * YR BI(I,K1) = U1 * XI - U1I * XR + U2 * YI BR(I,K) = U1 * YR - U1I * YI - U2 * XR BI(I,K) = U1 * YI + U1I * YR - U2 * XI 250 CONTINUE C BI(K1,K1) = 0.0 BR(K1,K) = 0.0 BI(K1,K) = 0.0 IF (.NOT. MATZ) GO TO 260 C DO 255 I = 1, N XR = ZR(I,K1) XI = ZI(I,K1) YR = ZR(I,K) YI = ZI(I,K) ZR(I,K1) = U1 * XR + U1I * XI + U2 * YR ZI(I,K1) = U1 * XI - U1I * XR + U2 * YI ZR(I,K) = U1 * YR - U1I * YI - U2 * XR ZI(I,K) = U1 * YI + U1I * YR - U2 * XI 255 CONTINUE C 260 CONTINUE C SET LAST A SUBDIAGONAL REAL AND END QZ STEP C IF (AI(EN,NA) .EQ. 0.0) GO TO 70 IF (AI(EN,NA) .NE. 0.0) THEN R = DSQRT(AR(EN,NA)**2+AI(EN,NA)**2) U1 = AR(EN,NA) / R U1I = AI(EN,NA) / R AR(EN,NA) = R AI(EN,NA) = 0.0 C DO 270 J = EN, ENORN XI = U1 * AI(EN,J) - U1I * AR(EN,J) AR(EN,J) = U1 * AR(EN,J) + U1I * AI(EN,J) AI(EN,J) = XI XI = U1 * BI(EN,J) - U1I * BR(EN,J) BR(EN,J) = U1 * BR(EN,J) + U1I * BI(EN,J) BI(EN,J) = XI 270 CONTINUE C next line added for "THEN" above ENDIF C write(6,*) 'GO TO 70, ITS,L,EN =', ITS, L, EN GO TO 70 C SET ERROR -- BOTTOM SUBDIAGONAL ELEMENT HAS NOT C BECOME NEGLIGIBLE AFTER 50 ITERATIONS 1000 IERR = EN C SAVE EPSB FOR USE BY CQZVEC 1001 IF (N .GT. 1) BR(N,1) = EPSB RETURN C LAST CARD OF CQZVAL END SUBROUTINE CQZVEC(NM,N,AR,AI,BR,BI,ALFR,ALFI,BETA,ZR,ZI) IMPLICIT NONE INTEGER I,J,K,M,N,EN,II,JJ,NA,NM,NN DOUBLE PRECISION AR(NM,N),AI(NM,N),BR(NM,N),BI(NM,N), 1 ALFR(N),ALFI(N), X BETA(N),ZR(NM,N),ZI(NM,N) DOUBLE PRECISION R,T,RI,TI,XI,ALMI,ALMR,BETM,EPSB DOUBLE PRECISION CDABS DOUBLE COMPLEX Z3 DOUBLE COMPLEX DCMPLX DOUBLE PRECISION DREAL,DIMAG C THIS SUBROUTINE IS A COMPLEX ANALOGUE OF THE FOURTH STEP OF THE C QZ ALGORITHM FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF COMPLEX MATRICES IN UPPER C TRIANGULAR FORM, WHERE ONE OF THEM FURTHER MUST HAVE REAL DIAGONAL C ELEMENTS. IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM C AND TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. C IT IS USUALLY PRECEDED BY CQZHES AND CQZVAL. 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 MATRICES, C C A=(AR,AI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX, C C B=(BR,BI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX WITH REAL C DIAGONAL ELEMENTS. IN ADDITION, LOCATION BR(N,1) CONTAINS C THE TOLERANCE QUANTITY (EPSB) COMPUTED AND SAVED IN CQZVAL, C C ALFR, ALFI, AND BETA ARE VECTORS WITH COMPONENTS WHOSE C RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED C EIGENVALUES. THEY ARE USUALLY OBTAINED FROM CQZVAL, C C Z=(ZR,ZI) CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTIONS BY CQZHES AND CQZVAL, IF PERFORMED. C IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE C DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX. C C ON OUTPUT- C C A IS UNALTERED, C C B HAS BEEN DESTROYED, C C ALFR, ALFI, AND BETA ARE UNALTERED, C C Z CONTAINS THE EIGENVECTORS. EACH EIGENVECTOR IS NORMALIZED C SO THAT THE MODULUS OF ITS LARGEST COMPONENT IS 1.0 . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY IF (N .LE. 1) GO TO 1001 EPSB = BR(N,1) C FOR EN=N STEP -1 UNTIL 2 DO -- DO 800 NN = 2, N EN = N + 2 - NN NA = EN - 1 ALMR = ALFR(EN) ALMI = ALFI(EN) BETM = BETA(EN) C FOR I=EN-1 STEP -1 UNTIL 1 DO -- DO 700 II = 1, NA I = EN - II R = 0.0 RI = 0.0 M = I + 1 C DO 610 J = M, EN T = BETM * AR(I,J) - ALMR * BR(I,J) + ALMI * BI(I,J) TI = BETM * AI(I,J) - ALMR * BI(I,J) - ALMI * BR(I,J) IF (J .NE. EN) then XI = T * BI(J,EN) + TI * BR(J,EN) T = T * BR(J,EN) - TI * BI(J,EN) TI = XI endif R = R + T RI = RI + TI 610 CONTINUE C T = ALMR * BETA(I) - BETM * ALFR(I) TI = ALMI * BETA(I) - BETM * ALFI(I) IF (T .EQ. 0.0 .AND. TI .EQ. 0.0) T = EPSB Z3 = DCMPLX(R,RI) / DCMPLX(T,TI) BR(I,EN) = DREAL(Z3) BI(I,EN) = DIMAG(Z3) 700 CONTINUE C 800 CONTINUE C END BACK SUBSTITUTION. C TRANSFORM TO ORIGINAL COORDINATE SYSTEM. C FOR J=N STEP -1 UNTIL 2 DO -- DO 880 JJ = 2, N J = N + 2 - JJ M = J - 1 C DO 881 I = 1, N C DO 860 K = 1, M ZR(I,J) = ZR(I,J) + ZR(I,K) * BR(K,J) - ZI(I,K) * BI(K,J) ZI(I,J) = ZI(I,J) + ZR(I,K) * BI(K,J) + ZI(I,K) * BR(K,J) 860 CONTINUE C 881 CONTINUE 880 CONTINUE C NORMALIZE SO THAT MODULUS OF LARGEST C COMPONENT OF EACH VECTOR IS 1 DO 950 J = 1, N T = 0.0 C DO 930 I = 1, N R = DSQRT(ZR(I,J)**2+ZI(I,J)**2) IF (R .GT. T) T = R 930 CONTINUE C DO 940 I = 1, N ZR(I,J) = ZR(I,J) / T ZI(I,J) = ZI(I,J) / T 940 CONTINUE C 950 CONTINUE C 1001 RETURN C LAST CARD OF CQZVEC END C THIS DRIVER TESTS QZ FOR THE CLASS OF COMPLEX GENERALIZED MATRIX C SYSTEMS EXHIBITING THE USE OF QZ TO FIND ALL THE EIGENVALUES C AND EIGENVECTORS FOR THE EIGENPROBLEM A*X = (LAMBDA)*B*X . C C THE DIMENSION OF A,AI,B,BI,Z, AND ZI SHOULD BE NM BY NM. C THE DIMENSION OF ALFR,ALFI,BETA, AND NORM SHOULD BE NM. C THE DIMENSION OF AHOLD AND BHOLD SHOULD BE NM BY NM. C THE DIMENSION OF AHOLDI AND BHOLDI SHOULD BE NM BY NM. C HERE NM = 4000. C program five35 parameter (nm = 4000) DOUBLE PRECISION A(nm,nm),B(nm,nm),Z(nm,nm), $ ALFR(nm),ALFI(nm),BETA(nm),NORM(nm), $ RESDUL,EPS1, AI(nm,nm),BI(nm,nm),ZI(nm,nm), $ AHOLD(nm,nm),BHOLD(nm,nm), $ AHOLDI(nm,nm),BHOLDI(nm,nm), NORM1 DOUBLE PRECISION DMAX1, DABS, DSQRT INTEGER ERROR 10 CALL RMATIN(NM,N,A,AI,AHOLD,AHOLDI,0) CALL RMATIN(NM,N,B,BI,BHOLD,BHOLDI,0) WRITE(6,20) N 20 FORMAT(30H1THE FULL MATRIX A OF ORDER, $ I4,22H IS (PRINTED BY ROWS)/) DO 30 I = 1,N DO 31 J = 1,N WRITE(6,*)'A',I,J,A(I,J),AI(I,J) 31 CONTINUE 30 CONTINUE WRITE(6,41) N 41 FORMAT(30H0THE FULL MATRIX B OF ORDER, $ I4,22H IS (PRINTED BY ROWS)/) DO 42 I = 1,N DO 43 J = 1,N WRITE(6,*)'B',I,J,B(I,J),BI(I,J) 43 CONTINUE 42 CONTINUE C EIGENVALUES AND EIGENVECTORS USING CQZVAL AND CQZVEC EPS1 = 0.0 CALL CQZHES(NM,N,A,AI,B,BI,.TRUE.,Z,ZI) WRITE(6,*) 'HESSENBERG RESULTS' DO 44 I=1,N DO 441 J=1,N WRITE(6,*)'A',I,J,A(I,J),AI(I,J) 441 CONTINUE 44 CONTINUE DO 45 I=1,N DO 451 J=1,N WRITE(6,*)'B',I,J,B(I,J),BI(I,J) 451 CONTINUE 45 CONTINUE DO 46 I=1,N DO 461 J=1,N WRITE(6,*)'Z',I,J,Z(I,J),ZI(I,J) 461 CONTINUE 46 CONTINUE CALL CQZVAL(NM,N,A,AI,B,BI,EPS1,ALFR,ALFI,BETA, $ .TRUE.,Z,ZI,ERROR) IF(ERROR .NE. 0.0) WRITE(6,*) 'ERROR = ', ERROR WRITE(6,*) 'CQZVAL RESULTS' DO 47 I=1,N DO 471 J=1,N WRITE(6,*)'A',I,J,A(I,J),AI(I,J) 471 CONTINUE 47 CONTINUE DO 48 I=1,N DO 481 J=1,N WRITE(6,*)'B',I,J,B(I,J),BI(I,J) 481 CONTINUE 48 CONTINUE DO 49 I=1,N DO 491 J=1,N WRITE(6,*)'Z',I,J,Z(I,J),ZI(I,J) 491 CONTINUE 49 CONTINUE CALL CQZVEC(NM,N,A,AI,B,BI,ALFR,ALFI,BETA,Z,ZI) WRITE(6,292) 292 FORMAT(/15X,7HALFR(I),19X,7HALFI(I),19X,7HBETA(I)) DO 295 I = 1,N WRITE(6,293) $ I,ALFR(I),ALFI(I),BETA(I) 293 FORMAT(I2,3(1PE23.6,3X)) 295 CONTINUE DO 296 I=1,N DO 297 J=1,N WRITE(6,*)'Z',I,J,Z(I,J),ZI(I,J) 297 CONTINUE 296 CONTINUE CALL RMATIN(NM,N,A,AI,AHOLD,AHOLDI,1) CALL RMATIN(NM,N,B,BI,BHOLD,BHOLDI,1) CALL RMATIN(NM,N,AHOLD,AHOLDI,Z,ZI,1) CALL CGGWZR(NM,N,A,AI,B,BI,ALFR,ALFI,BETA,AHOLD,AHOLDI, $ NORM,RESDUL) WRITE(6,300) 300 FORMAT(/14X,35HCOMPUTED EIGENVALUE AND EIGENVECTOR,20X,8HRESIDUAL) C NORMALIZE EIGENVECTORS DO 507 I = 1,N NORM1 = 0.0 DO 505 J = 1,N NORM1 = NORM1 + Z(J,I)*Z(J,I) 505 CONTINUE NORM1 = DSQRT(NORM1) DO 506 J = 1,N Z(J,I) = Z(J,I)/NORM1 506 CONTINUE 507 CONTINUE DO 510 K = 1, N BETA(K) = DMAX1(BETA(K),1.0D-20) ALFR(K) = ALFR(K) / BETA(K) ALFI(K) = ALFI(K) / BETA(K) 340 WRITE(6,350) $ K,ALFR(K),ALFI(K),NORM(K),(Z(I,K), $ ZI(I,K),I=1,N) 350 FORMAT(/I2,1P2E23.6,E29.2/(5X,2E23.6)) 510 CONTINUE GO TO 10 END SUBROUTINE RMATIN(NM,N,A,AI,AHOLD,AHOLDI,INITIL) C THIS INPUT SUBROUTINE READS TWO REAL MATRICES A AND AI FROM C SYSIN OF ORDER N. C TO GENERATE THE MATRICES INITIALLY, INITIL IS TO BE 0. C TO REGENERATE THE MATRICES FOR THE PURPOSE OF THE RESIDUAL C CALCULATION, INITIL IS TO BE 1. C C THIS ROUTINE IS CATALOGUED AS EISPDRV4(RSGREADI). IMPLICIT NONE INTEGER N,NM,M,I,J,INITIL DOUBLE PRECISION A(NM,NM),AI(NM,NM),AHOLD(NM,NM),AHOLDI(NM,NM) DOUBLE PRECISION DFLOAT INTEGER IA(nm), IB(nm) IF( INITIL .EQ. 1 ) GO TO 30 READ(5,5) N, M 5 FORMAT(I6,6X,I6) IF(M .EQ. -1) THEN CALL MKRAND(NM,N,A,AI,AHOLD,AHOLDI) return ENDIF IF(M .EQ. -2) THEN CALL MKIDENT(NM,N,A,AI,AHOLD,AHOLDI) return ENDIF IF(M .EQ. -3) THEN WRITE(6,*)'A=1/(i+j) n=',N DO 7 i=1,n DO 8 j=1,n A(i,j) = 1.0/(i+j) AHOLD(i,j) = A(i,j) AI(i,j) = 0.0 AHOLDI(i,j) = AI(i,j) 8 CONTINUE 7 CONTINUE return ENDIF IF( N .EQ. 0 ) GO TO 70 IF (M .NE. 1) GO TO 16 DO 10 I = 1,N READ(5,17) $ (IA(J), J=I,N) DO 9 J = I,N A(I,J) = DFLOAT(IA(J)) A(J,I) = A(I,J) 9 CONTINUE 10 CONTINUE 11 READ(5,5) $ N,M IF( M .NE. 1 ) GO TO 20 DO 15 I = 1,N READ(5,17) $ (IB(J), J=I,N) DO 14 J = I,N AI(I,J)=DFLOAT(IB(J)) AI(J,I)=AI(I,J) 14 CONTINUE 15 CONTINUE GO TO 22 16 CONTINUE DO 18 I = 1,N READ(5,17) $ (IA(J), J=1,N) 17 FORMAT(6I6) DO 19 J = 1,N A(I,J) = DFLOAT(IA(J)) 19 CONTINUE 18 CONTINUE GO TO 11 20 CONTINUE DO 25 I = 1,N READ(5,17) $ (IB(J),J=1,N) DO 251 J = 1,N AI(I,J) = DFLOAT(IB(J)) 251 CONTINUE 25 CONTINUE 22 DO 23 I = 1,N DO 231 J = 1,N AHOLDI(I,J) = AI(I,J) AHOLD(I,J) = A(I,J) 231 CONTINUE 23 CONTINUE RETURN 30 DO 40 I = 1,N DO 401 J = 1,N AI(I,J) = AHOLDI(I,J) A(I,J) = AHOLD(I,J) 401 CONTINUE 40 CONTINUE RETURN 70 WRITE(6,80) 80 FORMAT(47H0END OF DATA FOR SUBROUTINE RMATIN(RSGREADI). /1H1) STOP END double precision function drand(ix) integer a,p,ix,b15,b16,xhi,xalo,leftlo,fhi,k data a/16807/, b15/32768/, b16/65536/, p/2147483647/ xhi=ix/b16 xalo=(ix-xhi*b16)*a leftlo=xalo/b16 fhi=xhi*a+leftlo k=fhi/b15 ix=(((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k if (ix .lt. 0) ix=ix+p drand = dfloat(ix)*4.656612875e-10 return end SUBROUTINE MKRAND(NM,N,A,B,AHOLD,BHOLD) IMPLICIT NONE INTEGER N,NM,I,J DOUBLE PRECISION A(NM,NM),B(NM,NM),AHOLD(NM,NM),BHOLD(NM,NM) DOUBLE PRECISION DRAND INTEGER SEED SAVE SEED DATA SEED / 13 / WRITE(6,*)'mkrand n=',N,', seed=',seed DO 10 I = 1,N DO 11 J = 1,N A(I,J) = drand(SEED) AHOLD(I,J) = A(I,J) B(I,J) = DRAND(SEED) BHOLD(I,J) = B(I,J) 11 CONTINUE 10 CONTINUE WRITE(6,*)'mkrand n=',N,', A(1,1)=',A(1,1) RETURN END SUBROUTINE MKIDENT(NM,N,A,B,AHOLD,BHOLD) IMPLICIT NONE INTEGER N,NM,I,J DOUBLE PRECISION A(NM,NM),B(NM,NM),AHOLD(NM,NM),BHOLD(NM,NM) WRITE(6,*)'MKIDENT n=',N DO 10 I = 1,N DO 11 J = 1,N A(I,J) = 0.0D0 IF(I .EQ. J) A(I,J) = 1.0D0 AHOLD(I,J) = A(I,J) B(I,J) = 0.0D0 BHOLD(I,J) = B(I,J) 11 CONTINUE 10 CONTINUE RETURN END SUBROUTINE CGGWZR(NM,N,A,AI,B,BI,ALFR,ALFI,BETA,Z,ZI,NORM,RESDUL) IMPLICIT NONE INTEGER N,NM,I,K,L DOUBLE PRECISION NORM(N), ALFR(N), ALFI(N), BETA(N), $ CPART(2), A(NM,N), $ B(NM,N), Z(NM,N), XR, XI, S, SUMZ, SUMR, SUMI, $ RESDUL,NORMAB,SUMA,SUMB,NORMA,NORMB, $ SUMR2,SUMI2,SUMR3,SUMI3, AI(NM,N),BI(NM,N),ZI(NM,N) DOUBLE PRECISION CDABS,DABS,DFLOAT,DMAX1 DOUBLE COMPLEX C, C1 DOUBLE COMPLEX DCMPLX EQUIVALENCE(C1,CPART(1)) C THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX C A*Z-B*Z*DIAG(W) WHERE A AND B ARE COMPLEX GENERAL MATRICES, Z IS C A MATRIX WHICH CONTAINS THE EIGENVECTORS OF THE EIGENPROBLEM C A*Z - B*Z*DIAG(W), AND W STANDS FOR A VECTOR OF CORRESPONDING C EIGENVALUES OF THE EIGENPROBLEM OBTAINED FROM THE VECTORS ALFR, C ALFI, AND BETA BY THE CORRESPONDENCES C W(J) = (ALFR(J) + I*ALFI(J)) / BETA(J) . C ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS. C C INPUT- C C NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRICES A AND B, C C A(NM,N),AI(NM,N),B(NM,N), AND BI(NM,N) ARE ARRAYS WHICH C CONTAIN THE MATRICES OF THE SYSTEM, C C Z(NM,N) AND ZI(NM,N) ARE ARRAYS WHICH CONTAIN THE C EIGENVECTORS OF THE SYSTEM, C C ALFR(N), ALFI(N), AND BETA(N) ARE ARRAYS CONTAINING THE C COMPONENTS OF THE EIGENVALUES OF THE SYSTEM. C C OUTPUT- C C Z(NM,N) AND ZI(NM,N) ARE ARRAYS WHICH CONTAIN THE NORMALIZED C APPROXIMATE EIGENVECTORS OF THE SYSTEM. THE EIGENVECTORS C ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST C ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE C EIGENVECTOR DIVIDED BY N IS REAL AND POSITIVE, C C NORM(N) IS AN ARRAY SUCH THAT FOR EACH K, C C ..BETA(K)*A*Z(K)-ALFA*B*Z(K).. C NORM(K) = ----------------------------------------------- C ..Z(K)..*(BETA(K)*..A..+.ALFA(K).*..B..) C C WHERE Z(K) IS THE K-TH EIGENVECTOR AND ALFA = (ALFR + I*ALFI), C C RESDUL IS THE REAL NUMBER C ..BETA*A*Z-ALFA*B*Z../(..Z..*(BETA*..A..+.ALFA.*..B..)) NORMB = 0.0 NORMA = 0.0 RESDUL = 0.0 DO 20 I = 1,N SUMA = 0.0 SUMB = 0.0 DO 10 L = 1,N SUMA = SUMA + DABS(A(L,I)) + DABS(AI(L,I)) SUMB = SUMB + DABS(B(L,I)) + DABS(BI(L,I)) 10 CONTINUE NORMA = DMAX1(NORMA,SUMA) NORMB = DMAX1(NORMB,SUMB) 20 CONTINUE DO 160 I=1,N S = 0.0 SUMZ = 0.0 DO 110 L = 1,N SUMR = 0.0 SUMI = 0.0 SUMR2 = 0.0 SUMI2 = 0.0 SUMZ = SUMZ + DSQRT(Z(L,I)**2+ZI(L,I)**2) DO 100 K=1,N SUMR = SUMR + B(L,K)*Z(K,I) - BI(L,K)*ZI(K,I) SUMI = SUMI + B(L,K)*ZI(K,I) + BI(L,K)*Z(K,I) SUMR2 = SUMR2 + A(L,K)*Z(K,I) - AI(L,K)*ZI(K,I) SUMI2 = SUMI2 + A(L,K)*ZI(K,I) + AI(L,K)*Z(K,I) 100 continue SUMR3 = -ALFR(I)*SUMR + ALFI(I)*SUMI SUMI3 = -ALFI(I)*SUMR - ALFR(I)*SUMI S = S+CDABS(DCMPLX(SUMR3,SUMI3)+DCMPLX(SUMR2,SUMI2)*BETA(I)) 110 continue NORMAB = NORMA*BETA(I)+NORMB*DSQRT(ALFR(I)**2+ALFI(I)**2) IF( NORMAB .EQ. 0.0 ) NORMAB = 1.0 NORM(I) = SUMZ IF( SUMZ .EQ. 0.0 ) GO TO 150 C THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL C ALWAYS EXIST AN ELEMENT IN THE VECTOR (Z(I),ZI(I)) C LARGER THAN ..(Z(I),ZI(I)../N DO 120 L=1,N IF(DSQRT(Z(L,I)**2+ZI(L,I)**2) .GE. NORM(I)/DFLOAT(N)) $ GO TO 130 120 CONTINUE 130 XR = NORM(I)*Z(L,I)/DSQRT(Z(L,I)**2+ZI(L,I)**2) XI = NORM(I)*ZI(L,I)/DSQRT(Z(L,I)**2+ZI(L,I)**2) C = DCMPLX(XR,XI) DO 140 L= 1,N C1 = DCMPLX(Z(L,I),ZI(L,I))/C Z(L,I) = CPART(1) ZI(L,I) = CPART(2) 140 continue NORM(I) = S/(NORM(I)*NORMAB) 150 RESDUL = DMAX1(NORM(I),RESDUL) 160 CONTINUE RETURN END