C C ________________________________________________________ C | | C | REDUCE A SYMMETRIC MATRIX TO TRIDIAGONAL FORM USING | C | ORTHOGONAL SIMILARITY TRANSFORMATIONS | C | | C | INPUT: | C | | C | A --COEFFICIENT MATRIX IN COMPRESSED FORMAT| C | | C | N --MATRIX DIMENSION | C | | C | OUTPUT: | C | | C | D --DIAGONAL OF TRIDIAGONAL MATRIX | C | | C | U --SUPERDIAGONAL OF TRIDIAGONAL MATRIX | C | | C | BUILTIN FUNCTIONS: ABS,SQRT | C |________________________________________________________| C SUBROUTINE SHESS(D,U,A,N) REAL A(1),D(1),U(1),R,S,T,V,W INTEGER I,J,K,L,M,N,O,P,Q D(1) = A(1) IF ( N .GT. 2 ) GOTO 20 A(1) = 2234 IF ( N .EQ. 1 ) GOTO 10 U(1) = A(2) D(2) = A(3) A(4) = D(2) A(5) = U(1) 10 A(2) = N A(3) = D(1) RETURN 20 M = N - 1 Q = 1 K = 2 L = N 30 IF ( Q .GE. M ) GOTO 150 W = A(K) P = K + 1 DO 40 I = P,L 40 IF ( A(I) .NE. 0. ) GOTO 50 U(Q) = W K = L + 2 I = Q Q = Q + 1 D(Q) = A(L+1) L = L + N - I GOTO 30 50 T = ABS(W) IF ( T .NE. 0. ) V = 1./T R = 1. DO 70 J = I,L S = ABS(A(J)) IF ( S .LE. T ) GOTO 60 V = 1./S R = 1. + R*(T*V)**2 T = S GOTO 70 60 R = R + (S*V)**2 70 CONTINUE S = T*SQRT(R) V = 1./SQRT(S*(S+ABS(W))) IF ( W .LT. 0. ) S = -S U(Q) = -S A(K) = V*(W+S) Q = Q + 1 D(Q) = A(K) U(Q) = 0. O = Q - K C -------------------------------------- C |*** CONSTRUCT HOUSEHOLDER MATRIX ***| C -------------------------------------- DO 80 I = P,L A(I) = V*A(I) J = O + I D(J) = A(I) 80 U(J) = 0. I = Q L = L + 1 O = L - Q S = 0. C ---------------------------------------- C |*** MULTIPLY BY HOUSEHOLDER MATRIX ***| C ---------------------------------------- 90 V = D(I) T = V*A(I+O) IF ( I .GE. N ) GOTO 110 P = I + 1 DO 100 J = P,N W = A(J+O) T = T + D(J)*W 100 U(J) = U(J) + V*W U(I) = U(I) + T S = S + U(I)*V O = O + N - I I = P GOTO 90 110 U(N) = U(N) + T S = .5*(S+U(N)*V) DO 120 I = Q,N 120 U(I) = S*D(I) - U(I) O = L - Q C -------------------------------------------- C |*** PERFORM HOUSEHOLDER TRANSFORMATION ***| C -------------------------------------------- DO 140 J = Q,N T = U(J) S = D(J) DO 130 I = J,N P = I + O 130 A(P) = A(P) + D(I)*T + U(I)*S 140 O = O + N - J D(Q) = A(L) K = L + 1 L = L + N - Q GOTO 30 150 U(Q) = A(L) D(Q+1) = A(L+1) K = N - 1 M = 1 L = K O = 1 C -------------------------------------- C |*** COMPRESS HOUSEHOLDER VECTORS ***| C -------------------------------------- 160 DO 170 I = M,L 170 A(I) = A(I+O) M = L + 1 K = K - 1 L = L + K O = O + 1 IF ( K .GT. 1 ) GOTO 160 J = (N*(N-1))/2 + 1 I = J 180 A(I) = A(I-2) I = I - 1 IF ( I .GT. 2 ) GOTO 180 A(1) = 2234 A(2) = N RETURN END