      program jacobi
C             programma per la soluzione della esercitazione sui
C            sistemi lineari con il metodo di Jacobi 
C
C
C   variabili utilizzate
C  
C   variabili intere
C
C   n     dimensione  dei vettori  e delle matrici (dato di input)
C   itmax  massimo numero di iterazioni per la soluzione del sistema lineare
C         (dato di input)
C
C   variabili reali
C
C   tol   tolleranza per il ciclo do while (dato di input)
C   scarto norma euclidea degli scarti
C
C   array reali
C
C   A(20,20)      matrice del sistema
C   b(20)       termine noto del sistema
C   xnew(20)    soluzione approssimata alla iterazione corrente x_{k+1}
C   xold(20)    soluzione approssimata alla iterazione precedente x_k
C   sc(20)      vettore scarto d_{k+1}=x_{k+1}-x_k
C
C
      implicit none
      integer i,j, itmax, iter, n
      real*8 A(20,20), xold(20),xnew(20), b(20),sc(20)
      real*8 scarto, tol

      open(15,file='jacobi.dat')
      open(16,file='jacobi.ris')

      read(15,*) n, itmax, tol
      write(16,*) 'dimensione n', n 
      write(16,*) 'numero max di iterazioni per il metodo', itmax
      write(16,*) 'tolleranza ',  tol
      write(16,*)  'matrice A' 
      do i=1,n 
         read(15,*) (A(i,j),j=1,n)
         write(16,1) (A(i,j),j=1,n)
      end do
      read(15,*) (b(i),i=1,n)
      write(16,*) 'termine noto'
      write(16,1) (b(i),i=1,n)
C
C  inizializzazione del vettore approssimazione della soluzione 
C  alla iterazione 0
C  x_0
      do i=1,n
        xold(i)=0.0
      end do
C
C   soluzione del sistema lineare nella subroutine jacobi
C
      call soljacobi(n,A,b,xold,xnew,tol,itmax,iter,scarto,sc)
C
C   stampa della soluzione
C
      if (iter.ge.itmax) then
      write(16,*)
      write(16,*) 'raggiunto il numero massimo di iterazioni'
      else
      write(16,*) 'effettuate ', iter, ' iterazioni'
      end if
      write(16,*)
      write(16,*) 'Soluzione approssimata mediante il metodo di Jacobi'
      write(16,1) (xnew(i),i=1,n)

      stop
1     format(5e15.6)
      end
     
C-----------------------------------------------------------------------

      subroutine  soljacobi(n,A,b,xold,xnew,tol,itmax,iter,
     1                scarto,sc)
      implicit none
      integer i, j, n,iter,itmax
      real*8 A(20,20), b(20), xold(20),xnew(20)
      real*8 sc(20), scarto, scarto_old, tol
      real*8 som1,som2,asint,R,euc
  
C
C  inizializzazione delle iterate e dello scarto
C
      iter=0
      scarto=2.00*tol
      scarto_old=1.0
      write(16,*) '        metodo di Jacobi             '
      write(16,*) '   iter     ||scarto||             M            R'
      write(16,*)
      do while (iter.lt.itmax.and.scarto.gt.tol)
         iter=iter+1
C
C  implementazione dell'iterazione di Jacobi
         do i=1,n
             som1=0.0
             som2=0.0
             do j=1,i-1
               som1=som1+A(i,j)*xold(j)
             end do
             do j=i+1,n
               som2=som2+A(i,j)*xold(j)
             end do
             xnew(i) = (b(i)-som1-som2)/A(i,i)
         end do
C
C   calcolo della norma euclidea dello scarto
C
         do i=1,n
           sc(i)=xnew(i) - xold(i)
         end do
         scarto=euc(n,sc)
C
C  stima del fattore di convergenza - M -  come rapporto delle norme
C  degli scarti consecutivi
C
         asint=scarto/scarto_old
C
C  velocita' asintotica di convergenza
C  per k (in questo caso iter) sufficientemente grande, la velocita' asintotica
C  di convergenza si puo' ottenere facendo il logaritmo del rapporto tra
C  gli scarti alla (k+1)-sima e alla k-sima iterazione
C
         R=-log10(asint)
C
C  stampa dei risultati relativi all'iterazione in corso
C

        write(16,100) iter,scarto,asint,r
 100    format(1x,i4,2x,e18.8,2f14.8)
C
C  aggiornamento degli scarti e della soluzione approssimata prima
C  di passare all'iterata successiva
C
        scarto_old=scarto
        write(18,*) 'iter' , iter
        do i=1,n
          xold(i)=xnew(i)
          write(18,*) xold(i)
        end do
C
C  fine del ciclo do while
      end do

       return
       end

C-----------------------------------------------------------------------
                    
      real*8 function euc(n,x) 
      implicit none
      integer nmax,n, i
      real*8 x(20)

      euc=0.0
      do i=1,n
        euc = euc + x(i)**2
      end do
      euc = dsqrt(euc)

      return 
      end

