Solving a linear system with Gauss-Seidel's method

4 -1 0 -1 0 0
-1 4 -1 0 -1 0
0 -1 4 0 0 -1
-1 0 0 4 -1 0
0 -1 0 -1 4 -1
0 0 -1 0 -1 4

1
1
1
1
1
1

1
1
1
1
1
1

0.57024383544921875     
  0.71393871307373047     
  0.57024383544921875     
  0.57118320465087891     
  0.71261024475097656     
  0.57118320465087891

0.27777910232543945     
  0.62212115526199341     
  0.54369938373565674     
  0.48600935935974121     
  0.66520172357559204     
  0.55222527682781219

program hope
implicit none

integer :: p, n, n_max
real :: tol
98 write(*,*)"In order to solve the linear system, do you want to use Jacobi(1)'s method or Gauss-Seidel(2)'s one?"
read(*,*)p
write(*,*)"Enter the order n of the matrix A, a tolerance and a maximum number of iterations"
read(*,*)n, tol, n_max
if (p==1) then
call jacobi(n_max,n,tol)
else if (p==2) then
call gauss(n_max,n,tol)
else 
write(*,*)"Please choose a valid option"
goto 98
end if
end program


subroutine jacobi (n_max,n,tol)
implicit none
integer                             :: i,j,k
real                                :: sup_sum, low_sum, norm
integer, intent(in)                 :: n_max, n
real, dimension(n)                  :: b, u,x
real, intent(in)                    :: tol
real, dimension(n,n)                :: A  



open(unit=10,file="A_matrix.dat")
open(unit=15,file="b_vector.dat")
open(unit=20,file="initial_vector.dat")
do i=1,n
  read(10,*) (A(i,j),j=1,n)
  read(15,*) b(i)
  read(20,*) x(i)
end do
close(10) 
close(15) 
close(20)
k=1
norm=100
do while ( (k<=n_max).and.(tol<norm) )

  sup_sum = b(1)
  do j=2,n
    sup_sum = sup_sum - A(1,j)*x(j)
  end do
  u(1)=sup_sum/a(1,1)
  
  do i=2,n
    low_sum=dot_product( A(i,1:i-1), x(1:i-1))
    sup_sum=dot_product( A(i,i+1:n), x(i+1:n))
    u(i)=(b(i)-low_sum-sup_sum)/A(i,i)   
  end do
  
  
  norm=0
  do i=1,n
  norm=norm+abs(x(i)-u(i))
  end do
  
  do i=1,n
    x(i)=u(i)
  end do
  write(*,*) k, norm
  k=k+1
end do

open(unit=25,file="final_vectorjacobi.dat")
do i=1,n
  write(25,*) x(i)
end do

close(25)


end subroutine


subroutine gauss (n_max,n,tol)
implicit none
integer                             :: i, j, k
real                                :: low_sum, sup_sum, norm
real, dimension(n,n)                :: A
real, dimension(n)                  :: b, x, oldx
integer, intent(in)                 :: n_max
real, intent(in)                    :: tol
integer, intent(in)                 :: n
  


open(unit=10,file="A_matrix.dat")
open(unit=15,file="b_vector.dat")
open(unit=20,file="initial_vector.dat")

do i=1,n
  read(10,*) (A(i,j),j=1,n)
  read(15,*) b(i)
  read(20,*) x(i)
end do
close(10) 
close(15) 
close(20)

k=1
norm=100
x=0
do while ( (k<=n_max).and.(tol<norm) )
oldx=x

  do j=2,n
    sup_sum = sup_sum - A(1,j)*x(j)
  end do
  x(1)=sup_sum/a(1,1) 

  do i=2,n
    low_sum=dot_product( A(i,1:i-1), x(1:i-1))
    sup_sum=dot_product( A(i,i+1:n), x(i+1:n))
    x(i)=(b(i)-low_sum-sup_sum)/A(i,i)   
  end do
  
  norm=0

  
  norm=sum(abs(x-oldx))
    write(*,*) k, norm
    
k=k+1
end do
open(unit=25,file="final_vectorgauss.dat")
do i=1,n
  write(25,*) x(i)
end do
close(25)
end subroutine