C C ________________________________________________________ C | | C | REDUCE A REAL SYMMETRIC BAND MATRIX TO TRIDIAGONAL FORM| C | USING GIVENS ROTATIONS | C | | C | INPUT: | C | | C | A --ARRAY CONTAINING DIAGONAL AND SUBDIAG- | C | ONAL BANDS FROM COEFFICIENT MATRIX | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | | C | N --MATRIX DIMENSION | C | | C | H --HALF BANDWIDTH | C | | C | W --WORK ARRAY WHICH CAN BE IDENTIFIED WITH| C | ARRAY A ALTHOUGH THE ORIGINAL COEFFI- | C | CIENTS WILL BE DESTROYED | C | (LENGTH AT LEAST (H+1)(N-.5H+6) - 12) | C | | C | OUTPUT: | C | | C | D --DIAGONAL OF TRIDIAGONAL MATRIX | C | | C | U --SUPERDIAGONAL OF TRIDIAGONAL MATRIX | C | | C | A --THE ORIGINAL A ARRAY IS UNTOUCHED | C | UNLESS W IS IDENTIFIED WITH A | C | | C | BUILTIN FUNCTIONS: ABS,MAX0,MIN0,SQRT | C | PACKAGE SUBROUTINES: HTR | C |________________________________________________________| C SUBROUTINE HHESS(D,U,A,LA,N,H,W) INTEGER G,H,I,J,K,L,LA,M,N,O REAL A(LA,1),D(1),U(1),W(1) M = N + 1 K = N - 1 IF ( H .GT. 0 ) GOTO 20 D(N) = A(1,N) IF ( N .EQ. 1 ) RETURN DO 10 J = 1,K D(J) = A(1,J) 10 U(J) = 0. RETURN 20 IF ( H .GT. 1 ) GOTO 40 DO 30 J = 1,K D(J) = A(1,J) 30 U(J) = A(2,J) D(N) = A(1,N) RETURN 40 G = H + 1 K = 0 DO 60 J = 1,N L = MIN0(G,M-J) DO 50 I = 1,L 50 W(I+K) = A(I,J) 60 K = K + G L = N*G + 1 DO 70 I = 1,N L = L - G K = 0 O = L + MIN0(N-I,H) DO 70 J = L,O W(J) = W(J-K) 70 K = K + G J = 0 K = 0 L = 0 DO 90 I = 1,N M = L + 1 O = MIN0(I,G) L = L + O DO 80 J = M,L 80 W(J) = W(J+K) 90 K = K + G - O G = H - 1 I = L + 1 J = I + G K = J + G L = K + G M = L + G G = M + G CALL HTR(W,W(I),W(J),W(K),W(L),W(M),W(G),D,U,N,H) RETURN END C% SUBROUTINE HTR(A,B,C,Z,BN,CN,ZN,D,UU,N,H) REAL A(1),B(1),C(1),BN(1),CN(1),D(1),UU(1) REAL E,P,Q,QQ,R,S,T,X,Y INTEGER F,G,H,HG,HW,H1,H2,I,IX,IY,J,JL,JP,K,L,M,N INTEGER NH,O,U,V,W,Z(1),ZN(1) E = 65536.**(-3) DO 10 I = 1,N 10 D(I) = 1. G = H + 1 F = G + 1 HG = H*G H1 = H - 1 H2 = H1 + H1 O = H2 - N NH = N - H W = 1 20 HW = H + W V = W W = W + 1 IF ( W .GE. N ) GOTO 410 T = 1. DO 30 I = W,N 30 IF ( D(I) .LT. T ) T = D(I) IF ( T .GT. E ) GOTO 130 C --------------------------------------------- C |*** SCALE TO AVOID OVER AND UNDER FLOWS ***| C --------------------------------------------- DO 120 I = W,N T = SQRT(D(I)) D(I) = 1. IF ( I .LT. G ) GOTO 70 J = G*I - HG/2 K = J - H 40 A(J) = T*A(J) IF ( J .EQ. K ) GOTO 50 J = J - 1 GOTO 40 50 K = J + F*MIN0(H,N-I) 60 A(J) = T*A(J) J = J + F IF ( J .LE. K ) GOTO 60 GOTO 120 70 J = (I*(I+1))/2 K = J - I + 1 80 A(J) = T*A(J) IF ( J .EQ. K ) GOTO 90 J = J - 1 GOTO 80 90 K = I 100 A(J) = T*A(J) K = K + 1 J = J + K IF ( K .LE. G ) GOTO 100 K = G - MIN0(I,N-H) IF ( K .GE. H ) GOTO 120 110 A(J) = T*A(J) J = J + F K = K + 1 IF ( K .LT. H ) GOTO 110 120 CONTINUE C ----------------------------------------- C |*** START SIMILARITY TRANSFORMATION ***| C ----------------------------------------- 130 U = HG/2 + F*MIN0(V,NH) - V - F J = U K = U + F IX = MIN0(HW,N) IF ( IX .EQ. G ) J = J + 1 M = MAX0(0,V-NH) GOTO 330 140 IF ( IX .GT. NH ) GOTO 20 IX = IY + H2 U = U + HG J = U M = MAX0(0,IX-N) IF ( M .EQ. 0 ) GOTO 150 IX = N M = M - 1 J = J - M*F K = J - F 150 M = M + 1 S = A(J) JP = J + 1 DO 170 I = 1,M T = A(JP) IF ( Z(I) .GT. 0 ) GOTO 160 A(J) = T + S*B(I) S = T*C(I) - S J = JP JP = JP + 1 GOTO 170 160 A(J) = S - T*B(I) S = T + S*C(I) J = JP JP = JP + 1 170 CONTINUE A(J) = S IF ( IX .LT. N ) GOTO 190 IF ( M-O .GT. IY ) GOTO 190 IF ( M .EQ. H1 ) GOTO 20 J = K GOTO 150 190 K = J + G IF ( Z(M) .GT. 0 ) GOTO 200 T = -A(K) A(K) = -T*B(M) GOTO 210 200 T = A(K)*C(M) 210 L = J - H 220 IY = IX IX = IX - 1 X = D(IX) Y = D(IY) IF ( ABS(S) .GE. ABS(T) ) GOTO 230 R = X/Y P = -S/T Q = R*P*P IF ( Q .LE. 1. ) GOTO 250 P = T/S QQ = Q/(1.+Q) Q = P/R GOTO 290 230 IF ( S .EQ. 0. ) GOTO 240 R = Y/X P = T/S Q = R*P*P IF ( Q .LE. 1. ) GOTO 280 P = -S/T QQ = Q/(1.+Q) Q = P/R GOTO 260 240 P = 0. Q = 0. QQ = 1. GOTO 290 250 QQ = 1./(1.+Q) Q = R*P 260 ZN(M) = 0 A(J) = Q*S - T D(IX) = QQ*Y D(IY) = QQ*X K = K - J + L JL = J - 1 J = L L = K + 1 JP = J + 1 S = A(L) A(L) = A(J) + P*S A(J) = Q*A(J) - S T = Q*S - A(K) R = S + P*A(K) A(J) = Q*A(J) - T S = A(L) A(K) = S + P*R A(L) = Q*S - R IF ( JP .GT. JL ) GOTO 310 DO 270 I = JP,JL L = L + 1 S = A(L) A(L) = A(I) + P*S 270 A(I) = Q*A(I) - S GOTO 310 280 QQ = 1./(1.+Q) Q = R*P 290 ZN(M) = 1 A(J) = S + Q*T D(IX) = QQ*X D(IY) = QQ*Y K = K - J + L JL = J - 1 J = L L = K + 1 JP = J + 1 S = A(L) A(L) = S - P*A(J) A(J) = A(J) + Q*S T = S + Q*A(K) R = A(K) - P*S A(J) = A(J) + Q*T S = A(L) A(K) = R - P*S A(L) = S + Q*R IF ( JP. GT. JL ) GOTO 310 DO 300 I = JP,JL L = L + 1 S = A(L) A(L) = S - P*A(I) 300 A(I) = A(I) + Q*S 310 BN(M) = P CN(M) = Q IF ( M .EQ. H1 ) GOTO 350 IF ( IX .LT. HW ) GOTO 320 J = J - 2 - M GOTO 150 320 IF ( IX .GT. G ) GOTO 340 J = J - V K = J + IX 330 L = J + W - IX M = M + 1 S = A(J) T = A(K) K = K - 1 GOTO 220 340 K = J + IX - V J = K - F GOTO 330 350 DO 360 I = 1,H1 B(I) = BN(I) C(I) = CN(I) 360 Z(I) = ZN(I) J = J + F L = MAX0(IY-NH,2) K = G - H M = H 370 J = J + K + M IF ( IX .GT. M ) GOTO 380 J = IY + H - M J = 2 + (J*J+J)/2 380 M = M - 1 IF ( M .LT. L ) GOTO 140 S = A(J) JP = J + 1 DO 400 I = M,H1 T = A(JP) IF ( Z(I) .GT. 0 ) GOTO 390 A(J) = T + S*B(I) S = T*C(I) - S J = JP JP = JP + 1 GOTO 400 390 A(J) = S - T*B(I) S = T + S*C(I) J = JP JP = JP + 1 400 CONTINUE A(J) = S GOTO 370 410 T = 1. D(1) = A(1) J = 1 C ------------------------------------------ C |*** MULTIPLY MATRIX BY SCALE FACTORS ***| C ------------------------------------------ DO 440 I = 2,N K = I - 1 IF ( I .GT. G ) GOTO 420 J = J + K GOTO 430 420 J = J + G 430 S = D(I) D(I) = S*A(J) S = SQRT(S) UU(K) = S*T*A(J+1) T = S 440 CONTINUE RETURN END