***GAUSS-SEIDEL METHOD*************************************************
*         PROGRAM FOR SOLUTION OF SYSTEM OF LINEAR EQUATION BY        *
*                    GAUSS-SEIDEL ITERATION METHOD                    *
***********************************************************************
*
      program GAUSS_SEIDEL
      implicit doubleprecision(a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd),X(nd),C(nd)
      integer p
*
      print*, 'Enter number of equation, n:'
      read(*,*) n
      m = n + 1
      write(6,*)
      print*, 'Enter the augmented coefficient matrix row-wise:'
      read(*,*) ((A(i,j), j = 1,m), i = 1,n)
      write(6,*)
      print*, 'Maximum number of iterations:'
      read(*,*) itermax
      write(6,*)
      print*, 'Acceptable error (maximum norm of residual):'
      read(*,*) epsil
      write(6,*)'-------------------------------------------------------
     $---------------------'
c
c     checking for diagonal dominance
      p = 0
      do i = 1, n
         summ = 0
         do j = 1, n
            if (i.ne.j) then
               summ = summ + dabs(A(i,j))
            endif
         enddo
         if (dabs(A(i,i)) < summ) then 
            p = 1                           
         endif     
      enddo
*             
      if (p == 1) then
         write(6,*) 'The coefficient matrix is not diagonally dominant.
     $Hence convergence may not be attained.'
      elseif(p == 0) then
         write(6,*) 'The coefficient matrix is diagonally dominant. Henc
     $e the solution must converge'
      endif
c
c     extraction of right-hand-side column vector
      do i = 1, n
         C(i) = A(i,m)
      enddo
c
c     calculation of maximum norm of rhs vector
      rhsnorm = 0
      do i = 1, n
         rhsnorm = max(rhsnorm, dabs(C(i)))
      enddo
      rhsnorm = rhsnorm + 1d-30
c
c     initial guess of the solution
      do i = 1, n
         X(i) = 0
      enddo
*
      iter = 0
      relresidual = 1d9
10    if ((relresidual > epsil).and.(iter < itermax)) then
         do i = 1, n
            sum = 0
            do j = 1, n
               if (i.ne.j) then
                  sum = sum + A(i,j)*X(j)
               endif
            enddo
            X(i) = (C(i) - sum) /A(i,i)
         enddo
c
c        compute maximum norm of the residual vector
         residnorm = 0
         do i = 1, n
            resid = C(i)
            do j = 1, n
               resid = resid - A(i,j)*X(j)
            enddo
            residnorm = max(residnorm, dabs(resid))
         enddo
         relresidual = residnorm /rhsnorm      ! relative residual norm
         iter = iter + 1
         go to 10
      endif 
*
      if (relresidual > epsil) then
         write(6,*) '-------------------------------------------------'
         write(6,*)
         write(6,*) 'The convergence not attained. Try increasing the ma 
     $ximum number of iterations.'
         write(6,*)
         write(6,*) '-------------------------------------------------'
      else
         write(6,*) '              -----------------------------------'
         write(6,*) '                 The converged solution vector is:'
         write(6,*)
         write(6,101) (X(i),i = 1,n)
101      format('',29x,f12.4)
         write(6,*) '--------------------------------------------------'
         write(6,*) 'The number of iteration performed before convergenc
     $e =', iter
         write(6,*)'The maximum norm of relative residual =',relresidual
         write(6,*) '--------------------------------------------------'
      endif
*
      stop
      end