C     PROGRAM GAUSE 
C
C     GAUSS-SEIDEL METHOD
C
C      OPEN(UNIT=6,FILE='GAUSE.RES')
C      OPEN(UNIT=6,FILE='LPT1')
      PARAMETER(ND=4)
      DOUBLE PRECISION EPS,A,B,C,INVQ,X,Y,E
      INTEGER IPVT
      DIMENSION A(ND,ND),B(ND),C(ND,ND),INVQ(ND,ND),X(ND),Y(ND),IPVT(ND)
      OPEN(UNIT=5,FILE='GAUSE.DAT')
      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='
    5 FORMAT(1H+,F6.1,3X)
      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,C,INVQ,NIT,EPS,X,KIT,E,ITEST,Y,IPVT)
      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,C,INVQ,NIT,EPS,X,KIT,E,ITEST,Y,IPVT)
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     C: GAUSS-SEIDEL METHOD ITERATIVE MATRIX
C     INVQ: ΚΑΤΩ ΤΡΙΓΩΝΙΚΟ ΚΟΜΜΑΤΙ ΤΟΥ ΠΙΝΑΚΑ Α
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
      INTEGER IPVT
      DOUBLE PRECISION EPS,A,B,C,INVQ,X,Y,CN,S,E
      DIMENSION A(ND,N),C(ND,N),INVQ(ND,N),B(ND),X(ND),Y(ND),IPVT(ND)
      ITEST=1
C
C     NORM OF THE GAUSS-SEIDEL ITERATION MATRIX
C     
      
C     ARXIKOPOIHSH TOY Q (XRHSIMOPOIV TON C VS ENDIAMESO PINAKA) 
      DO 3 J=1,N
        DO 1 I=1,J-1
    1     C(I,J)=0
        DO 2 I=J,N
    2     C(I,J)=A(I,J)
    3 CONTINUE    
   
C     YΠΟΛΟΓΙΣΜΟΣ ΤΟΥ INV(Q)
      CALL INVERT(ND,N,C,INVQ,CN,IPVT,Y)

C     ΠΟΛΛΑΠΛΑΣΙΑΣΜΟΣ ΤΟΥ (Q^-1)*P ΜΕ ΤΑΥΤΟΧΡΟΝΟ ΜΗΔΕΝΙΣΜΟ ΤΟΥ C
      
       DO 5 I=1,NNN
        DO 6 J=1,NN
          C(I,J)=0D0
          DO 7 K=1,NN-1
    7        C(I,J)=C(I,J)-INVQ(I,K)*A(K,J)
    6   CONTINUE
    5 CONTINUE
      
C     ΥΠΟΛΟΓΙΣΜΟΣ ΤΗΣ ΝΟΡΜΑΣ ΓΡΑΜΜΗΣ ΤΟΥ C
      CN=0.D0
      DO 10 I=1,N
      S=0.D0
         DO 20 J=1,N
            S=S+ABS(C(I,J))
   20    CONTINUE
      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,45) (X(I),I=1,N)
   45 FORMAT(' X=',4(D15.10,1X))
C
C     ERROR TEST
C
      E=E*CN/(1.D0-CN)
      IF(E.LE.EPS) RETURN
   30 CONTINUE
      ITEST=2
      RETURN
      END        
      

C
C
C
      SUBROUTINE INVERT (NDIM,N,A,X,COND,IPVT,WORK)

C     ΔΗΛΩΣΕΙΣ ΜΕΤΑΒΛΗΤΩΝ       
      INTEGER  I,J,N,NDIM
      DOUBLE PRECISION A(NDIM,N),X(NDIM,N),WORK(N),COND
      INTEGER  IPVT(N)

C     ΓΙΝΕΤΑΙ Ο ΠΙΝΑΚΑΣ Χ ΙΣΟΣ ΜΕ ΤΟΝ ΜΟΝΑΣΙΑΙΟ Ι(ΝxΝ)
      DO 20 J=1,N
       DO 10 I=1,N
        X(I,J)=0
   10  CONTINUE
      X(J,J)=1
   20 CONTINUE

C     Τρίγωνοποιείται ο Α και υπολογίζεται ο δείκτης κατάστασής του (COND)
C     μέσω της υπορουτίνας DECOMP
      CALL DECOMP (NDIM,N,A,COND,IPVT,WORK)

C     Εάν ο πίνακας A εχεί άσχημο δείκτη κατάστασης τοτε βγαζεί μήνυμα
C     στην ο8όνη και γυρίζει στο κυρίως πρόγραμμα.
      IF((COND+1).EQ.COND) RETURN

C     Εδώ καλείται η SOLVE μία φορά για κάθε στήλη του Χ 
      DO 40 I=1,N
   40 CALL SOLVE(NDIM,N,A,X(1,I),IPVT)

C     Εχω τον αντίστροφο του Α στον Χ και επιστρέφω στο κυρίως πρόγραμμα
      RETURN
      END