Solving a system of linear equations using Gauss Siedel method, where the only difference from Jacobi's method is that the new values of x are applied to the subsequent equations in the same iteration

gauss-siedel.py

# -*- coding: utf-8 -*-
"""
Created on Sun Jan 19 12:47:08 2020

@author: DASA0
"""
""" Gauss- Sidel's method"""

import numpy as np

a = np.array([[4, 1, 2, -1],
           [3, 6, -1, 2],
           [2, -1, 5, -3],
           [4, 1, -3, -8]])
b = np.array([2, -1, 3, 2])
(n,) = np.shape(b)
x = np.full(n, 1.0, float) # initial guess
xdiff = np.empty(n, float)
iterlimit = 100
tolerance = 1.0e-4

# iterations:
for iter in range(iterlimit):
    for i in range(n):
        s = 0
        for j in range(n):
            if j != i:
                s += a[i, j]*x[j]
        xnew = -1/a[i,i] * (s - b[i])
        xdiff = abs(xnew - x[i])
        x[i] = xnew
    if (xdiff < tolerance).all():
        break

print('Number of iterations: %d '% (iter + 1))
print('The solution for the system of linear equations is:')
print(x)