conjugate_gradient.py
c35860ee0260master

conjugate_gradient.py
View Options

			#!/usr/bin/env python3
			# -- coding: utf-8 --
			"""
			Created on Mon Oct 15 14:29:42 2018

			@author: sp4e
			"""

			import numpy as np
			from scipy.sparse.linalg import cg
			import time
			import math
			import numpy

			def conjugate_grad(A, b, x):
			n = len(b)
			r = b-np.dot(A, x)
			p = r
			rsold = np.dot(numpy.transpose(r), r)
			x=np.array([x])
			for i in range (0,10):
			Ap = np.dot(A, p)
			alpha = rsold / np.dot(numpy.transpose(p), Ap)
			#x = x+alpha * p
			x=np.append(x,[x[i]+alpha*p],axis=0)
			r = r-alpha * Ap
			rsnew = np.dot(numpy.transpose(r), r)
			beta = rsnew / rsold
			rsold= rsnew
			if abs(rsnew) < 1e-10:
			break
			p = beta * p + r
			x=np.einsum('ijk->ij',x)
			return x