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