program sistema
  implicit none

  integer :: scelta,k
  real(8),dimension(3) :: b,x
  real(8),dimension(3,3) :: a

  do 
     write(*,*) "Seleziona l'esempio (1-4, 0 per uscire) "
     read(*,*) scelta

     select case (scelta)
     case(0)
        exit
     case(1)
        ! Esempio in cui GS converge ma non J (pag. 244)
        b=(/ 7,13,-4 /)
        a(1,:)=(/ 3,0,4 /)
        a(2,:)=(/ 7,4,2 /)
        a(3,:)=(/ -1,-1,-2 /)

     case(2)
        ! Esempio in cui J converge ma non GS (pag. 245)
        b=(/ -6,-5,3 /)
        a(1,:)=(/ -3,3,-6 /)
        a(2,:)=(/ -4,7,-8 /)
        a(3,:)=(/ 5,7,-9 /)

     case(3)
        ! Esempio in cui J converge ma non GS (pag. 245)
        b=(/ -6,-5,3 /)
        a(1,:)=(/ -3,3,-6 /)
        a(2,:)=(/ -4,7,-8 /)
        a(3,:)=(/ 5,7,-9 /)

     case(4)
        ! Esempio in cui J converge meglio di GS (pag. 246)
        b=(/ 22,5,-2 /)
        a(1,:)=(/ 7,6,9 /)
        a(2,:)=(/ 4,5,-4 /)
        a(3,:)=(/ -7,-3,8 /)
     end select

     write(*,*) "matrice a"
     do k=1,3
        write(*,*) a(k,1:3)
     end do
     write(*,*) "b = ",b
     write(*,*) " "

     call jacobi(3,a,b,x)
     write (*,*) "La soluzione e` ",x
     call gs(3,a,b,x)

     write(*,*) "La soluzione e` ",x
  end do

end program sistema

subroutine jacobi (n,a,b,x)
  implicit none

  integer :: n,k,i,j,maxiter=1000
  real(8),dimension(n,n) :: a
  real(8),dimension(n)   :: b,x,xnew,xold
  real(8) :: s



  xold=0.d0
  do k=1,maxiter
     if (k==maxiter) then
        write(*,'(a)') "Il metodo non converge"
        exit
     end if
     do i=1,n ! Questo ciclo calcola x^{(k+1)}
        s=0.d0
        do j=1,i-1
           s=s+a(i,j)*xold(j)
        end do
        do j=i+1,n
           s=s+a(i,j)*xold(j)
        end do
        xnew(i)=(b(i)-s)/a(i,i)
     end do

     if (maxval(abs(xnew-xold))<1.d-5) exit
     xold=xnew
  end do
  
  x=xnew
  write(*,*) "Iterazioni del metodo di Jacobi       = ",k
end subroutine jacobi

subroutine gs (n,a,b,x)
  implicit none

  integer :: n,k,i,j,maxiter=1000
  real(8),dimension(n,n) :: a
  real(8),dimension(n)   :: b,x,xnew,xold
  real(8) :: s



  xold=2.d0
  do k=1,maxiter
     if (k==maxiter) then
        write(*,'(a)') "Il metodo non converge"
        exit
     end if
     do i=1,n ! Questo ciclo calcola x^{(k+1)}
        s=0.d0
        do j=1,i-1
           s=s+a(i,j)*xnew(j)
        end do
        do j=i+1,n
           s=s+a(i,j)*xold(j)
        end do
        xnew(i)=(b(i)-s)/a(i,i)
     end do

     if (maxval(abs(xnew-xold))<1.d-12) exit
     xold=xnew
  end do
  
  x=xnew
  write(*,*) "Iterazioni del metodo di Gauss-Seidel       = ",k
end subroutine gs
