Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# tkralphs/gaussian_elim.py

Forked from num3ric/gaussian_elim.py
Last active Apr 6, 2021
 import numpy as np def GEPP(A, b, doPricing = True): ''' Gaussian elimination with partial pivoting. input: A is an n x n numpy matrix b is an n x 1 numpy array output: x is the solution of Ax=b with the entries permuted in accordance with the pivoting done by the algorithm post-condition: A and b have been modified. ''' n = len(A) if b.size != n: raise ValueError("Invalid argument: incompatible sizes between"+ "A & b.", b.size, n) # k represents the current pivot row. Since GE traverses the matrix in the # upper right triangle, we also use k for indicating the k-th diagonal # column index. # Elimination for k in range(n-1): if doPricing: # Pivot maxindex = abs(A[k:,k]).argmax() + k if A[maxindex, k] == 0: raise ValueError("Matrix is singular.") # Swap if maxindex != k: A[[k,maxindex]] = A[[maxindex, k]] b[[k,maxindex]] = b[[maxindex, k]] else: if A[k, k] == 0: raise ValueError("Pivot element is zero. Try setting doPricing to True.") #Eliminate for row in range(k+1, n): multiplier = A[row,k]/A[k,k] A[row, k:] = A[row, k:] - multiplier*A[k, k:] b[row] = b[row] - multiplier*b[k] # Back Substitution x = np.zeros(n) for k in range(n-1, -1, -1): x[k] = (b[k] - np.dot(A[k,k+1:],x[k+1:]))/A[k,k] return x if __name__ == "__main__": A = np.array([[1.,-1.,1.,-1.],[1.,0.,0.,0.],[1.,1.,1.,1.],[1.,2.,4.,8.]]) b = np.array([[14.],[4.],[2.],[2.]]) print GEPP(np.copy(A), np.copy(b), doPricing = False) print GEPP(A,b)

### alipars10 commented Feb 19, 2016

 Thanks for this, just a quick question. You mentioned that you fixed a bug. What exactly was the bug?

### jgcastro89 commented Oct 17, 2017

 Awesome job! I have forked your code and modularized it. I needed to import it and apply it to an existing project. https://gist.github.com/jgcastro89/49090cc69a499a129413597433b9baab
to join this conversation on GitHub. Already have an account? Sign in to comment