# Gauss-Seidel method with python from wikipedia
# https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method
import numpy as np
iteration_limit = 1000
# initialize the matrix
A = np.array([[10., -1., 2., 0.],
[-1., 11., -1., 3.],
[2., -1., 10., -1.],
[0., 3., -1., 8.]])
# initialize the RHS vector
b = np.array([6.0, 25.0, -11.0, 15.0])
print("System of equations:")
for i in range(A.shape[0]):
row = ["{0:3g}*x{1}".format(A[i, j], j+1) for j in range(A.shape[1])]
print("[{0}] = [{1:3g}]".format("+".join(row), b[i]))
x = np.zeros_like(b)
for it_count in range(1, iteration_limit):
x_new = np.zeros_like(x)
print("Iteration {0}:{1}".format(it_count, x))
for i in range(A.shape[0]):
s1 = np.dot(A[i, :i], x_new[:i])
s2 = np.dot(A[i, i+1:], x[i+1:])
x_new[i] = (b[i] - s1 - s2) / A[i, i]
if np.allclose(x, x_new, rtol=1e-8):
break
x = x_new
print("Solution: {0}".format(x))
error = np.dot(A, x) - b
print("Error: {0}".format(error))
TODO: wavefront parallel method with openmp