C
C      ________________________________________________________
C     |                                                        |
C     |         REDUCE A REAL MATRIX TO HESSENBERG FORM        |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         A     --REAL ARRAY CONTAINING MATRIX           |
C     |                 (LENGTH AT LEAST 2 + N(N+1))           |
C     |                                                        |
C     |         LA    --LEADING (ROW) DIMENSION OF ARRAY A     |
C     |                                                        |
C     |         N     --DIMENSION OF MATRIX STORED IN A        |
C     |                                                        |
C     |         W     --WORK ARRAY WITH AT LEAST N ELEMENTS    |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         A     --HESSENBERG MATRIX                      |
C     |                                                        |
C     |   BUILTIN FUNCTIONS: ABS,SQRT                          |
C     |   PACKAGE SUBROUTINES: PACK                            |
C     |________________________________________________________|
C
      SUBROUTINE HESS(A,LA,N,W)
      REAL A(1),W(1),R,S,T,U
      INTEGER C,D,E,G,H,I,J,K,L,LA,M,N,O,P,Q
      IF ( LA .GT. N ) CALL PACK(A,LA,N)
      I = N*N
      M = N + 1
      O = M + 1
      E = I + M
      J = M
      K = I
C     ---------------------------
C     |*** SHIFT MATRIX DOWN ***|
C     ---------------------------
10    K = K - N
20    A(I+J) = A(I)
      I = I - 1
      IF ( I .GT. K ) GOTO 20
      J = J - 1
      IF ( K .GT. 0 ) GOTO 10
      A(1) = 2230
      A(2) = N
      K = 4
      L = O
      D = 1
      C = 2
30    IF ( C .GE. N ) RETURN
      P = K + 1
      DO 40 I = P,L
40         IF ( A(I) .NE. 0. ) GOTO 50
      A(L+1) = 0.
      GOTO 180
50    T = ABS(A(K))
      IF ( T .NE. 0. ) U = 1./T
      R = 1.
      DO 70 J = I,L
           S = ABS(A(J))
           IF ( S .LE. T ) GOTO 60
           U = 1./S
           R = 1. + R*(T*U)**2
           T = S
           GOTO 70
60         R = R + (S*U)**2
70    CONTINUE
      S = T*SQRT(R)
      R = A(K)
      U = 1./SQRT(S*(S+ABS(R)))
      IF ( R .LT. 0. ) S = -S
      I = L
C     ------------------------------------
C     |*** COMPUTE HOUSEHOLDER MATRIX ***|
C     ------------------------------------
80    A(I+1) = U*A(I)
      I = I - 1
      IF ( I .GT. K ) GOTO 80
      A(K) = -S
      A(P) = U*(R+S)
      H = L
      DO 90 I = 1,N
90         W(I) = 0.
C     --------------------------------------
C     |*** MULTIPLY FROM RIGHT AND LEFT ***|
C     --------------------------------------
100   H = H + M
      S = A(P)
      P = P + 1
      Q = H - N
      DO 110 I = 1,D
110        W(I) = W(I) + S*A(I+Q)
      J = K - D
      T = 0.
      DO 120 I = C,N
           R = A(I+Q)
           T = T + R*A(I+J)
120        W(I) = W(I) + R*S
      A(H+1) = T
      IF ( H .LT. E ) GOTO 100
      T = 0.
      H = L + 1
      P = K + 1
      J = C - P
      DO 130 I = P,H
130        T = T + W(I+J)*A(I)
      DO 140 I = C,N
140        W(I) = W(I) - T*A(I-J)
      H = L
C     -----------------------------------
C     |*** UPDATE COEFFICIENT MATRIX ***|
C     -----------------------------------
150   G = H + 2
      Q = H + M
      H = H + C
      T = A(Q+1)
      S = A(P)
      P = P + 1
      J = 1 - G
      DO 160 I = G,H
160        A(I) = A(I) - W(I+J)*S
      I = H
      H = Q
      Q = K - I
      G = I + 1
      DO 170 I = G,H
170        A(I) = A(I) - A(I+Q)*T - W(I+J)*S
      IF ( H .LT. E ) GOTO 150
180   K = K + O
      L = L + M
      D = C
      C = C + 1
      GOTO 30
      END