C     PROGRAM GAUSE 
C
C     GAUSS-SEIDEL METHOD
C
C      OPEN(UNIT=5,FILE='GAUSE.DAT')
C      OPEN(UNIT=6,FILE='GAUSE.RES')
C      OPEN(UNIT=6,FILE='LPT1')
      PARAMETER(ND=10)
      DOUBLE PRECISION EPS,A,B,X,Y,E
      DIMENSION A(ND,ND),B(ND),X(ND),Y(ND)
      READ(5,*) N,NIT,EPS
      READ(5,*) ((A(I,J),J=1,N),I=1,N)
      READ(5,*) (B(I),I=1,N)
      READ(5,*) (X(I),I=1,N)
      WRITE(6,*) 'INPUT:'
      WRITE(6,*) 'N=',N,' NIT=',NIT,' EPS=',EPS
      WRITE(6,*) 'A='
      DO 10 I=1,N
      WRITE(6,*) (A(I,J),J=1,N)
   10 WRITE(6,*) ' '
      WRITE(6,*) 'B=',(B(I),I=1,N)
      WRITE(6,*) 'X=',(X(I),I=1,N)
      WRITE(6,*) ' '
      WRITE(6,*) 'OUTPUT:'
      CALL GAUSEID(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y)
      WRITE(6,*) 'ITEST=',ITEST,' KIT=',KIT
      WRITE(6,*) 'X=',(X(I),I=1,N)
      WRITE(6,*) 'E=',E
      STOP
      END
C-----------------------------------------------------------------
      SUBROUTINE GAUSEID(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y)
C
C     SOLUTION OF THE LINEAR SYSTEM AX=B BY THE GAUSS-SEIDEL METHOD
C
C     INPUT:
C     ND: MAXIMUM COLUMN DIMENSION OF ARRAYS
C     N: DIMENSION OF ARRAYS
C     A: SYSTEM MATRIX (NXN)
C     B: RIGHT-HAND SIDE VECTOR (N) 
C     NIT: MAXIMUM NUMBER OF ITERATIONS
C     EPS: CONVERGENCE TEST NUMBER
C
C     OUTPUT:
C     X: APPROXIMATE SOLUTION VECTOR (N)
C     KIT: NUMBER OF EXECUTED ITERATIONS
C     E: ERROR BOUND
C     ITEST = 1, CONVERGENCE TEST SATISFIED
C           = 2, CONVERGENCE TEST NOT SATISFIED
C
      DOUBLE PRECISION EPS,A,B,X,Y,CN,S,E
      DIMENSION A(ND,N),B(ND),X(ND),Y(ND)
      ITEST=1
C
C     NORM OF THE JACOBI ITERATION MATRIX
C     
      CN=0.D0
      DO 10 I=1,N
      S=0.D0
         DO 20 J=1,N
         IF(J.EQ.I) GO TO 20
         S=S+ABS(A(I,J))
   20    CONTINUE
      S=S/DABS(A(I,I))
      IF(S.GT.CN) CN=S
   10 CONTINUE
C
C     ITERATIONS
C
      DO 30 K=1,NIT
      KIT=K
         DO 25 I=1,N
   25    Y(I)=X(I)
      E=0.D0
         DO 40 I=1,N
         X(I)=B(I)
            DO 50 J=1,N
            IF(J.EQ.I) GO TO 50
            X(I)=X(I)-A(I,J)*X(J)
   50       CONTINUE
         X(I)=X(I)/A(I,I)
         S=DABS(X(I)-Y(I))
         IF(S.GT.E) E=S
   40    CONTINUE
      WRITE(6,*) '  X=',(X(I),I=1,N)
C
C     ERROR TEST
C
      E=E*CN/(1.D0-CN)
      IF(E.LE.EPS) RETURN
   30 CONTINUE
      ITEST=2
      RETURN
      END