C     PROGRAM GAUSE 
C	  Iwannou Spyros,Kontopodhs Dhmhtrhs
C	  To compilation apaitei thn decsold (me double precision)
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,Q,IQ,C,CNORM,P
	  INTEGER I,J,K
      DIMENSION A(ND,N),B(ND),X(ND),Y(ND),Q(ND,N),IQ(ND,N),C(ND,N)
	  DIMENSION P(ND,N)
      ITEST=1
C

C	  Eyresh tou Q (katw trigwnikou pou prokuptei apo ton A
	  DO 15 I=1,N
	  
	  DO 20 J=1,N

	  IF (I.GE.J) THEN
	  Q(I,J)=A(I,J)
	  ELSE
	  Q(I,J)=0
	  ENDIF
20	  CONTINUE
15	  CONTINUE

C	  Eyresh tou P (anw trigonikou, me stoixeia antitheta tou A)

	  DO 22 I=1,N
	  DO 23 J=1,N

      P(I,J)=(Q(I,J)-A(I,J))

23	  CONTINUE
22	  CONTINUE
	  CALL INVERT(Q,IQ,N,ND)


C	  Eyresh tou C=IQ*P (IQ=Antistrofos tou Q)

	  DO 73 I=1,N
	  DO 24 J=1,N
	  C(I,J)=0
	  DO 25 K=1,N
	  C(I,J)=C(I,J)+IQ(I,K)*P(K,I)
25	  CONTINUE
24	  CONTINUE
73	  CONTINUE
	  CALL NORM(N,C,CNORM)
C	  Thetoume gia logous omoiomorfias CN=CNORM

	  CN=CNORM






C     ITERATIONS
C
      DO 30 K=1,NIT
      KIT=K
         DO 28 I=1,N
   28    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        

        SUBROUTINE NORM(N,A,ANORM)
        DOUBLE PRECISION A(N,N), ANORM,T

        ANORM=0.0
        DO 65 J=1,N
         T=0.0
        DO 60 I=1,N
         T=T+ABS(A(I,J))
60      CONTINUE
        IF (T.GT.ANORM) ANORM=T
65      CONTINUE
        RETURN
        END

        SUBROUTINE INVERT(A,R,N,NDIM)
        INTEGER N,IPVT(N),I,J,NDIM
        DOUBLE PRECISION R(NDIM,N),A(NDIM,N),WORK(N),COND,T(N,N),M(N)
        DO 80 J=1,N
        DO 80 I=1,N
            T(I,J)=A(I,J)
80      CONTINUE
        CALL DECOMP(NDIM,N,A,COND,IPVT,WORK)
        IF ((COND+1).EQ.COND) THEN
        WRITE(6,*) 'INVERT:PINAKAS MH ANTISTREPSIMOS '
        RETURN
        ENDIF

        DO 82 J=1,N
        DO 81 I=1,N
        IF (I .EQ. J)   THEN
            M(I)=1
        ELSE
            M(I)=0
        ENDIF
81      CONTINUE
        CALL SOLVE(NDIM,N,A,M,IPVT)
        DO 83 I=1,N
        R(I,J)=M(I)
83      CONTINUE
82      CONTINUE
        DO 85 J=1,N
        DO 85 I=1,N

C       Epanaferei ton A pou prohgoumenws eixe krathsei ston pinaka T

        A(I,J)=T(I,J)
85      CONTINUE
        RETURN
        END