!! File  : guss-seidel.f90
!!
!! Creation Date: 03.11.2015; Time: 14:48
!! Time-stamp: <2015-11-03 15:54:54 (andrzejbaran)> 
!!
!! Andrzej Baran
!! IF UMCS, ul. Radziszewskiego 10, 20-031 Lublin, PL
!! baran@kft.umcs.lublin.pl
!!
program test_gauss_seidel
  IMPLICIT NONE
  integer, parameter    :: n = 3
  integer, parameter    :: maxit =  20
  real, parameter       :: tol = 1e-4  !  tolerance
  real, dimension(n,n)  :: a
  real, dimension(n)    :: b
  real, dimension(n)    :: x

  data a /12., 1., 3., 3., 5., 7., -5., 3., 13./, &   ! kolumnowo
       b /1., 28., 76./
  ! Solution /1., 3., 4/

  ! solve ax=b
  call gauss_seidel(a, b, n, x, tol, maxit)

CONTAINS

  Subroutine gauss_seidel(a, b, n, x, tol, maxit)
    IMPLICIT NONE
    real, dimension(n,n)    :: a
    real, dimension(n)      :: b
    real, dimension(n)      :: x
    integer                 :: i, j, k, n, iter, maxit
    real                    :: old, sum, rel, s, aux, tol
    ! ================================================================
    ! The routine solves the system Ax = b by Gauss-Seidel iterations
    ! ================================================================
    ! INPUT:
    !    n - number of equations
    !    a - n x n matrix of coeff.
    !    b - n vector of the right side
    !    maxit - max number of iterations
    ! Output:
    !    x - solution (vector)
    ! ================================================================    
    ! check diagonal dominant
    do i = 1,n
       s = 0
       do j = 1,n
          if(i .ne. j) s = s + abs(a(i,j))
       end do
       if(abs(a(i,i)) .lt. s) then
          write(*,*) "The A matrix is not diagonally-dominant."
          stop
       end if
    end do
    do i = 1,n
       aux = a(i,i)
       do j = 1,n
          a(i,j) = a(i,j)/aux
       end do
       b(i) = b(i)/aux
    end do
    iter = 0
    k = 0
    do while (iter < maxit .and. k == 0)
       iter = iter+1
       print "(/, a, i4)", "Iteration:", iter
       k = 1
       do i = 1,n
          old = x(i)
          sum = b(i)
          do j = 1,n
             if(i .ne. j) then
                sum = sum - a(i,j)*x(j)
             end if
          end do
          x(i) = sum
          print *, x(i)
          if(x(i) .ne. 0 .and. k==1) rel = (abs(x(i)-old)/x(i))*100 
          if(rel .gt. tol) k = 0
       end do
    end do
  end Subroutine gauss_seidel

END program test_gauss_seidel


