C C ________________________________________________________ C | | C | REDUCE A REAL SYMMETRIC BAND MATRIX TO TRIDIAGONAL FORM| C | USING GIVENS ROTATIONS AND FORM THE SIMILARITY | C | TRANSFORMATION USED IN THE REDUCTION PROCESS | C | ORIGINAL A = P TIMES TRIDIAGONAL MATRIX TIMES P INVERSE| C | | C | INPUT: | C | | C | LP --LEADING (ROW) DIMENSION OF ARRAY P | 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 | P --SIMILARITY TRANSFORMATION | 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: HTR2 | C |________________________________________________________| C SUBROUTINE HSIM(P,LP,D,U,A,LA,N,H,W) INTEGER G,H,I,J,K,L,LA,LP,M,N,O REAL A(LA,1),D(1),P(LP,1),U(1),W(1) M = N + 1 K = N - 1 IF ( H .GT. 0 ) GOTO 50 IF ( N .EQ. 1 ) GOTO 40 DO 20 J = 1,K DO 10 I = 1,N 10 P(I,J) = 0. P(J,J) = 1. D(J) = A(1,J) 20 U(J) = 0. DO 30 I = 1,N 30 P(I,N) = 0. 40 P(N,N) = 1. D(N) = A(1,N) RETURN 50 IF ( H .GT. 1 ) GOTO 90 DO 70 J = 1,K DO 60 I = 1,N 60 P(I,J) = 0. P(J,J) = 1. D(J) = A(1,J) 70 U(J) = A(2,J) DO 80 I = 1,N 80 P(I,N) = 0. P(N,N) = 1. D(N) = A(1,N) RETURN 90 G = H + 1 K = 0 DO 120 J = 1,N L = MIN0(G,M-J) DO 100 I = 1,L 100 W(I+K) = A(I,J) DO 110 I = 1,N 110 P(I,J) = 0. P(J,J) = 1. 120 K = K + G L = N*G + 1 DO 130 I = 1,N L = L - G K = 0 O = L + MIN0(N-I,H) DO 130 J = L,O W(J) = W(J-K) 130 K = K + G J = 0 K = 0 L = 0 DO 150 I = 1,N M = L + 1 O = MIN0(I,G) L = L + O DO 140 J = M,L 140 W(J) = W(J+K) 150 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 HTR2(W,W(I),W(J),W(K),W(L),W(M),W(G),D,U,N,H,P,LP) RETURN END C% SUBROUTINE HTR2(A,B,C,Z,BN,CN,ZN,D,UU,N,H,V,LV) INTEGER F,G,H,HG,HW,H1,H2,I,IX,IY,J,K,L,LV,M,N,NH INTEGER O,U,W,W0,Z(1),ZN(1) REAL A(1),B(1),C(1),BN(1),CN(1),D(1),UU(1),V(LV,1) REAL E,P,Q,QQ,R,S,T,X,Y 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 W0 = W W = W + 1 IF ( W .GE. N ) GOTO 450 T = 1. DO 30 I = W,N 30 IF ( D(I) .LT. T ) T = D(I) IF ( T .GT. E ) GOTO 140 C --------------------------------------------- C |*** SCALE TO AVOID OVER AND UNDER FLOWS ***| C --------------------------------------------- DO 130 J = W,N T = SQRT(D(J)) D(J) = 1. DO 40 I = 1,N 40 V(I,J) = T*V(I,J) IF ( J .LT. G ) GOTO 80 K = G*J - HG/2 L = K - H 50 A(K) = T*A(K) IF ( K .EQ. L ) GOTO 60 K = K - 1 GOTO 50 60 L = K + F*MIN0(H,N-J) 70 A(K) = T*A(K) K = K + F IF ( K .LE. L ) GOTO 70 GOTO 130 80 K = (J*(J+1))/2 L = K - J + 1 90 A(K) = T*A(K) IF ( K .EQ. L ) GOTO 100 K = K - 1 GOTO 90 100 L = J 110 A(K) = T*A(K) L = L + 1 K = K + L IF ( L .LE. G ) GOTO 110 L = G - MIN0(J,N-H) IF ( L .GE. H ) GOTO 130 120 A(K) = T*A(K) K = K + F L = L + 1 IF ( L .LT. H ) GOTO 120 130 CONTINUE C ----------------------------------------- C |*** START SIMILARITY TRANSFORMATION ***| C ----------------------------------------- 140 U = HG/2 + F*MIN0(W0,NH) - W0 - F J = U K = U + F IX = MIN0(HW,N) IF ( IX .EQ. G ) J = J + 1 M = MAX0(0,W0-NH) GOTO 360 150 IF ( IX .GT. NH ) GOTO 20 IX = IY + H2 U = U + HG J = U M = MAX0(0,IX-N) IF ( M .EQ. 0 ) GOTO 160 IX = N M = M - 1 J = J - M*F K = J - F 160 M = M + 1 I = 0 S = A(J) 170 IF ( I .EQ. M ) GOTO 190 I = I + 1 T = A(J+1) IF ( Z(I) .GT. 0 ) GOTO 180 A(J) = T + S*B(I) S = T*C(I) - S J = J + 1 GOTO 170 180 A(J) = S - T*B(I) S = T + S*C(I) J = J + 1 GOTO 170 190 A(J) = S IF ( IX .LT. N ) GOTO 200 IF ( M-O .GT. IY ) GOTO 200 IF ( M .EQ. H1 ) GOTO 20 J = K GOTO 160 200 K = J + G IF ( Z(I) .GT. 0 ) GOTO 210 T = -A(K) A(K) = -T*B(I) GOTO 220 210 T = A(K)*C(I) 220 L = J - H 230 IY = IX IX = IX - 1 X = D(IX) Y = D(IY) IF ( ABS(S) .GE. ABS(T) ) GOTO 240 R = X/Y P = -S/T Q = R*P*P IF ( Q .LE. 1. ) GOTO 260 P = T/S QQ = Q/(1.+Q) Q = P/R GOTO 310 240 IF ( S .EQ. 0. ) GOTO 250 R = Y/X P = T/S Q = R*P*P IF ( Q .LE. 1. ) GOTO 300 P = -S/T QQ = Q/(1.+Q) Q = P/R GOTO 270 250 P = 0. Q = 0. QQ = 1. GOTO 310 260 QQ = 1./(1.+Q) Q = R*P 270 ZN(M) = 0 A(J) = Q*S - T D(IX) = QQ*Y D(IY) = QQ*X DO 280 I = 1,N T = V(I,IX) S = V(I,IY) V(I,IX) = Q*T - S 280 V(I,IY) = T + P*S 290 J = J - 1 S = A(K) A(K) = A(J) + P*S A(J) = Q*A(J) - S K = K - 1 IF ( J .GT. L ) GOTO 290 T = Q*S - A(K) R = S + P*A(K) A(J) = Q*A(J) - T I = K + 1 S = A(I) A(K) = S + P*R A(I) = Q*S - R GOTO 340 300 QQ = 1./(1.+Q) Q = R*P 310 ZN(M) = 1 A(J) = S + Q*T D(IX) = QQ*X D(IY) = QQ*Y DO 320 I = 1,N T = V(I,IX) S = V(I,IY) V(I,IX) = T + Q*S 320 V(I,IY) = S - P*T 330 J = J - 1 S = A(K) A(K) = S - P*A(J) A(J) = A(J) + Q*S K = K - 1 IF ( J .GT. L ) GOTO 330 T = S + Q*A(K) R = A(K) - P*S A(J) = A(J) + Q*T I = K + 1 S = A(I) A(K) = R - P*S A(I) = S + Q*R 340 BN(M) = P CN(M) = Q IF ( M .EQ. H1 ) GOTO 380 IF ( IX .LT. HW ) GOTO 350 J = J - 2 - M GOTO 160 350 IF ( IX .GT. G ) GOTO 370 J = J - W0 K = J + IX 360 L = J + W - IX M = M + 1 S = A(J) T = A(K) K = K - 1 GOTO 230 370 K = J + IX - W0 J = K - F GOTO 360 380 I = 1 390 B(I) = BN(I) C(I) = CN(I) Z(I) = ZN(I) I = I + 1 IF ( I .LT. H ) GOTO 390 J = J + F L = MAX0(IY-NH,2) K = G - H M = H 400 J = J + K + M IF ( IX .GT. M ) GOTO 410 J = IY + H - M J = 2 + (J*J+J)/2 410 M = M - 1 IF ( M .LT. L ) GOTO 150 I = M S = A(J) 420 T = A(J+1) IF ( Z(I) .GT. 0 ) GOTO 430 A(J) = T + S*B(I) S = T*C(I) - S GOTO 440 430 A(J) = S - T*B(I) S = T + S*C(I) 440 J = J + 1 I = I + 1 IF ( I .LT. H ) GOTO 420 A(J) = S GOTO 400 450 T = 1. D(1) = A(1) L = 1 C ------------------------------------------ C |*** MULTIPLY MATRIX BY SCALE FACTORS ***| C ------------------------------------------ DO 480 J = 2,N K = J - 1 IF ( J .GT. G ) GOTO 460 L = L + K GOTO 470 460 L = L + G 470 S = D(J) D(J) = S*A(L) S = SQRT(S) UU(K) = S*T*A(L+1) T = S DO 480 I = 1,N 480 V(I,J) = V(I,J)*T RETURN END