! 535_double.f90 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 ! THIS SUBROUTINE IS A COMPLEX ANALOGUE OF THE FIRST STEP OF THE ! QZ ALGORITHM FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, ! SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. ! ! THIS SUBROUTINE ACCEPTS A PAIR OF COMPLEX GENERAL MATRICES AND ! REDUCES ONE OF THEM TO UPPER HESSENBERG FORM WITH REAL (AND NON- ! NEGATIVE) SUBDIAGONAL ELEMENTS AND THE OTHER TO UPPER TRIANGULAR ! FORM USING UNITARY TRANSFORMATIONS. IT IS USUALLY FOLLOWED BY ! CQZVAL AND POSSIBLY CQZVEC. ! ! ON INPUT- ! ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ! ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM ! DIMENSION STATEMENT, ! ! N IS THE ORDER OF THE MATRICES, ! ! A=(AR,AI) CONTAINS A COMPLEX GENERAL MATRIX, ! ! B=(BR,BI) CONTAINS A COMPLEX GENERAL MATRIX, ! ! MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS ! ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING ! EIGENVECTORS, AND TO .FALSE. OTHERWISE. ! ! ON OUTPUT- ! ! A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS ! BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO, AND THE ! SUBDIAGONAL ELEMENTS HAVE BEEN MADE REAL (AND NON-NEGATIVE), ! ! B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS ! BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO, ! ! Z=(ZR,ZI) CONTAINS THE PRODUCT OF THE RIGHT HAND ! TRANSFORMATIONS IF MATZ HAS BEEN SET TO .TRUE. ! OTHERWISE, Z IS NOT REFERENCED. ! ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY ! 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 ! REDUCE B TO UPPER TRIANGULAR FORM WITH ! 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 ! REDUCE A TO UPPER HESSENBERG FORM WITH REAL SUBDIAGONAL ! ELEMENTS, WHILE KEEPING B TRIANGULAR DO 160 K = 1, NM1 K1 = K + 1 ! 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 ! FOR L=N-1 STEP -1 UNTIL K+1 DO -- DO 150 LB = 1, NK1 L = N - LB L1 = L + 1 ! 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 ! 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, & MATZ,ZR,ZI,IERR) IMPLICIT NONE INTEGER I,J,K,L,N,EN,K1,K2,LL,L1,NA,NM,ITS,KM1,LM1, & ENM2,IERR,LOR1,ENORN DOUBLE PRECISION AR(NM,N),AI(NM,N),BR(NM,N),BI(NM,N), & ALFR(N),ALFI(N),BETA(N),ZR(NM,N),ZI(NM,N) DOUBLE PRECISION R,S,A1,A2,EP,SH,U1,U2,XI,XR,YI,YR,ANI,A1I,& A33,A34,A43,A44,BNI,B11,B33,B44,SHI,U1I,& A33I,A34I,A43I,A44I,B33I,B44I,& 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 ! THIS SUBROUTINE IS A COMPLEX ANALOGUE OF STEPS 2 AND 3 OF THE ! QZ ALGORITHM FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, ! SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, ! AS MODIFIED IN TECHNICAL NOTE NASA TN E-7305(1973) BY WARD. ! ! THIS SUBROUTINE ACCEPTS A PAIR OF COMPLEX MATRICES, ONE OF THEM ! IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM, ! THE HESSENBERG MATRIX MUST FURTHER HAVE REAL SUBDIAGONAL ELEMENTS. ! IT REDUCES THE HESSENBERG MATRIX TO TRIANGULAR FORM USING ! UNITARY TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM ! OF THE OTHER MATRIX AND FURTHER MAKING ITS DIAGONAL ELEMENTS ! REAL AND NON-NEGATIVE. IT THEN RETURNS QUANTITIES WHOSE RATIOS ! GIVE THE GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY ! CQZHES AND POSSIBLY FOLLOWED BY CQZVEC. ! ! ON INPUT- ! ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ! ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM ! DIMENSION STATEMENT, ! ! N IS THE ORDER OF THE MATRICES, ! ! A=(AR,AI) CONTAINS A COMPLEX UPPER HESSENBERG MATRIX ! WITH REAL SUBDIAGONAL ELEMENTS, ! ! B=(BR,BI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX, ! ! EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. ! EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN ! ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF ! ERROR TIMES THE NORM OF ITS MATRIX. IF THE INPUT EPS1 IS ! POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE ! IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX. A ! POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, ! BUT LESS ACCURATE RESULTS, ! ! MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS ! ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING ! EIGENVECTORS, AND TO .FALSE. OTHERWISE, ! ! Z=(ZR,ZI) CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE ! TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION ! BY CQZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. ! IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. ! ! ON OUTPUT- ! ! A HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS ! BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO, ! ! B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS ! HAVE BEEN ALTERED. IN PARTICULAR, ITS DIAGONAL HAS BEEN SET ! REAL AND NON-NEGATIVE. THE LOCATION BR(N,1) IS USED TO ! STORE EPS1 TIMES THE NORM OF B FOR LATER USE BY CQZVEC, ! ! ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE ! DIAGONAL ELEMENTS OF THE TRIANGULARIZED A MATRIX, ! ! BETA CONTAINS THE REAL NON-NEGATIVE DIAGONAL ELEMENTS OF THE ! CORRESPONDING B. THE GENERALIZED EIGENVALUES ARE THEN ! THE RATIOS ((ALFR+I*ALFI)/BETA), ! ! Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS ! (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE., ! ! IERR IS SET TO ! ZERO FOR NORMAL RETURN, ! J IF AR(J,J-1) HAS NOT BECOME ! ZERO AFTER 50 ITERATIONS. ! ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY IERR = 0 ! COMPUTE EPSA,EPSB ANORM = 0.0 BNORM = 0.0 ! DO 30 I = 1, N ANI = 0.0 IF (I .NE. 1) ANI = DABS(AR(I,I-1)) BNI = 0.0 ! 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 ! IF (ANI .GT. ANORM) ANORM = ANI IF (BNI .GT. BNORM) BNORM = BNI 30 CONTINUE ! 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 ! 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 ! REDUCE A TO TRIANGULAR FORM, WHILE ! KEEPING B TRIANGULAR LOR1 = 1 ENORN = N EN = N ! BEGIN QZ STEP 60 IF (EN .EQ. 0) GO TO 1001 IF (.NOT. MATZ) ENORN = EN ITS = 0 NA = EN - 1 ENM2 = NA - 1 ! CHECK FOR CONVERGENCE OR REDUCIBILITY. ! 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 ! 90 AR(L,LM1) = 0.0 ! 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 ! 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 ! BI(L,L) = 0.0 98 IF (L .NE. EN) GO TO 100 ! 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 ! 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 ! 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 ! LM1 = L L = L1 write(6,*) 'GO TO 90, L = ', L GO TO 90 ! ITERATION STRATEGY 120 CONTINUE write(6,*) 'ITS', ITS IF (ITS .EQ. 50) GO TO 1000 IF (ITS .EQ. 10) GO TO 135 ! 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 & - AR(NA,NA) * BR(NA,EN) + AI(NA,NA) * BI(NA,EN) A34I = AR(NA,EN) * B33I + AI(NA,EN) * B33 & - 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)) & / DCMPLX(B3344,B3344I) SH = DREAL(Z3) SHI = DIMAG(Z3) GO TO 140 ! AD HOC SHIFT 135 SH = AR(EN,NA) + AR(NA,ENM2) SHI = 0.0 ! 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 ! MAIN LOOP DO 260 K = L, NA K1 = K + 1 K2 = K + 2 KM1 = MAX0(K-1,L) ! 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 ! 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 ! IF (K .EQ. L) GO TO 240 AI(K,KM1) = 0.0 AR(K1,KM1) = 0.0 AI(K1,KM1) = 0.0 ! 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 ! 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 ! BI(K1,K1) = 0.0 BR(K1,K) = 0.0 BI(K1,K) = 0.0 IF (.NOT. MATZ) GO TO 260 ! 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 ! 260 CONTINUE ! SET LAST A SUBDIAGONAL REAL AND END QZ STEP ! 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 ! 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 ! next line added for "THEN" above ENDIF ! write(6,*) 'GO TO 70, ITS,L,EN =', ITS, L, EN GO TO 70 ! SET ERROR -- BOTTOM SUBDIAGONAL ELEMENT HAS NOT ! BECOME NEGLIGIBLE AFTER 50 ITERATIONS 1000 IERR = EN ! SAVE EPSB FOR USE BY CQZVEC 1001 IF (N .GT. 1) BR(N,1) = EPSB RETURN ! 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), & ALFR(N),ALFI(N),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 ! THIS SUBROUTINE IS A COMPLEX ANALOGUE OF THE FOURTH STEP OF THE ! QZ ALGORITHM FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, ! SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. ! ! THIS SUBROUTINE ACCEPTS A PAIR OF COMPLEX MATRICES IN UPPER ! TRIANGULAR FORM, WHERE ONE OF THEM FURTHER MUST HAVE REAL DIAGONAL ! ELEMENTS. IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ! AND TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. ! IT IS USUALLY PRECEDED BY CQZHES AND CQZVAL. ! ! ON INPUT- ! ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ! ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM ! DIMENSION STATEMENT, ! ! N IS THE ORDER OF THE MATRICES, ! ! A=(AR,AI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX, ! ! B=(BR,BI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX WITH REAL ! DIAGONAL ELEMENTS. IN ADDITION, LOCATION BR(N,1) CONTAINS ! THE TOLERANCE QUANTITY (EPSB) COMPUTED AND SAVED IN CQZVAL, ! ! ALFR, ALFI, AND BETA ARE VECTORS WITH COMPONENTS WHOSE ! RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED ! EIGENVALUES. THEY ARE USUALLY OBTAINED FROM CQZVAL, ! ! Z=(ZR,ZI) CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE ! REDUCTIONS BY CQZHES AND CQZVAL, IF PERFORMED. ! IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE ! DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX. ! ! ON OUTPUT- ! ! A IS UNALTERED, ! ! B HAS BEEN DESTROYED, ! ! ALFR, ALFI, AND BETA ARE UNALTERED, ! ! Z CONTAINS THE EIGENVECTORS. EACH EIGENVECTOR IS NORMALIZED ! SO THAT THE MODULUS OF ITS LARGEST COMPONENT IS 1.0 . ! ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY IF (N .LE. 1) GO TO 1001 EPSB = BR(N,1) ! 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) ! 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 ! 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 ! 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 ! 800 CONTINUE ! END BACK SUBSTITUTION. ! TRANSFORM TO ORIGINAL COORDINATE SYSTEM. ! FOR J=N STEP -1 UNTIL 2 DO -- DO 880 JJ = 2, N J = N + 2 - JJ M = J - 1 ! DO 880 I = 1, N ! 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 ! 880 CONTINUE ! NORMALIZE SO THAT MODULUS OF LARGEST ! COMPONENT OF EACH VECTOR IS 1 DO 950 J = 1, N T = 0.0 ! DO 930 I = 1, N R = DSQRT(ZR(I,J)**2+ZI(I,J)**2) IF (R .GT. T) T = R 930 CONTINUE ! DO 940 I = 1, N ZR(I,J) = ZR(I,J) / T ZI(I,J) = ZI(I,J) / T 940 CONTINUE ! 950 CONTINUE ! 1001 RETURN ! LAST CARD OF CQZVEC END ! THIS DRIVER TESTS QZ FOR THE CLASS OF COMPLEX GENERALIZED MATRIX ! SYSTEMS EXHIBITING THE USE OF QZ TO FIND ALL THE EIGENVALUES ! AND EIGENVECTORS FOR THE EIGENPROBLEM A*X = (LAMBDA)*B*X . ! ! THE DIMENSION OF A,AI,B,BI,Z, AND ZI SHOULD BE NM BY NM. ! THE DIMENSION OF ALFR,ALFI,BETA, AND NORM SHOULD BE NM. ! THE DIMENSION OF AHOLD AND BHOLD SHOULD BE NM BY NM. ! THE DIMENSION OF AHOLDI AND BHOLDI SHOULD BE NM BY NM. ! HERE NM = 4000. ! 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) DOUBLE PRECISION DMAX1 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 30 J = 1,N WRITE(6,*)'A',I,J,A(I,J),AI(I,J) 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 42 J = 1,N WRITE(6,*)'B',I,J,B(I,J),BI(I,J) 42 CONTINUE ! 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 44 J=1,N WRITE(6,*)'A',I,J,A(I,J),AI(I,J) 44 CONTINUE DO 45 I=1,N DO 45 J=1,N WRITE(6,*)'B',I,J,B(I,J),BI(I,J) 45 CONTINUE DO 46 I=1,N DO 46 J=1,N WRITE(6,*)'Z',I,J,Z(I,J),ZI(I,J) 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 47 J=1,N WRITE(6,*)'A',I,J,A(I,J),AI(I,J) 47 CONTINUE DO 48 I=1,N DO 48 J=1,N WRITE(6,*)'B',I,J,B(I,J),BI(I,J) 48 CONTINUE DO 49 I=1,N DO 49 J=1,N WRITE(6,*)'Z',I,J,Z(I,J),ZI(I,J) 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 296 J=1,N WRITE(6,*)'Z',I,J,Z(I,J),ZI(I,J) 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) 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) ! THIS INPUT SUBROUTINE READS TWO REAL MATRICES A AND AI FROM ! SYSIN OF ORDER N. ! TO GENERATE THE MATRICES INITIALLY, INITIL IS TO BE 0. ! TO REGENERATE THE MATRICES FOR THE PURPOSE OF THE RESIDUAL ! CALCULATION, INITIL IS TO BE 1. ! ! 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( 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)) 9 A(J,I) = A(I,J) 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)) 14 AI(J,I)=AI(I,J) 15 CONTINUE GO TO 22 16 DO 18 I = 1,N READ(5,17) (IA(J), J=1,N) 17 FORMAT(6I6) DO 18 J = 1,N 18 A(I,J) = DFLOAT(IA(J)) GO TO 11 20 DO 25 I = 1,N READ(5,17) (IB(J),J=1,N) DO 25 J = 1,N 25 AI(I,J) = DFLOAT(IB(J)) 22 DO 23 I = 1,N DO 23 J = 1,N AHOLDI(I,J) = AI(I,J) 23 AHOLD(I,J) = A(I,J) RETURN 30 DO 40 I = 1,N DO 40 J = 1,N AI(I,J) = AHOLDI(I,J) 40 A(I,J) = AHOLD(I,J) 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 10 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) 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 10 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) 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)) ! THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX ! A*Z-B*Z*DIAG(W) WHERE A AND B ARE COMPLEX GENERAL MATRICES, Z IS ! A MATRIX WHICH CONTAINS THE EIGENVECTORS OF THE EIGENPROBLEM ! A*Z - B*Z*DIAG(W), AND W STANDS FOR A VECTOR OF CORRESPONDING ! EIGENVALUES OF THE EIGENPROBLEM OBTAINED FROM THE VECTORS ALFR, ! ALFI, AND BETA BY THE CORRESPONDENCES ! W(J) = (ALFR(J) + I*ALFI(J)) / BETA(J) . ! ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS. ! ! INPUT- ! ! NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS ! AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT, ! ! N IS THE ORDER OF THE MATRICES A AND B, ! ! A(NM,N),AI(NM,N),B(NM,N), AND BI(NM,N) ARE ARRAYS WHICH ! CONTAIN THE MATRICES OF THE SYSTEM, ! ! Z(NM,N) AND ZI(NM,N) ARE ARRAYS WHICH CONTAIN THE ! EIGENVECTORS OF THE SYSTEM, ! ! ALFR(N), ALFI(N), AND BETA(N) ARE ARRAYS CONTAINING THE ! COMPONENTS OF THE EIGENVALUES OF THE SYSTEM. ! ! OUTPUT- ! ! Z(NM,N) AND ZI(NM,N) ARE ARRAYS WHICH CONTAIN THE NORMALIZED ! APPROXIMATE EIGENVECTORS OF THE SYSTEM. THE EIGENVECTORS ! ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST ! ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE ! EIGENVECTOR DIVIDED BY N IS REAL AND POSITIVE, ! ! NORM(N) IS AN ARRAY SUCH THAT FOR EACH K, ! ! ..BETA(K)*A*Z(K)-ALFA*B*Z(K).. ! NORM(K) = ----------------------------------------------- ! ..Z(K)..*(BETA(K)*..A..+.ALFA(K).*..B..) ! ! WHERE Z(K) IS THE K-TH EIGENVECTOR AND ALFA = (ALFR + I*ALFI), ! ! RESDUL IS THE REAL NUMBER ! ..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)) 10 SUMB = SUMB + DABS(B(L,I)) + DABS(BI(L,I)) NORMA = DMAX1(NORMA,SUMA) 20 NORMB = DMAX1(NORMB,SUMB) 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 ! THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL ! ALWAYS EXIST AN ELEMENT IN THE VECTOR (Z(I),ZI(I)) ! 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