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