import numpy as np

A = np.array([[7, 1, 2], [1, 8, 3], [2, 3, 9]])
b = np.array([10, 8, 6])

LU = np.zeros(A.shape)
D = np.zeros(A.shape)
for y in range(A.shape[1]):
    for x in range(A.shape[0]):
        if x == y:
            D[x][y] = A[x][y]
        else:
            LU[x][y] = A[x][y]

x = np.array([1000, 2000, 3000])

for n in range(100):
    # 1. calculate error
    e = np.linalg.norm(b - A @ x)
    print(n, x, e)

    # 2. update
    x = np.linalg.inv(D) @ (b - LU @ x)

e = np.linalg.norm(b - A @ x)
print("answer is %s (error = %s)" % (x, e))