C     ===== HILBERT MAIN PROGRAM =====
C
C     INVERT A 5-BY-5 HILBERT MATRIX AND THEN INVERT THE INVERSE.
C
C     TEST SUBROUTINES THAT SOLVE SYSTEMS OF LINEAR EQUATIONS
C     USING LU DECOMPOSITION (WITHOUT ITERATIVE IMPROVEMENT).
C
C     AN EXAMPLE FORTRAN IV PROGRAM WRITTEN BY RON MAK FOR CS 152.
C     REFERENCE: COMPUTER SOLUTION OF LINEAR ALGEBRAIC SYSTEMS
C                BY GEORGE FORSYTHE AND CLEVE B. MOLER
C
      DOUBLE PRECISION H(5, 5)
      INTEGER N
      DOUBLE PRECISION DENOM
C
      N = 5
C
C     COMPUTE THE HILBERT MATRIX.
      DO 200 I = 1, N
          DO 100 J = 1, N
              DENOM = I + J - 1
              H(I, J) = 1.0/DENOM
  100     CONTINUE
  200 CONTINUE
C
      WRITE (6, 300)
  300 FORMAT (/'HILBERT MATRIX'/)
      CALL PRINT (H, N)
C
      CALL INVERT(N, H, H)
C
      WRITE (6, 400)
  400 FORMAT (/'HILBERT MATRIX INVERTED:'/)
      CALL PRINT (H, N)
C
      CALL INVERT (N, H, H)
C
      WRITE (6, 500)
  500 FORMAT (/'INVERSE HILBERT MATRIX INVERTED:'/)
      CALL PRINT (H, N)
C
      PAUSE
      STOP
      END
C
C     ===== SUBROUTINE DECOMP =====
C
      SUBROUTINE DECOMP (N, A, LU)
      INTEGER N
      DOUBLE PRECISION A(5, 5), LU(5, 5), SCALES(5), PS(5)
C
      INTEGER I, J, K, PIVOTX, NM1, PSI, PSK, PSN, KP1
      DOUBLE PRECISION NRMROW, PIVOT, SIZE, BIGGST, MULT
      COMMON PS
C
C     INITIALIZE PS, LU, AND SCALES.
      DO 200 I = 1, N
          PS(I) = I
          NRMROW = 0.0
C
          DO 100 J = 1, N
              LU(I, J) = A(I, J)
C
C             FIND THE LARGEST ROW ELEMENT.
              IF (NRMROW .GE. DABS(LU(I, J))) GO TO 100
                  NRMROW = DABS(LU(I, J))
  100     CONTINUE
C
C         SET THE SCALING FACTOR FOR ROW EQUILIBRATION.
          IF (NRMROW .NE. 0) GO TO 110
              SCALES(I) = 0.0
              CALL SINGLR(0)
              GO TO 200
C
  110     SCALES(I) = 1.0/NRMROW
  200 CONTINUE
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING.
C     PIVOT ROW K.
      NM1 = N-1
      DO 400 K = 1, NM1
          PIVOTX = 0
          BIGGST = 0.0
C
C         GO DOWN ROWS FROM ROW K.
          DO 300 I = K, N
              PSI = PS(I)
C
C             DIVIDE BY THE BIGGEST ROW ELEMENT.
              SIZE = DABS(LU(PSI, K))*SCALES(PSI)
C
              IF (BIGGST .GE. SIZE) GO TO 300
                  BIGGST = SIZE
                  PIVOTX = I
  300     CONTINUE
C
          IF (BIGGST .NE. 0.0) GO TO 310
              CALL SINGLR(1)
              GO TO 400
C
  310     IF (PIVOTX .EQ. K) GO TO 350
C
C             EXCHANGE ROWS.
              J = PS(K)
              PS(K) = PS(PIVOTX)
              PS(PIVOTX) = J
C
C         PIVOT ELEMENT.
  350     PSK = PS(K)
          PIVOT = LU(PSK, K)
