C C ________________________________________________________ C | | C | REDUCE A REAL MATRIX TO HESSENBERG FORM | C | | C | INPUT: | C | | C | A --REAL 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 ELEMENTS | C | | C | OUTPUT: | C | | C | A --HESSENBERG MATRIX | C | | C | BUILTIN FUNCTIONS: ABS,SQRT | C | PACKAGE SUBROUTINES: PACK | C |________________________________________________________| C SUBROUTINE HESS(A,LA,N,W) REAL A(1),W(1),R,S,T,U INTEGER C,D,E,G,H,I,J,K,L,LA,M,N,O,P,Q IF ( LA .GT. N ) CALL PACK(A,LA,N) I = N*N 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) = 2230 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 ( A(I) .NE. 0. ) GOTO 50 A(L+1) = 0. GOTO 180 50 T = ABS(A(K)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 70 J = I,L S = ABS(A(J)) IF ( S .LE. T ) GOTO 60 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 70 60 R = R + (S*U)**2 70 CONTINUE S = T*SQRT(R) R = A(K) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S I = L C ------------------------------------ C |*** COMPUTE HOUSEHOLDER MATRIX ***| C ------------------------------------ 80 A(I+1) = U*A(I) I = I - 1 IF ( I .GT. K ) GOTO 80 A(K) = -S A(P) = U*(R+S) H = L DO 90 I = 1,N 90 W(I) = 0. C -------------------------------------- C |*** MULTIPLY FROM RIGHT AND LEFT ***| C -------------------------------------- 100 H = H + M S = A(P) P = P + 1 Q = H - N DO 110 I = 1,D 110 W(I) = W(I) + S*A(I+Q) J = K - D T = 0. DO 120 I = C,N R = A(I+Q) T = T + R*A(I+J) 120 W(I) = W(I) + R*S A(H+1) = T IF ( H .LT. E ) GOTO 100 T = 0. H = L + 1 P = K + 1 J = C - P DO 130 I = P,H 130 T = T + W(I+J)*A(I) DO 140 I = C,N 140 W(I) = W(I) - T*A(I-J) H = L C ----------------------------------- C |*** UPDATE COEFFICIENT MATRIX ***| C ----------------------------------- 150 G = H + 2 Q = H + M H = H + C T = A(Q+1) S = A(P) P = P + 1 J = 1 - G DO 160 I = G,H 160 A(I) = A(I) - W(I+J)*S I = H H = Q Q = K - I G = I + 1 DO 170 I = G,H 170 A(I) = A(I) - A(I+Q)*T - W(I+J)*S IF ( H .LT. E ) GOTO 150 180 K = K + O L = L + M D = C C = C + 1 GOTO 30 END