Skip to content

Instantly share code, notes, and snippets.

@glederrey
Created May 28, 2018 13:17
Show Gist options
  • Save glederrey/7fe6e03bbc85c81ed60f3585eea2e073 to your computer and use it in GitHub Desktop.
Save glederrey/7fe6e03bbc85c81ed60f3585eea2e073 to your computer and use it in GitHub Desktop.
Conjugate Gradient in Python
def conjgrad(A, b, x):
"""
A function to solve [A]{x} = {b} linear equation system with the
conjugate gradient method.
More at: http://en.wikipedia.org/wiki/Conjugate_gradient_method
========== Parameters ==========
A : matrix
A real symmetric positive definite matrix.
b : vector
The right hand side (RHS) vector of the system.
x : vector
The starting guess for the solution.
"""
r = b - np.dot(A, x)
p = r
rsold = np.dot(np.transpose(r), r)
for i in range(len(b)):
Ap = np.dot(A, p)
alpha = rsold / np.dot(np.transpose(p), Ap)
x = x + np.dot(alpha, p)
r = r - np.dot(alpha, Ap)
rsnew = np.dot(np.transpose(r), r)
if np.sqrt(rsnew) < 1e-8:
break
p = r + (rsnew/rsold)*p
rsold = rsnew
return x
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment