*DECK QZVEC SUBROUTINE QZVEC (NM, N, A, B, ALFR, ALFI, BETA, Z) C***BEGIN PROLOGUE QZVEC C***PURPOSE The optional fourth step of the QZ algorithm for C generalized eigenproblems. Accepts a matrix in C quasi-triangular form and another in upper triangular C and computes the eigenvectors of the triangular problem C and transforms them back to the original coordinates C Usually preceded by QZHES, QZIT, and QZVAL. C***LIBRARY SLATEC (EISPACK) C***CATEGORY D4C3 C***TYPE SINGLE PRECISION (QZVEC-S) C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK C***AUTHOR Smith, B. T., et al. C***DESCRIPTION C C This subroutine is the optional fourth step of the QZ algorithm C 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 REAL matrices, one of them in C quasi-triangular form (in which each 2-by-2 block corresponds to C a pair of complex eigenvalues) and the other in upper triangular C form. It computes the eigenvectors of the triangular problem and C transforms the results back to the original coordinate system. C It is usually preceded by QZHES, QZIT, and QZVAL. C C On Input C C NM must be set to the row dimension of the two-dimensional C array parameters, A, B, and Z, as declared in the calling C program dimension statement. NM is an INTEGER variable. C C N is the order of the matrices A and B. N is an INTEGER C variable. N must be less than or equal to NM. C C A contains a real upper quasi-triangular matrix. A is a two- C dimensional REAL array, dimensioned A(NM,N). C C B contains a real upper triangular matrix. In addition, C location B(N,1) contains the tolerance quantity (EPSB) C computed and saved in QZIT. B is a two-dimensional REAL C array, dimensioned B(NM,N). C C ALFR, ALFI, and BETA are one-dimensional REAL arrays with C components whose ratios ((ALFR+I*ALFI)/BETA) are the C generalized eigenvalues. They are usually obtained from C QZVAL. They are dimensioned ALFR(N), ALFI(N), and BETA(N). C C Z contains the transformation matrix produced in the reductions C by QZHES, QZIT, and QZVAL, if performed. If the C eigenvectors of the triangular problem are desired, Z must C contain the identity matrix. Z is a two-dimensional REAL C array, dimensioned Z(NM,N). C C On Output C C A is unaltered. Its subdiagonal elements provide information C about the storage of the complex eigenvectors. C C B has been destroyed. C C ALFR, ALFI, and BETA are unaltered. C C Z contains the real and imaginary parts of the eigenvectors. C If ALFI(J) .EQ. 0.0, the J-th eigenvalue is real and C the J-th column of Z contains its eigenvector. C If ALFI(J) .NE. 0.0, the J-th eigenvalue is complex. C If ALFI(J) .GT. 0.0, the eigenvalue is the first of C a complex pair and the J-th and (J+1)-th columns C of Z contain its eigenvector. C If ALFI(J) .LT. 0.0, the eigenvalue is the second of C a complex pair and the (J-1)-th and J-th columns C of Z contain the conjugate of its eigenvector. C Each eigenvector is normalized so that the modulus C 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***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- C system Routines - EISPACK Guide, Springer-Verlag, C 1976. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 760101 DATE WRITTEN C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE QZVEC C INTEGER I,J,K,M,N,EN,II,JJ,NA,NM,NN,ISW,ENM2 REAL A(NM,*),B(NM,*),ALFR(*),ALFI(*),BETA(*),Z(NM,*) REAL D,Q,R,S,T,W,X,Y,DI,DR,RA,RR,SA,TI,TR,T1,T2 REAL W1,X1,ZZ,Z1,ALFM,ALMI,ALMR,BETM,EPSB C C***FIRST EXECUTABLE STATEMENT QZVEC EPSB = B(N,1) ISW = 1 C .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... DO 800 NN = 1, N EN = N + 1 - NN NA = EN - 1 IF (ISW .EQ. 2) GO TO 795 IF (ALFI(EN) .NE. 0.0E0) GO TO 710 C .......... REAL VECTOR .......... M = EN B(EN,EN) = 1.0E0 IF (NA .EQ. 0) GO TO 800 ALFM = ALFR(M) BETM = BETA(M) C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... DO 700 II = 1, NA I = EN - II W = BETM * A(I,I) - ALFM * B(I,I) R = 0.0E0 C DO 610 J = M, EN 610 R = R + (BETM * A(I,J) - ALFM * B(I,J)) * B(J,EN) C IF (I .EQ. 1 .OR. ISW .EQ. 2) GO TO 630 IF (BETM * A(I,I-1) .EQ. 0.0E0) GO TO 630 ZZ = W S = R GO TO 690 630 M = I IF (ISW .EQ. 2) GO TO 640 C .......... REAL 1-BY-1 BLOCK .......... T = W IF (W .EQ. 0.0E0) T = EPSB B(I,EN) = -R / T GO TO 700 C .......... REAL 2-BY-2 BLOCK .......... 640 X = BETM * A(I,I+1) - ALFM * B(I,I+1) Y = BETM * A(I+1,I) Q = W * ZZ - X * Y T = (X * S - ZZ * R) / Q B(I,EN) = T IF (ABS(X) .LE. ABS(ZZ)) GO TO 650 B(I+1,EN) = (-R - W * T) / X GO TO 690 650 B(I+1,EN) = (-S - Y * T) / ZZ 690 ISW = 3 - ISW 700 CONTINUE C .......... END REAL VECTOR .......... GO TO 800 C .......... COMPLEX VECTOR .......... 710 M = NA ALMR = ALFR(M) ALMI = ALFI(M) BETM = BETA(M) C .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT C EIGENVECTOR MATRIX IS TRIANGULAR .......... Y = BETM * A(EN,NA) B(NA,NA) = -ALMI * B(EN,EN) / Y B(NA,EN) = (ALMR * B(EN,EN) - BETM * A(EN,EN)) / Y B(EN,NA) = 0.0E0 B(EN,EN) = 1.0E0 ENM2 = NA - 1 IF (ENM2 .EQ. 0) GO TO 795 C .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... DO 790 II = 1, ENM2 I = NA - II W = BETM * A(I,I) - ALMR * B(I,I) W1 = -ALMI * B(I,I) RA = 0.0E0 SA = 0.0E0 C DO 760 J = M, EN X = BETM * A(I,J) - ALMR * B(I,J) X1 = -ALMI * B(I,J) RA = RA + X * B(J,NA) - X1 * B(J,EN) SA = SA + X * B(J,EN) + X1 * B(J,NA) 760 CONTINUE C IF (I .EQ. 1 .OR. ISW .EQ. 2) GO TO 770 IF (BETM * A(I,I-1) .EQ. 0.0E0) GO TO 770 ZZ = W Z1 = W1 R = RA S = SA ISW = 2 GO TO 790 770 M = I IF (ISW .EQ. 2) GO TO 780 C .......... COMPLEX 1-BY-1 BLOCK .......... TR = -RA TI = -SA 773 DR = W DI = W1 C .......... COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR,DI) .......... 775 IF (ABS(DI) .GT. ABS(DR)) GO TO 777 RR = DI / DR D = DR + DI * RR T1 = (TR + TI * RR) / D T2 = (TI - TR * RR) / D GO TO (787,782), ISW 777 RR = DR / DI D = DR * RR + DI T1 = (TR * RR + TI) / D T2 = (TI * RR - TR) / D GO TO (787,782), ISW C .......... COMPLEX 2-BY-2 BLOCK .......... 780 X = BETM * A(I,I+1) - ALMR * B(I,I+1) X1 = -ALMI * B(I,I+1) Y = BETM * A(I+1,I) TR = Y * RA - W * R + W1 * S TI = Y * SA - W * S - W1 * R DR = W * ZZ - W1 * Z1 - X * Y DI = W * Z1 + W1 * ZZ - X1 * Y IF (DR .EQ. 0.0E0 .AND. DI .EQ. 0.0E0) DR = EPSB GO TO 775 782 B(I+1,NA) = T1 B(I+1,EN) = T2 ISW = 1 IF (ABS(Y) .GT. ABS(W) + ABS(W1)) GO TO 785 TR = -RA - X * B(I+1,NA) + X1 * B(I+1,EN) TI = -SA - X * B(I+1,EN) - X1 * B(I+1,NA) GO TO 773 785 T1 = (-R - ZZ * B(I+1,NA) + Z1 * B(I+1,EN)) / Y T2 = (-S - ZZ * B(I+1,EN) - Z1 * B(I+1,NA)) / Y 787 B(I,NA) = T1 B(I,EN) = T2 790 CONTINUE C .......... END COMPLEX VECTOR .......... 795 ISW = 3 - ISW 800 CONTINUE C .......... END BACK SUBSTITUTION. C TRANSFORM TO ORIGINAL COORDINATE SYSTEM. C FOR J=N STEP -1 UNTIL 1 DO -- .......... DO 880 JJ = 1, N J = N + 1 - JJ C DO 880 I = 1, N ZZ = 0.0E0 C DO 860 K = 1, J 860 ZZ = ZZ + Z(I,K) * B(K,J) C Z(I,J) = ZZ 880 CONTINUE C .......... NORMALIZE SO THAT MODULUS OF LARGEST C COMPONENT OF EACH VECTOR IS 1. C (ISW IS 1 INITIALLY FROM BEFORE) .......... DO 950 J = 1, N D = 0.0E0 IF (ISW .EQ. 2) GO TO 920 IF (ALFI(J) .NE. 0.0E0) GO TO 945 C DO 890 I = 1, N IF (ABS(Z(I,J)) .GT. D) D = ABS(Z(I,J)) 890 CONTINUE C DO 900 I = 1, N 900 Z(I,J) = Z(I,J) / D C GO TO 950 C 920 DO 930 I = 1, N R = ABS(Z(I,J-1)) + ABS(Z(I,J)) IF (R .NE. 0.0E0) R = R * SQRT((Z(I,J-1)/R)**2 1 +(Z(I,J)/R)**2) IF (R .GT. D) D = R 930 CONTINUE C DO 940 I = 1, N Z(I,J-1) = Z(I,J-1) / D Z(I,J) = Z(I,J) / D 940 CONTINUE C 945 ISW = 3 - ISW 950 CONTINUE C RETURN END