C C ------------------------------------------------------------------ C SUBROUTINE CQZHES(NM,N,AR,AI,BR,BI,MATZ,ZR,ZI) C INTEGER I,J,K,L,N,K1,LB,L1,NM,NK1,NM1 REAL AR(NM,N),AI(NM,N),BR(NM,N),BI(NM,N),ZR(NM,N),ZI(NM,N) REAL R,S,T,TI,U1,U2,XI,XR,YI,YR,RHO,U1I REAL SQRT,CABS,ABS LOGICAL MATZ COMPLEX CMPLX C 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 C ------------------------------------------------------------------ C C ********** INITIALIZE Z ********** IF (.NOT. MATZ) GO TO 10 C DO 3 I = 1, N C DO 2 J = 1, N ZR(I,J) = 0.0 ZI(I,J) = 0.0 2 CONTINUE C ZR(I,I) = 1.0 3 CONTINUE C ********** REDUCE B TO UPPER TRIANGULAR FORM WITH C TEMPORARILY REAL DIAGONAL ELEMENTS ********** 10 IF (N .LE. 1) GO TO 170 NM1 = N - 1 C DO 100 L = 1, NM1 L1 = L + 1 S = 0.0 C DO 20 I = L, N S = S + ABS(BR(I,L)) + ABS(BI(I,L)) 20 CONTINUE C IF (S .EQ. 0.0) GO TO 100 RHO = 0.0 C 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 C R = SQRT(RHO) XR = CABS(CMPLX(BR(L,L),BI(L,L))) IF (XR .EQ. 0.0) GO TO 27 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) GO TO 28 C 27 BR(L,L) = R U1 = -1.0 U1I = 0.0 C 28 DO 50 J = L1, N T = 0.0 TI = 0.0 C 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 C T = T / RHO TI = TI / RHO C 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 C 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 C DO 80 J = 1, N T = 0.0 TI = 0.0 C 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 C T = T / RHO TI = TI / RHO C 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 C 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 C BR(L,L) = R * S BI(L,L) = 0.0 C DO 90 I = L1, N BR(I,L) = 0.0 BI(I,L) = 0.0 90 CONTINUE C 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 ********** R = CABS(CMPLX(AR(N,K),AI(N,K))) IF (R .EQ. 0.0) GO TO 105 U1 = AR(N,K) / R U1I = AI(N,K) / R AR(N,K) = R AI(N,K) = 0.0 C 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 C XI = U1 * BI(N,N) - U1I * BR(N,N) BR(N,N) = U1 * BR(N,N) + U1I * BI(N,N) BI(N,N) = XI 105 IF (K .EQ. NM1) GO TO 170 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 = ABS(AR(L,K)) + ABS(AI(L,K)) + AR(L1,K) IF (S .EQ. 0.0) GO TO 150 U1 = AR(L,K) / S U1I = AI(L,K) / S U2 = AR(L1,K) / S R = SQRT(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 C 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 C XR = BR(L,L) BR(L,L) = U1 * XR BI(L,L) = -U1I * XR BR(L1,L) = -U2 * XR C 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 = ABS(BR(L1,L1)) + ABS(BI(L1,L1)) + ABS(BR(L1,L)) IF (S .EQ. 0.0) GO TO 150 U1 = BR(L1,L1) / S U1I = BI(L1,L1) / S U2 = BR(L1,L) / S R = SQRT(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 C 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 C 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 C IF (.NOT. MATZ) GO TO 150 C 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 C 150 CONTINUE C 160 CONTINUE C 170 RETURN C ********** LAST CARD OF CQZHES ********** END C C ------------------------------------------------------------------ C SUBROUTINE CQZVAL(NM,N,AR,AI,BR,BI,EPS1,ALFR,ALFI,BETA, X MATZ,ZR,ZI,IERR) C INTEGER I,J,K,L,N,EN,K1,K2,LL,L1,NA,NM,ITS,KM1,LM1, X ENM2,IERR,LOR1,ENORN REAL AR(NM,N),AI(NM,N),BR(NM,N),BI(NM,N),ALFR(N),ALFI(N), X BETA(N),ZR(NM,N),ZI(NM,N) REAL R,S,A1,A2,EP,SH,U1,U2,XI,XR,YI,YR,ANI,A1I,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 REAL SQRT,CABS,ABS INTEGER MAX0 LOGICAL MATZ COMPLEX Z3 COMPLEX CSQRT,CMPLX REAL REAL,AIMAG C C C C C 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 C C ------------------------------------------------------------------ C 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 = ABS(AR(I,I-1)) BNI = 0.0 C DO 20 J = I, N ANI = ANI + ABS(AR(I,J)) + ABS(AI(I,J)) BNI = BNI + ABS(BR(I,J)) + ABS(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 (ABS(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 = CABS(CMPLX(BR(L,L),BI(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 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 = ABS(AR(L,L)) + ABS(AI(L,L)) + ABS(AR(L1,L)) U1 = AR(L,L) / S U1I = AI(L,L) / S U2 = AR(L1,L) / S R = SQRT(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 GO TO 90 C ********** ITERATION STRATEGY ********** 120 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 (CABS(CMPLX(B33,B33I)) .GE. EPSB) GO TO 122 B33 = EPSB B33I = 0.0 122 B44 = BR(EN,EN) B44I = BI(EN,EN) IF (CABS(CMPLX(B44,B44I)) .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 = CSQRT(CMPLX(YR**2-YI**2+XR,2.0*YR*YI+XI)) U1 = REAL(Z3) U1I = AIMAG(Z3) IF (YR * U1 + YI * U1I .GE. 0.0) GO TO 125 U1 = -U1 U1I = -U1I 125 Z3 = (CMPLX(SH,SHI) - CMPLX(XR,XI) / CMPLX(YR+U1,YI+U1I)) X / CMPLX(B3344,B3344I) SH = REAL(Z3) SHI = AIMAG(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 = ABS(A1) + ABS(A1I) + ABS(A2) U1 = A1 / S U1I = A1I / S U2 = A2 / S R = SQRT(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 = ABS(BR(K1,K1)) + ABS(BI(K1,K1)) + ABS(BR(K1,K)) U1 = BR(K1,K1) / S U1I = BI(K1,K1) / S U2 = BR(K1,K) / S R = SQRT(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 ********** IF (AI(EN,NA) .EQ. 0.0) GO TO 70 R = CABS(CMPLX(AR(EN,NA),AI(EN,NA))) 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 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 C C ------------------------------------------------------------------ C SUBROUTINE CQZVEC(NM,N,AR,AI,BR,BI,ALFR,ALFI,BETA,ZR,ZI) C INTEGER I,J,K,M,N,EN,II,JJ,NA,NM,NN REAL AR(NM,N),AI(NM,N),BR(NM,N),BI(NM,N),ALFR(N),ALFI(N), X BETA(N),ZR(NM,N),ZI(NM,N) REAL R,T,RI,TI,XI,ALMI,ALMR,BETM,EPSB REAL CABS COMPLEX Z3 COMPLEX CMPLX REAL REAL,AIMAG C C C C C 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 C C ------------------------------------------------------------------ C 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 .EQ. EN) GO TO 605 XI = T * BI(J,EN) + TI * BR(J,EN) T = T * BR(J,EN) - TI * BI(J,EN) TI = XI 605 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 = CMPLX(R,RI) / CMPLX(T,TI) BR(I,EN) = REAL(Z3) BI(I,EN) = AIMAG(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 880 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 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 = CABS(CMPLX(ZR(I,J),ZI(I,J))) 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 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 = 20. C REAL A( 20, 20),B( 20, 20),Z( 20, 20), $ ALFR( 20),ALFI( 20),BETA( 20),NORM( 20), $ RESDUL,EPS1, AI(20,20),BI(20,20),ZI(20,20), $ AHOLD( 20, 20),BHOLD( 20, 20), AHOLDI(20,20),BHOLDI(20,20) REAL AMAX1 INTEGER ERROR DATA IWRITE/6/ C NM = 20 10 CALL RMATIN(NM,N,A,B,AHOLD,BHOLD,0) CALL RMATIN(NM,N,AI,BI,AHOLDI,BHOLDI,0) WRITE(IWRITE,20) $ N 20 FORMAT(30H1THE FULL MATRIX A OF ORDER, $ I4,22H IS (PRINTED BY ROWS)/) DO 30 I = 1,N 30 WRITE(IWRITE,40) $ (A(I,J),AI(I,J),J=1,N) 40 FORMAT(5(2F6.0,1HI,3X)) WRITE(IWRITE,41) $ N 41 FORMAT(30H0THE FULL MATRIX B OF ORDER, $ I4,22H IS (PRINTED BY ROWS)/) DO 42 I = 1,N 42 WRITE(IWRITE,43) $ (B(I,J),BI(I,J),J=1,N) 43 FORMAT(5(2F6.0,1HI,3X)) C C EIGENVALUES AND EIGENVECTORS USING CQZVAL AND CQZVEC C EPS1 = 0.0 CALL CQZHES(NM,N,A,AI,B,BI,.TRUE.,Z,ZI) CALL CQZVAL(NM,N,A,AI,B,BI,EPS1,ALFR,ALFI,BETA, $ .TRUE.,Z,ZI,ERROR) CALL CQZVEC(NM,N,A,AI,B,BI,ALFR,ALFI,BETA,Z,ZI) WRITE(IWRITE,292) 292 FORMAT(/15X,7HALFR(I),19X,7HALFI(I),19X,7HBETA(I)) DO 295 I = 1,N WRITE(IWRITE,293) $ I,ALFR(I),ALFI(I),BETA(I) 293 FORMAT(I2,3(1PE23.6,3X)) 295 CONTINUE CALL RMATIN(NM,N,A,B,AHOLD,BHOLD,1) CALL RMATIN(NM,N,AI,BI,AHOLDI,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(IWRITE,300) 300 FORMAT(/14X,35HCOMPUTED EIGENVALUE AND EIGENVECTOR,20X,8HRESIDUAL) C DO 510 K = 1, N BETA(K) = AMAX1(BETA(K),1.0E-50) ALFR(K) = ALFR(K) / BETA(K) ALFI(K) = ALFI(K) / BETA(K) C C ONE EIGENVECTOR. C 340 WRITE(IWRITE,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,B,AHOLD,BHOLD,INITIL) C C THIS INPUT SUBROUTINE READS TWO REAL MATRICES A AND B 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). C REAL A(NM,NM),B(NM,NM),AHOLD(NM,NM),BHOLD(NM,NM) REAL FLOAT INTEGER IA( 20), IB( 20) DATA IREADA/1/,IREADB/2/,IWRITE/6/ C IF( INITIL .EQ. 1 ) GO TO 30 READ(IREADA,5) $ N, M 5 FORMAT(I6,6X,I6) IF( N .EQ. 0 ) GO TO 70 IF (M .NE. 1) GO TO 16 DO 10 I = 1,N READ(IREADA,17) $ (IA(J), J=I,N) DO 9 J = I,N A(I,J) = FLOAT(IA(J)) 9 A(J,I) = A(I,J) 10 CONTINUE 11 READ(IREADB,5) $ N,M IF( M .NE. 1 ) GO TO 20 DO 15 I = 1,N READ(IREADB,17) $ (IB(J), J=I,N) DO 14 J = I,N B(I,J)=FLOAT(IB(J)) 14 B(J,I)=B(I,J) 15 CONTINUE GO TO 22 16 DO 18 I = 1,N READ(IREADA,17) $ (IA(J), J=1,N) 17 FORMAT(6I12) DO 18 J = 1,N 18 A(I,J) = FLOAT(IA(J)) GO TO 11 20 DO 25 I = 1,N READ(IREADB,17) $ (IB(J),J=1,N) DO 25 J = 1,N 25 B(I,J) = FLOAT(IB(J)) 22 DO 23 I = 1,N DO 23 J = 1,N BHOLD(I,J) = B(I,J) 23 AHOLD(I,J) = A(I,J) RETURN 30 DO 40 I = 1,N DO 40 J = 1,N B(I,J) = BHOLD(I,J) 40 A(I,J) = AHOLD(I,J) RETURN 70 WRITE(IWRITE,80) 80 FORMAT(47H0END OF DATA FOR SUBROUTINE RMATIN(RSGREADI). /1H1) STOP END C C ---------------------------------------------------------------- C SUBROUTINE CGGWZR(NM,N,A,AI,B,BI,ALFR,ALFI,BETA,Z,ZI,NORM,RESDUL) C REAL 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) REAL CABS,ABS,FLOAT,AMAX1 COMPLEX C, C1 COMPLEX CMPLX EQUIVALENCE(C1,CPART(1)) C 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..)) C C ---------------------------------------------------------------- C NORMB = 0.0 NORMA = 0.0 RESDUL = 0.0 C DO 20 I = 1,N SUMA = 0.0 SUMB = 0.0 DO 10 L = 1,N SUMA = SUMA + ABS(A(L,I)) + ABS(AI(L,I)) 10 SUMB = SUMB + ABS(B(L,I)) + ABS(BI(L,I)) NORMA = AMAX1(NORMA,SUMA) 20 NORMB = AMAX1(NORMB,SUMB) C 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 + CABS(CMPLX(Z(L,I),ZI(L,I))) C 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) 100 SUMI2 = SUMI2 + A(L,K)*ZI(K,I) + AI(L,K)*Z(K,I) SUMR3 = -ALFR(I)*SUMR + ALFI(I)*SUMI SUMI3 = -ALFI(I)*SUMR - ALFR(I)*SUMI C 110 S = S+CABS(CMPLX(SUMR3,SUMI3)+CMPLX(SUMR2,SUMI2)*BETA(I)) NORMAB = NORMA*BETA(I)+NORMB*CABS(CMPLX(ALFR(I),ALFI(I))) IF( NORMAB .EQ. 0.0 ) NORMAB = 1.0 C 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(CABS(CMPLX(Z(L,I),ZI(L,I))) .GE. NORM(I)/FLOAT(N)) $ GO TO 130 120 CONTINUE C 130 XR = NORM(I)*Z(L,I)/CABS(CMPLX(Z(L,I),ZI(L,I))) XI = NORM(I)*ZI(L,I)/CABS(CMPLX(Z(L,I),ZI(L,I))) C = CMPLX(XR,XI) C DO 140 L= 1,N C1 = CMPLX(Z(L,I),ZI(L,I))/C Z(L,I) = CPART(1) 140 ZI(L,I) = CPART(2) C NORM(I) = S/(NORM(I)*NORMAB) 150 RESDUL = AMAX1(NORM(I),RESDUL) 160 CONTINUE C RETURN END 5 DATA1 1 -238 86 164 -166 56 DATA1 2 76 -96 40 60 -60 DATA1 3 118 55 -13 34 -176 DATA1 4 -314 132 114 -90 -424 DATA1 5 -54 -205 109 158 -38 DATA1 6 5 DATA1 7 -344 178 240 -308 158 DATA1 8 152 -128 -32 184 -136 DATA1 9 284 -182 460 -192 -214 DATA1 10 -160 78 296 -164 -374 DATA1 11 -24 -400 148 312 -96 DATA1 12 5 DATA1 13 41 -143 -20 20 104 DATA1 14 148 144 -6 -78 8 DATA1 15 -19 87 4 -56 -164 DATA1 16 -60 -81 99 34 84 DATA1 17 1 133 132 -46 -12 DATA1 18 5 DATA1 19 -369 -747 -1368 486 -432 DATA1 20 261 666 -1152 45 -540 DATA1 21 819 243 1548 -954 180 DATA1 22 -945 -279 171 441 -144 DATA1 23 -468 747 774 -45 -216 DATA1 24 5 DATA1 25 -15 -143 -83 41 55 DATA1 26 100 144 -60 -60 -34 DATA1 27 -19 87 4 -56 -164 DATA1 28 -116 -81 36 55 35 DATA1 29 -39 133 87 -31 -47 DATA1 30 5 DATA1 31 -1635 -173 -1142 1012 -566 DATA1 32 -829 1128 -1050 495 -710 DATA1 33 971 187 1558 -1016 218 DATA1 34 -562 -944 -1310 238 -937 DATA1 35 356 -149 -875 -434 -1015 DATA1 36 0 DATA1 37 5 DATA2 1 388 -386 -250 556 -396 DATA2 2 -304 384 -160 -240 240 DATA2 3 -658 -73 -109 -118 406 DATA2 4 -640 204 -692 288 -192 DATA2 5 -162 631 131 -758 278 DATA2 6 5 DATA2 7 94 -122 -14 130 -62 DATA2 8 -76 64 16 -92 68 DATA2 9 -136 100 -250 100 96 DATA2 10 -10 -42 -90 66 154 DATA2 11 -72 158 52 -184 76 DATA2 12 5 DATA2 13 90 180 36 -90 -72 DATA2 14 -105 -210 -42 105 84 DATA2 15 -90 -180 -36 90 72 DATA2 16 75 150 30 -75 -60 DATA2 17 -75 -150 -30 75 60 DATA2 18 5 DATA2 19 161 335 182 -162 -36 DATA2 20 -169 -322 24 167 204 DATA2 21 -211 -307 -160 186 36 DATA2 22 205 215 45 -165 -80 DATA2 23 -48 -299 -102 89 88 DATA2 24 5 DATA2 25 90 180 36 -90 -72 DATA2 26 -105 -210 -42 105 84 DATA2 27 -90 -180 -36 90 72 DATA2 28 75 150 30 -75 -60 DATA2 29 -75 -150 -30 75 60 DATA2 30 5 DATA2 31 307 253 163 -235 -39 DATA2 32 -49 -410 46 85 182 DATA2 33 -229 -253 -200 234 84 DATA2 34 171 301 184 -130 25 DATA2 35 -134 -185 59 146 195 DATA2 36