      PROGRAM JACOBI
************************************************************************
* Program to find an approximate solution of linear system AX = B      *
* using the Jacobi method.  Iteration continues until either the limit *
* NUMITS on the number of iterations is exceeded  or two successive    *
* approximations differ in each component by at most EPSIL.            *
* Identifiers used are:                                                *
*     LIM     : limit on number of equations                           *
*     NUMEQN  : the number of equations (= the number of unknowns)     *
*     A       : the matrix of coefficients                             *
*     B       : the constant vector                                    *
*     NUMITS  : limit on the number of iterations                      *
*     EPSIL   : used in termination criterion                          *
*     DONE    : signals termination condition is satisfied             *
*     XOLD    : the old approximation                                  *
*     XNEW    : the new approximation                                  *
*     SUM     : summation used in iteration formula                    *
*     N       : number of iterations                                   *
*     I, J    : subscripts                                             *
*                                                                      *
* Input:   NUMEQN, A, B, NUMITS, EPSIL, initial approximation XOLD     *
* Output:  A sequence of approximate solutions or a message            *
*          indicating the limit on the number of iterations was        *
*          exceeded                                                    *
************************************************************************


      INTEGER LIM, NUMEQN, NUMITS, N, I, J
      PARAMETER (LIM = 10)
      REAL A(LIM,LIM), B(LIM), XOLD(LIM), XNEW(LIM), EPSIL, SUM
      LOGICAL DONE

      PRINT *, 'ENTER THE NUMBER OF EQUATIONS:'
      READ *, NUMEQN

      PRINT *, 'ENTER THE COEFFICIENT MATRIX (ROWWISE):'
      DO 10 I = 1, NUMEQN
          READ *, (A(I,J), J = 1, NUMEQN)
10    CONTINUE

      PRINT *, 'ENTER THE CONSTANT VECTOR:'
      READ *, (B(I), I = 1, NUMEQN)

      PRINT *, 'ENTER THE TOLERANCE AND ITERATION LIMIT:'
      READ *, EPSIL, NUMITS

      PRINT *, 'ENTER AN INITIAL APPROXIMATE SOLUTION:'
      READ *, (XOLD(I), I = 1, NUMEQN)

      DONE = .FALSE.
      N = 0

* While NOT DONE do the following:
20    IF (.NOT. DONE) THEN

*     Generate new approximate solution XNEW
          N = N + 1
          DO 50 I = 1, NUMEQN
              SUM = B(I)
              DO 30 J = 1, I - 1
                  SUM = SUM - A(I,J) * XOLD(J)
30            CONTINUE
              DO 40 J = I + 1, NUMEQN
                  SUM = SUM - A(I,J) * XOLD(J)
40            CONTINUE
              XNEW(I) = SUM / A(I,I)
50        CONTINUE

*     Check if termination condition satisfied
          I = 1
          DONE = (ABS(XOLD(I) - XNEW(I)) .LT. EPSIL)
*         While DONE and I < NUMEQN do:
60        IF (DONE .AND. I .LT. NUMEQN) THEN
              I = I + 1
              DONE = (ABS(XOLD(I) - XNEW(I)) .LT. EPSIL)
              GO TO 60
          END IF
          DONE = DONE .OR. (N .GT. NUMITS)

*     Display the new approximate solution
          PRINT *
          PRINT *, 'ITERATION #:', N
          DO 70 I = 1, NUMEQN
              PRINT *, 'X(', I, ') = ', XNEW(I)
70        CONTINUE

*     Copy XOLD to XNEW to prepare for next iteration
          DO 80 I = 1, NUMEQN
              XOLD(I) = XNEW(I)
80        CONTINUE

          GO TO 20
      END IF

      END  

