julia# Implementing matrix term version of Gauss-seidel

I am trying to implement the below description from Ch. 11 of Heath's "Scientific Computing An Introductory Survey" of the Gauss-Seidel iterative method for solving a system of linear equations

however, I must be missing something because writing the equation exactly as is fails to produce the correct results. Please help!

```
using LinearAlgebra: lu
# From dicretization of Laplace equation
A = [
4.0 -1.0 -1.0 0.0
-1.0 4.0 0.0 -1.0
-1.0 0.0 4.0 -1.0
0.0 -1.0 -1.0 4.0]
# block diagonal matrix of A
D = [
4.0 -1.0 0.0 0.0
-1.0 4.0 0.0 0.0
0.0 0.0 4.0 -1.0
0.0 0.0 -1.0 4.0]
b = [0, 0, 1, 1] # rhs for laplace w/ Dirichlet boundary conditions
L, U = lu(A)
# gauss seidel using matrix terms
niters = 100
D_plus_L_inv = inv(D + L)
x = zeros(length(b))
for k in 1:niters
x = D_plus_L_inv*(b - U*x)
end
# clearly not equal
@show A \ b
# A \ b = [0.12499999999999999, 0.12499999999999999, #0.37499999999999994, 0.37499999999999994]
@show x
# x = [0.021795262674411238, 0.023610129063946845, 0.14893710589872386, 0.14211027447080438]
```

Solution

Per @Ian Bush comments, I was misinterpreting what L, U, and D were in the equations. L and U are strict lower and upper triangular matrices of A, while D is the diagonal entries of A. Here is the corrected function matching the described equation (even though computing the inverse D + L may not be done in practice).

```
using LinearAlgebra
"""
gauss_seidel_matrix_form_solver(
A_::Matrix, b::Vector, x0::Vector, niters::Int)
Return solution to linear system `Ax = b` via matrix Gauss-Seidel method.
# References
[1] : Ch. 11.5.3 Heath.
"""
function gauss_seidel_matrix_form_solver(
A_::Matrix, b::Vector, x0::Vector, niters::Int)
A = copy(A_)
# Get diagonal entries of A, store in matrix
D_vec = diag(A)
D = diagm(D_vec)
# Get the strict lower and upper triangular matrices
L = LowerTriangular(A)
L[diagind(L)] .= 0
U = UpperTriangular(A)
# Compute the inverse of D + L ahead of time
D_plus_L_inv = inv(D + L)
# initialize the guess
x = x0
# gauss seidel iterations
for k in 1:niters
x = D_plus_L_inv*(b - U*x)
end
return x
end
A = [
4.0 -1.0 -1.0 0.0
-1.0 4.0 0.0 -1.0
-1.0 0.0 4.0 -1.0
0.0 -1.0 -1.0 4.0]
b = [0, 0, 1, 1]
u_gs = gauss_seidel_matrix_form_solver(A, b, zeros(length(b)), 100)
u_builtin = A \ b
@show u_gs
# u_gs = [0.125, 0.125, 0.375, 0.375]
@show u_builtin
# u_builtin = [0.12499999999999999, 0.12499999999999999, 0.37499999999999994, 0.37499999999999994]
```

- How to get equally scaled axes with Plots in Julia
- Scatter plot of two rows of a DataFrame in Julia using Plotly
- Adding constraints to jump model from dict
- Non-iterable argument to a function called by Julia `map`?
- Compute row sums and column sums efficiently in Julia
- Julia copying folder into an existing folder
- How to force Julia to use multiple threads for matrix multiplication?
- Can Revise.jl handle `ERROR: LoadError: invalid redefinition of constant`?
- Define piecewise function with automatic broadcasting in Julia
- julia Handling time difference in dataframe
- Fast tensor-dot on sparse arrays with GPU in any programming language?
- Why do allocations occur during broadcasting assignment to a preallocated array?
- Comparing RK4 to Euler method
- How to put even numbers from matrix in separate vector in Julia?
- Number of iterations performed by a for-loop in Julia
- Can I redefine a function, but still use the old definition within the new definition?
- how do I use analytical form as gradient with ! function?
- Julia equivalent to R `as.numeric()`
- Why does Float16(1.1)-Float16(1)=Float16(0.0996)?
- Julia manual and defining an infix operator
- In Julia, how to convert a unsigned number to a signed number like in C?
- Julia Package DataFrames 1.6.1 doesn't recognize old version DataFrames 1.3.4 file
- Convert Julia DataFrame to an array of bytes for compression
- How to suppress svg output by IJulia in Jupyter notebooks
- how to get the documentation of a julia function in a vscode notebook?
- Passing command-line arguments to a Pluto Notebook
- Julia gives "MethodError: rand!" when trying to optimize over a certain manifold using Manopt
- Why doesn't the loss calculated by Flux `withgradient` match what I have calculated?
- How to import Julia packages into Python
- Obtain the number of CPU cores in Julia