C C ________________________________________________________ C | | C | REDUCE A COMPLEX MATRIX TO HESSENBERG FORM | C | | C | INPUT: | C | | C | A --COMPLEX ARRAY CONTAINING MATRIX | C | (LENGTH AT LEAST 2 + N(N+1)) | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | | C | N --DIMENSION OF MATRIX STORED IN A | C | | C | W --WORK ARRAY WITH AT LEAST N COMPLEX | C | ELEMENTS | C | | C | OUTPUT: | C | | C | A --HESSENBERG MATRIX | C | | C | BUILTIN FUNCTIONS: CABS,CONJG,SQRT | C | PACKAGE FUNCTIONS: MAG,SQR | C | PACKAGE SUBROUTINES: CPACK | C |________________________________________________________| C SUBROUTINE CHESS(A,LA,N,W) REAL R,S,SQR,T,U,MAG COMPLEX A(1),W(1),X,Y,Z INTEGER C,D,E,G,H,I,J,K,L,LA,M,N,O,P,Q IF ( LA .GT. N ) CALL CPACK(A,LA,N) I = N*N U = 0. M = N + 1 O = M + 1 E = I + M J = M K = I C --------------------------- C |*** SHIFT MATRIX DOWN ***| C --------------------------- 10 K = K - N 20 A(I+J) = A(I) I = I - 1 IF ( I .GT. K ) GOTO 20 J = J - 1 IF ( K .GT. 0 ) GOTO 10 A(1) = 2232 A(2) = N K = 4 L = O D = 1 C = 2 30 IF ( C .GE. N ) RETURN P = K + 1 DO 40 I = P,L 40 IF ( MAG(A(I)) .NE. 0. ) GOTO 50 A(L+1) = (0.,0.) GOTO 180 50 T = MAG(A(K)) IF ( T .NE. 0. ) U = 1./T R = SQR(A(K),U) DO 70 J = I,L S = MAG(A(J)) IF ( S .LE. T ) GOTO 60 U = 1./S R = SQR(A(J),U) + R*(T*U)**2 T = S GOTO 70 60 R = R + SQR(A(J),U) 70 CONTINUE S = T*SQRT(R) Z = A(K) T = CABS(Z) U = 1./SQRT(S*(S+T)) IF ( T .NE. 0. ) Z = Z/T IF ( T .EQ. 0. ) Z = (1.,0.) I = L C ------------------------------------ C |*** COMPUTE HOUSEHOLDER MATRIX ***| C ------------------------------------ 80 A(I+1) = U*CONJG(A(I)) I = I - 1 IF ( I .GT. K ) GOTO 80 A(K) = -Z*S A(P) = CONJG(Z)*U*(T+S) H = L DO 90 I = 1,N 90 W(I) = (0.,0.) C -------------------------------------- C |*** MULTIPLY FROM RIGHT AND LEFT ***| C -------------------------------------- 100 H = H + M Y = CONJG(A(P)) P = P + 1 Q = H - N DO 110 I = 1,D 110 W(I) = W(I) + Y*A(I+Q) J = K - D Z = (0.,0.) DO 120 I = C,N X = A(I+Q) Z = Z + X*A(I+J) 120 W(I) = W(I) + X*Y A(H+1) = Z IF ( H .LT. E ) GOTO 100 Z = (0.,0.) H = L + 1 P = K + 1 J = C - P DO 130 I = P,H Z = Z + W(I+J)*A(I) 130 A(I) = CONJG(A(I)) DO 140 I = C,N 140 W(I) = W(I) - Z*A(I-J) H = L C ----------------------------------- C |*** UPDATE COEFFICIENT MATRIX ***| C ----------------------------------- 150 G = H + 2 Q = H + M H = H + C Z = A(Q+1) Y = CONJG(A(P)) P = P + 1 J = 1 - G DO 160 I = G,H 160 A(I) = A(I) - W(I+J)*Y I = H H = Q Q = K - I G = I + 1 DO 170 I = G,H 170 A(I) = A(I) - A(I+Q)*Z - W(I+J)*Y IF ( H .LT. E ) GOTO 150 180 K = K + O L = L + M D = C C = C + 1 GOTO 30 END