C     MAIN PROGRAM FOR SOLVING AX=B USING GAUSS SEIDEL
      DIMENSION A(4,4),B(4), X(4)
      DATA N,MAX /4,35/
      DATA (A(1,J), J=1,4)/-40.0,28.5,0.0,0./
      DATA (A(2,J), J=1,4)/21.5,-60.0,28.5,0./
      DATA (A(3,J), J=1,4)/0.,21.5,-60.,28.5/
      DATA (A(4,J), J=1,4)/0.,0.,21.5,-40.0/
      DATA B(1),B(2),B(3),B(4)/-11.5,-10.,-10.,-18.5/
      PRINT * ,'A MATRIX AND VECTOR B'
      PRINT *
      DO I=1,N
       PRINT 3, (A (I,J), J=1,4)
        PRINT  4, B(I)
       END DO
      CALL SEIDEL(A,B,X,N,MAX,DIFF,ITER)
      ITER=ITER+1
      GO TO (50,51),ITER
50    PRINT 5, MAX,DIFF
      GO TO 52
51    PRINT 7,  MAX,DIFF
C     WRITE THE SOLUTION VECTOR
52    DO 53 I=1,N
53    PRINT 6, I, X(I)
2     FORMAT(1H1,/////,5X,'GAUSS SEIDEL SOLUTION OF A LINEAR SYSTEM OF
     *EQUATIONS',///,5X,'MATRIX A:'/)
3     FORMAT(1X, 9F12.6)
4     FORMAT (5X, F12.6)
5     FORMAT(//5X,'AFTER ',I3,3X, 'ITERATIONS MAX. ERROR =',/E15.6//)
6     FORMAT(5X,'X(',I2,')= ',F12.6)
7     FORMAT(//5X,'AFTER ',I3,'ITERATIONS THERE IS NO CONVERGENCE',
     */5X,'THE FINAL SOLUTION VECTOR IS:'/)
      STOP
      END
C     SUBROUTINE FOR SOLVING AX=B USING GAUSS-SEIDEL METHOD
      SUBROUTINE SEIDEL(A,B,X,N,MAX,DIFF,ITER)
      DIMENSION A(N,N),B(N),X(N)
      EPSI= 0.000001
      DO 50 I=1,N
50    X(I)=0.0
      NN=0
51    DIFF=0.0
      DO 56 I=1,N
      SUM=0.0
      DO 53 J=1,N
      IF(J-I) 52,53,52
52    SUM=SUM+A(I,J)*X(J)
53    CONTINUE
      XNEW=(-SUM+B(I))/A(I,I)
      CHECK=ABS((X(I)-XNEW)/XNEW)
      IF(CHECK-DIFF) 55,55,54
54    DIFF=CHECK
55    X(I)=XNEW
56    CONTINUE
      NN=NN+1
      IF(NN-MAX) 57,58,58
57    IF(DIFF-EPSI) 59,59,51
58    ITER=1
      MAX=NN
      RETURN
59    ITER=0
      MAX=NN
      RETURN
      END

