text: #---------------------------------------------#
  #  Gauss-Seidel Linear System Numeric Solver  #
  #                  09/05/21                   #
  #---------------------------------------------#
  
  import numpy as np
  
  #Setup 
  n = 3
  aMat = np.array([[4,-2,-1],[-5,-15,-9],[2,2,-5]])
  xMat = np.zeros(n)
  xMat_prev = np.zeros(n)
  bMat = np.array([1,-9,-8])
  betaMat = np.zeros(n)
  errorVal = 1e-3
  errorCalc = 0.0
  itNum = 0
  itMax = 1000
  
  #Matrix check
  for i in range(0, n):
      if aMat[i, i] == 0:
          print("Matrix has at least one element in main diagonal which equals 0")
  
  #Sassenfeld Criteria
  for i in range(0, n):
      for j in range(0 , i):
          betaMat[i] = betaMat[i] + np.abs(aMat[i, j])*betaMat[j]
      for j in range(1+i, n):
          betaMat[i] = betaMat[i] + np.abs(aMat[i, j])
      betaMat[i] = betaMat[i] / np.abs(aMat[i, i])
  if np.max(betaMat) < 1.:
      print("Sassenfeld Criteria fulfilled. Method will converge. Max beta = " + str(np.max(betaMat)))
      multiplier = (np.max(betaMat)/ 1 + np.max(betaMat))
  else:
      print("Sassenfeld Criteria not fulfilled. Cannot guarantee convergence.")
      multiplier = 1000
      
  #Iterative process
  for k in range(0, itMax):
      print(str(xMat[:]))
      for i in range(0, n):
          xMat[i] = bMat[i]
          for j in range(0, i):
              xMat[i] = xMat[i] - aMat[i, j]*xMat[j]
          for j in range(i+1, n):
              xMat[i] = xMat[i] - aMat[i, j]*xMat_prev[j]
          xMat[i] = xMat[i] / aMat[i, i] 
      #Error calculation
      errorCalc = multiplier * np.max(np.abs(xMat - xMat_prev))
      itNum += 1
      print("Iteration nº" + str(itNum) + " // Current errorValue = " + str(errorCalc))
      if(errorCalc < errorVal):
          break
      xMat_prev = np.copy(xMat)
  
  #Final verification
  finalVal = np.matmul(aMat, xMat)
  print("Final results are: " + str(xMat[:]))
  print("With residue: " + str(np.max(np.abs(bMat - finalVal))))
  