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 ) ) }