ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  MAIN: Jacobi Iteration
c 
c  N: number of equations 
c  ALOT: just a large parameter
c  A,x,f: Ax=f arrays
c  old* : values from previous iteration
c  xtrue: the actual solution... from Gaussian Elim. algorithm
c  E :  error of x guess from actual solution
c  del : error of current x from previous x   
c  *ratio: ratio of errors
c  epsil : small parameter to test for convergence
c  i,j,k : iterator integers
c  sumAx : gives a sum necessary to calculate x
c  mag* : temp variable
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      PROGRAM MAIN

      implicit none

      INTEGER N,ALOT
      PARAMETER(N=4)
      PARAMETER(ALOT=30)
      INTEGER i,j,k
      REAL A(N,N),x(N,ALOT),f(N),xtrue(N),magdel,magE,epsil
      REAL del(ALOT),olddel,delratio(ALOT),olddelratio
      REAL E(ALOT),oldE,Eratio(ALOT),sumAx
 
      epsil = 1.E-3
cccc ignore next seven initialization steps
      del(1) = 1.5
      olddel = 2.
      delratio(1) =1.5
      delratio(2) = 2.23
      olddelratio = 1.21
      Eratio(1) = 2.78 
      oldE = 1.1
      

cccc assigning values to A,f matrices
      CALL GETAfx(A,f,xtrue,N)
      
cccc guessing the initial values of x
      DO 10 i=1,N
         x(i,1) = 0.0
 10   CONTINUE


cccc Jacobi iteration
      DO 20 k=1,ALOT
            
cccc new x are calculated
         DO 30 i=1,N
            sumAx = 0.0
            CALL SUM1(A,x,i,k,N,ALOT,sumAx)
            x(i,k+1) = x(i,k) + 1./A(i,i) * (f(i) - sumAx)
 30      CONTINUE

cccc getting delta(k), E(k)
         magE = 0.0
         DO 40 i=1,N
            magE = magE + ( x(i,k+1) - xtrue(i) )**2
 40      CONTINUE
         E(k) = sqrt(magE)

         IF(k.GT.1) THEN
            magdel = 0.0
            DO 45 i=1,N
               magdel = magdel + ( x(i,k+1) - x(i,k) )**2
 45         CONTINUE
            del(k) = sqrt(magdel) 
         END IF

cccc getting delta, E ratios
         IF(k.GT.1) THEN
            Eratio(k) = E(k)/oldE

            IF(k.GT.2) THEN
               delratio(k) = del(k)/olddel
cccc test for convergence
               IF(abs(E(k)-oldE).LE.epsil) THEN
                  PRINT*,x(1,k)
                  GOTO 100
               END IF
            ENDIF

         END IF
c         PRINT*,del(k),E(k)

cccc update variables
         olddel = del(k)
         olddelratio = delratio(k) 
         oldE = E(k)

         PRINT*,x(1,k)

 20   CONTINUE

c!c!c! end Jacobi iteration
 100  print*,k

cccc printing results to file
      CALL PRINTVALS(A,x,f,del,delratio,E,Eratio,k,N,ALOT) 


      END       
ccccc MAIN ccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


cccccccccccccccccccccccccccccccccccc
c
c Function SUM: 
c sums A(i,j)*oldx(j) from j=1->N
c
cccccccccccccccccccccccccccccccccccc

      SUBROUTINE SUM1(A,x,i,k,N,ALOT,sumAx)

         implicit none
         INTEGER i,j,k,N,ALOT
         REAL A(N,N),x(N,ALOT),sumAx

         DO 200 j=1,N 
            sumAx = sumAx + A(i,j)*x(j,k)
 200     CONTINUE

         RETURN
         END
cccccccc SUM cccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccc      

cccc initialize A,xtrue,f
         SUBROUTINE GETAfx(A,f,xtrue,N)
cccc initialize A,xtrue,f

         implicit none
         INTEGER N
         REAL A(N,N),f(N),xtrue(N)

         A(1,1) = -4.
         A(1,2) = 1.
         A(1,3) = 1.
         A(1,4) = 1.

         A(2,1) = 1.
         A(2,2) = -4.
         A(2,3) = 1.
         A(2,4) = 1.

         A(3,1) = 1.
         A(3,2) = 1.
         A(3,3) = -4.
         A(3,4) = 1.

         A(4,1) = 1.
         A(4,2) = 1.
         A(4,3) = 1.
         A(4,4) = -4.

         f(1) = 1.
         f(2) = 1.
         f(3) = 1.
         f(4) = 1.

cccc from gaussian elimination algorithm
         xtrue(1) = -1.
         xtrue(2) = -1.
         xtrue(3) = -1.
         xtrue(4) = -1.
      
      RETURN
      END
ccccc GETAfx cccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccc


cccc output
      SUBROUTINE PRINTVALS(A,x,f,del,delratio,E,Eratio,k,N,ALOT)
cccc output

      implicit none

      INTEGER k,i,j,N,ALOT
      REAL A(N,N),x(N,ALOT),f(N),del(ALOT),delratio(ALOT)
      REAL E(ALOT),Eratio(ALOT)

      OPEN(UNIT=99,FILE='hw4_data',STATUS='OLD')
      WRITE(99,*) 'Jacobi Iteration'
      WRITE(99,*) ' '
      WRITE(99,*)'MATRIX A = '

      DO 300, i=1,N
         WRITE(99,'(4(F5.2,1x))') ( A(i,j),j=1,N)
 300  CONTINUE

      WRITE(99,*) ' '
      WRITE(99,*) 'MATRIX f = ' 
      DO 305, i=1,N
         WRITE(99,'(F8.2)') f(i)
 305  CONTINUE
      WRITE(99,*) ' '

cccc k
cccc x,delta,deltaratio,E,Eratio
      WRITE(99,*) 'k      x1     x2      x3      x4'

      DO 310 j=1,k
         WRITE(99,311) j,(x(i,j),i=1,4)
 310  CONTINUE   
 311  FORMAT(I2,1X,(4(F8.4)))

      WRITE(99,*) ' '
      WRITE(99,*) 'k      delta    delratio    E   Eratio  '

      DO 315 i=1,k
         WRITE(99,399) i,del(i),delratio(i),E(i),Eratio(i) 
 315  CONTINUE
 399  FORMAT(I2,1X,(4(F10.5)) )

      CLOSE(99)
      
      
      RETURN
      END
