#Gaussian Elimination with pivoting
import numpy as np

def gauss(a,b,x):
  N=len(a)
  for i in range(N):
#pivot
    maxr=i
    maxa=abs(a[i,i])
    for j in range(i+1,N):
      if abs(a[j,i]) > maxa:
        maxa=abs(a[j,i])
        maxr=j
    if maxa == 0.0:  print( "matrix is singular")
    if maxr != i:
      rtemp=a[maxr].copy()
      a[maxr]=a[i]
      a[i]=rtemp
      temp=b[maxr]
      b[maxr]=b[i]
      b[i]=temp

#reduce row
    for j in range(i+1,N):
      for k in range(i+1,N):
        a[j,k]-=a[j,i]*a[i,k]/a[i,i]
      b[j]-=b[i]*a[j,i]/a[i,i]

# back substitution
  temp=range(N)
  temp.reverse()
  for i in temp:
    for j in range(i+1,N): 
      b[i]-=a[i,j]*x[j]
    x[i]=b[i]/a[i,i]




a=np.loadtxt("matrix")
b=np.loadtxt("vector")
x=np.zeros(len(a))
gauss(a,b,x)
print(x)