PROGRAM iterativi
  IMPLICIT NONE
  INTERFACE
     SUBROUTINE jacobi(nit,a,b,x)
       INTEGER :: nit
       REAL(KIND(0.d0)),DIMENSION(:,:)::a
       REAL(KIND(0.d0)),DIMENSION(:)::b,x
     END SUBROUTINE jacobi
     SUBROUTINE gauss_seidel(nit,a,b,x)
       INTEGER :: nit
       REAL(KIND(0.d0)),DIMENSION(:,:)::a
       REAL(KIND(0.d0)),DIMENSION(:)::b,x
     END SUBROUTINE gauss_seidel
  END INTERFACE

  REAL(KIND(0.d0)),DIMENSION(:,:),ALLOCATABLE::a
  REAL(KIND(0.d0)),DIMENSION(:),ALLOCATABLE::b,x
  INTEGER :: n,i,j,nit

  WRITE(*,*)"n="
  READ(*,*)n
  ALLOCATE (a(n,n),b(n),x(n))
  a=1
!  a(n,1:n-1)=2
!  a(n-1,1:n-2)=2
  DO i=1,n
     a(i,i)=n+i
  END DO


  DO i=1,n
     b(i)=0
     DO j=1,n
        b(i)=b(i)+a(i,j)
     END DO
  END DO

  x=0  
  WRITE(*,*)"numero di iterazioni"
  READ(*,*)nit
  CALL jacobi(nit,a,b,x)
  WRITE(*,*)"errore nel metodo di Jacobi",MAXVAL(ABS(x-1))
  x=0
  CALL gauss_seidel(nit,a,b,x)
  WRITE(*,*)"errore nel metodo di Gauss-Seidel",MAXVAL(ABS(x-1))
END PROGRAM iterativi

SUBROUTINE jacobi(nit,a,b,x)
  ! metodo di Jacobi applicato al sistema Ax=b
  ! si applicano nit iterazioni
  IMPLICIT NONE
  REAL(KIND(0.d0)),DIMENSION(:,:)::a
  REAL(KIND(0.d0)),DIMENSION(:)::b,x
  REAL(KIND(0.d0)),DIMENSION(:),ALLOCATABLE::y
  INTEGER :: n,i,j,k,nit

  n=SIZE(x)
  ALLOCATE(y(n))
  DO k=1,nit
     DO i=1,n
        y(i)=0
        DO j=1,n
           IF(i/=j) y(i)=y(i)+a(i,j)*x(j)
        END DO
        y(i)=(b(i)-y(i))/a(i,i)
     END DO
     x=y
  END DO
END SUBROUTINE jacobi

SUBROUTINE gauss_seidel(nit,a,b,x)
  ! metodo di Gauss-Seidel applicato al sistema Ax=b
  ! si applicano nit iterazioni
  IMPLICIT NONE
  REAL(KIND(0.d0)),DIMENSION(:,:)::a
  REAL(KIND(0.d0)),DIMENSION(:)::b,x

  INTEGER :: n,i,j,k,nit

  n=SIZE(x)

  DO k=1,nit
     DO i=1,n
        x(i)=0
        DO j=1,n
           IF(i/=j) x(i)=x(i)+a(i,j)*x(j)
        END DO
        x(i)=(b(i)-x(i))/a(i,i)
     END DO
  END DO
END SUBROUTINE gauss_seidel






