C
C      ________________________________________________________
C     |                                                        |
C     | BALANCE A REAL MATRIX AND REDUCE IT TO HESSENBERG FORM |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         A     --REAL ARRAY CONTAINING MATRIX           |
C     |                 (LENGTH AT LEAST 1 + N(N+2))           |
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: BAL,PACK                        |
C     |________________________________________________________|
C
      SUBROUTINE AHESS(A,LA,N,W)
      REAL A(1),W(1),R,S,T,U,V
      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 = I + 1
      CALL BAL(A,N,N,W,A(J))
      DO 10 K = 1,N
10         A(E+K) = W(K)
      V = W(1)
      J = M
      K = I
C     ---------------------------
C     |*** SHIFT MATRIX DOWN ***|
C     ---------------------------
20    K = K - N
30    A(I+J) = A(I)
      I = I - 1
      IF ( I .GT. K ) GOTO 30
      J = J - 1
      IF ( K .GT. 0 ) GOTO 20
      A(1) = 2231
      A(2) = N
      K = 4
      L = O
      D = 1
      C = 2
40    IF ( C .GE. N ) GOTO 200
      P = K + 1
      DO 50 I = P,L
50         IF ( A(I) .NE. 0. ) GOTO 60
      A(L+1) = 0.
      GOTO 190
60    T = ABS(A(K))
      IF ( T .NE. 0. ) U = 1./T
      R = 1.
      DO 80 J = I,L
           S = ABS(A(J))
           IF ( S .LE. T ) GOTO 70
           U = 1./S
           R = 1. + R*(T*U)**2
           T = S
           GOTO 80
70         R = R + (S*U)**2
80    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     ------------------------------------
90    A(I+1) = U*A(I)
      I = I - 1
      IF ( I .GT. K ) GOTO 90
      A(K) = -S
      A(P) = U*(R+S)
      H = L
      DO 100 I = 1,N
100        W(I) = 0.
C     --------------------------------------
C     |*** MULTIPLY FROM RIGHT AND LEFT ***|
C     --------------------------------------
110   H = H + M
      S = A(P)
      P = P + 1
      Q = H - N
      DO 120 I = 1,D
120        W(I) = W(I) + S*A(I+Q)
      J = K - D
      T = 0.
      DO 130 I = C,N
           R = A(I+Q)
           T = T + R*A(I+J)
130        W(I) = W(I) + R*S
      A(H+1) = T
      IF ( H .LT. E ) GOTO 110
      T = 0.
      H = L + 1
      P = K + 1
      J = C - P
      DO 140 I = P,H
140        T = T + W(I+J)*A(I)
      DO 150 I = C,N
150        W(I) = W(I) - T*A(I-J)
      H = L
C     -----------------------------------
C     |*** UPDATE COEFFICIENT MATRIX ***|
C     -----------------------------------
160   G = H + 2
      Q = H + M
      H = H + C
      T = A(Q+1)
      S = A(P)
      P = P + 1
      J = 1 - G
      DO 170 I = G,H
170        A(I) = A(I) - W(I+J)*S
      I = H
      H = Q
      Q = K - I
      G = I + 1
      DO 180 I = G,H
180        A(I) = A(I) - A(I+Q)*T - W(I+J)*S
      IF ( H .LT. E ) GOTO 160
190   K = K + O
      L = L + M
      D = C
      C = C + 1
      GOTO 40
200   A(E+1) = V
      RETURN
      END