-
-
Save tkralphs/7554375 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
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
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Thanks for this, just a quick question. You mentioned that you fixed a bug. What exactly was the bug?