C
C         GO DOWN REST OF ROWS.
          KP1 = K+1
          DO 380 I = KP1, N
              PSI = PS(I)
              MULT = LU(PSI, K)/PIVOT
              LU(PSI, K) = MULT
C
              IF (MULT .EQ. 0.0) GO TO 380
C
C                 INNER LOOP. ONLY THE COLUMN SUBSCRIPT VARIES.
                  DO 360 J = KP1, N
                      LU(PSI, J) = LU(PSI, J) - MULT*LU(PSK, J)
  360             CONTINUE
  380     CONTINUE
  400 CONTINUE
C
C     CHECK THE BOTTOM RIGHT ELEMENT OF THE PERMUTED MATRIX.
      PSN = PS(N)
      IF (LU(PSN, N) .EQ. 0.0) CALL SINGLR(2)
c
      RETURN 
      END
C
C     ===== SUBROUTINE SOLVE =====
C
      SUBROUTINE SOLVE(N, LU, B, X)
      INTEGER N
      DOUBLE PRECISION LU(5, 5), B(5), X(5), PS(5)
C
      INTEGER I, J, IM1, PSI, NP1, IBACK, IP1
      DOUBLE PRECISION DOT
      COMMON PS
C
C     LY = B : SOLVE FOR Y.
      DO 200 I = 1, N
          IM1 = I-1
          PSI = PS(I)
C
          DOT = 0.0
          DO 100 J = 1, IM1
              DOT = DOT + LU(PSI, J)*X(J)
  100     CONTINUE
          X(I) = B(PSI) - DOT
  200 CONTINUE
C
C     UX = Y : SOLVE FOR X
      NP1 = N+1
      DO 400 IBACK = 1, N
          I = NP1-IBACK
C         I = N, N-1, N-2, ..., 2, 1
C
          IP1 = I+1
          PSI = PS(I)
C
          DOT = 0.0
          DO 300 J = IP1, N
              DOT = DOT + LU(PSI, J)*X(J)
  300     CONTINUE
          X(I) = (X(I) - DOT)/LU(PSI, I)
  400 CONTINUE
C
      RETURN
      END
C
C     ===== SUBROUTINE INVERT =====
C
      SUBROUTINE INVERT(N, A, AINV)
      INTEGER N
      DOUBLE PRECISION A(5, 5), AINV(5, 5)
C
C     COMPUTE AINV = INVERSE(A). NOTE THAT A AND AINV
C     ARE OFTEN PASSED THE SAME MATRIX.
C
      DOUBLE PRECISION LU(5, 5), B(5), X(5)
      INTEGER I, J
C
      CALL DECOMP(N, A, LU)
C
      DO 300 J = 1, N
          DO 100 I = 1, N
              B(I) = 0.0
              IF (I .EQ. J) B(I) = 1.0
  100     CONTINUE
c
          CALL SOLVE(N, LU, B, X)
C
          DO 200 I = 1, N
              AINV(I, J) = X(I)
  200     CONTINUE
  300 CONTINUE
C
      RETURN
      END
C
C     ===== SUBROUTINE SINGLR =====
C
      SUBROUTINE SINGLR(WHY)
      INTEGER WHY
C
      GO TO (1, 2, 3), WHY
C
    1 WRITE (6, 10)
   10 FORMAT ('MATRIX WITH ZERO ROW IN DECOMPOSE.')
      RETURN
C
    2 WRITE (6, 20)
   20 FORMAT ('SINGULAR MATRIX IN DECOMPOSE.')
      RETURN
C
    3 WRITE (6, 30)
   30 FORMAT ('NO CONVERGENCE IN IMPROVE.')
      RETURN
C
      END
C
C     ===== SUBROUTINE PRINT =====
C
      SUBROUTINE PRINT (A, N)
      INTEGER N
      DOUBLE PRECISION A(5, 5)
C
      INTEGER I, J
C
      DO 100 I = 1, N
          WRITE (6, 10) (A(I, J), J = 1, N)
   10     FORMAT(5F15.6)
  100 CONTINUE
C
      RETURN
      END