Gauss

# 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