# Solving a linear system with Gauss-Seidel's method

1. Jun 27, 2011

### fluidistic

Edit: Problem solved. Look at the bottom of this post.
Ok so far I've made a working program that can solve for x in Ax=b using Jacobi method.
I've tried to use Gauss-Seidel method too, but I failed in doing so.
I took as an example, the following A matrix:
Code (Text):
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

with the following b vector:
Code (Text):
1
1
1
1
1
1
I use the initial "guess" or vector:
Code (Text):
1
1
1
1
1
1
.
As a result, I should get a value for the vector x close to
Code (Text):
0.57024383544921875
0.71393871307373047
0.57024383544921875
0.57118320465087891
0.71261024475097656
0.57118320465087891
but I get the values
Code (Text):
0.27777910232543945
0.62212115526199341
0.54369938373565674
0.48600935935974121
0.66520172357559204
0.55222527682781219
Here comes my Jacobi and then Gauss-Seidel codes inside the same program:
Code (Text):
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?"
write(*,*)"Enter the order n of the matrix A, a tolerance and a maximum number of iterations"
if (p==1) then
call jacobi(n_max,n,tol)
else if (p==2) then
call gauss(n_max,n,tol)
else
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
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
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
I must say I'm a bit confused in Gauss-Seidel subroutine when I treat x and oldx as arrays rather than real numbers. I know I do get wrong values for x(i) for Gauss-Seidel method... but I don't see the error.
I'd appreciate any help to search the error(s). Thanking you in advance.

EDIT: Sorry guys! I think I've found out the error by a pure lucky guess. I changed the 2 starting values for the loops from 1 to n instead of from 2 to n. Well I won't erase this thread, so that everyone can use this code :)