program main
!---------------------------------------------------------------------------------
! Programma principale per la chiamata della subroutine Jacobi.                  |
! Il programma main legge una matrice A, il termine noto b, tolleranza di        |
! uscita e numero massimo di iterazioni. Risolve il sistema "Ax = b" chiamando   |
! la subroutine Jacobi.                                                          |
!---------------------------------------------------------------------------------
implicit none

! Dimensione massima
integer, parameter :: nmax = 100

! Dichiarazione delle variabili
integer            :: n,nitmax
integer            :: i,j
real(kind=8)       :: A(nmax,nmax),b(nmax),xnew(nmax),xold(nmax),s(nmax)
real(kind=8)       :: toll

! Apertura dell'unita' di Input
open(1,file='jacobi.dat',status='old')
! Leggi numero massimo di iterazioni e tolleranza
read(1,*) nitmax,toll
! Leggi dimensione del sistema
read(1,*) n
! Leggi la matrice A
do i = 1,n
   read(1,*) (A(i,j),j=1,n)
enddo
! Leggi il termine noto b
read(1,*) (b(i),i=1,n)
! Chiudi l'unita' di Input
close(1)

! Apri l'unita' di Output
open(2,file='jacobi.ris')
! Stampa numero massimo iterazioni e tolleranza di uscita
write(2,*) 'nitmax, toll:',nitmax,toll
! Stampa dimensione del sistema
write(2,*) 'Dimensione:',n
! Stampa A
write(2,*) 'A:'
do i = 1,n
   write(2,'(100e15.6)') (A(i,j),j=1,n)
enddo
! Stampa b
write(2,*) 'b:'
write(2,'(e15.6)') (b(i),i=1,n)

! Risolvi il sistema con il metodo di Jacobi
call jacobi(n,nmax,A,b,nitmax,toll,xnew,xold,s)

! Stampa la soluzione
write(2,*) 'x:'
write(2,'(e15.6)') (xnew(i),i=1,n)

! Chiudi l'unita' di Output
close(2)

end program main

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

subroutine jacobi(n,nmax,A,b,nitmax,toll,xnew,xold,s)
!---------------------------------------------------------------------------
! Risolve un sistema lineare con il metodo di Jacobi chiamando la function
! EUCL per il calcolo della norma euclidea di un vettore
!---------------------------------------------------------------------------
implicit none

! Dichiarazione delle variabili in Input
........
! Dichiarazione delle variabili in Output
........
! Dichiarazione delle variabili locali
........
! Dichiarazione della function
........

! Inizializzazione
iter = 0
scarto = 2.d0*toll
do i = 1,n
   xold(i) = 0.d0
enddo

! Ciclo di Jacobi
ciclo_princ: do while ( (scarto .ge. toll) .and. (iter .le. nitmax) )
   ! Incremento l'indice di iterazione
   iter = iter + 1
   ! Calcolo la nuova approssimazione componente "i" per componente
   do i = 1,n
      ! Calcolo SOM1
      ......
      ! Calcolo SOM2	  
      ......
	  ! Calcolo xnew
	  xnew(i) = .....
	  ! Calcolo s
	  s(i) = .....
   enddo
   ! Calcolo la norma di s
   scarto = EUCL(....)
   ! Aggiorno la soluzione xold prima di passare alla prossima iterazione
   do i = 1,n
      xold(i) = xnew(i)
   enddo
enddo ciclo_princ

end subroutine jacobi

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

function EUCL(....) result(norma)
!---------------------------------------------------------------------------
! Function per il calcolo della norma euclidea di un vettore.
!---------------------------------------------------------------------------
implicit none
! Dichiarazione delle variabili in Input
........
! Dichiarazione del risultato
........
! Dichiarazione delle variabili locali
........

! Calcolo della norma
........

end function EUCL
