program gauss_seidel_iterative
      implicit double precision (a-h,o-z)
      parameter (kdim=20)
      dimension a(kdim,kdim),b(kdim),x(kdim)
c
      open(1,file='matrix.dat')
      read(1,*)n
      do i=1,n
        read(1,*) (a(i,j),j=1,n),b(i)
        x(i)=0.d0  !  initialization
      enddo
      close(1)
      write(*,*)' Enter Number of Maximum Iterations '
      read(*,*)niter
      write(*,*)' Enter Relaxation Factor'
      read(*,*)omega
c
      do 100 iter=1,niter
        do i=1,n
          sum=0.d0
          do j=1,n
          if(i.ne.j) sum=sum+a(i,j)*x(j)
          enddo
          x(i)=(b(i)-sum)/a(i,i)*omega+x(i)*(1.d0-omega)
        enddo
c-compute residual, for the convergence criterion (if any...)
c-attention: this computation introduces extra CPU cost
        residual=0.d0
        do i=1,n
          sum=0.d0
          do j=1,n
          sum=sum+a(i,j)*x(j)
          enddo
          residual=residual+(sum-b(i))**2
        enddo
      write(*,*)' Iteration ',iter,'  Residual ',dsqrt(residual)
 100  continue
c
      open(1,file='solution.dat')
      do i=1,n
      write(1,*)i,x(i)
      enddo
      close(1)
c
      stop
      end