gaussseidel <- function ( A, b, tol = 10e-7, maxiter = 100 )

#*****************************************************************************80
#
## gaussseidel uses Gauss Seidel iteration to solve a linear system.
#
#  Licensing:
#
#    Copyright 2016 James P. Howard, II
#
#    The computer code and data files on this web page are distributed under
#    https://opensource.org/licenses/BSD-2-Clause, the BSD-2-Clause license.
#
#  Modified:
#
#    18 February 2020
#
#  Author:
#
#    Original R code by James Howard;
#    Modifications by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017,
#    ISBN13: 978-1-4987-2363-3.
#
{
  n <- length ( b )
  iter <- 0

  L <- A
  L[upper.tri(A, diag = FALSE)] <- 0
  U <- A
  U[lower.tri(A, diag = TRUE)] <- 0
  Linv <- solve ( L )

  x <- rep ( 0.0, n )
  newx <- rep ( tol * 10, n )

  while ( tol < vecnorm ( newx - x ) )
  {
    if ( maxiter < iter )
    {
      warning ( "gaussseidel: Iteration maximum exceeded." )
      break
    }
    x <- newx
    newx <- Linv %*% ( b - U %*% x )
    iter <- iter + 1
  }

  return ( as.vector ( newx ) )
}