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