Conjugate_gradient.cpp
1dc59d732f73master

Conjugate_gradient.cpp
View Options

			#include "Conjugate_gradient.hpp"
			#include <iostream>

			//Copy constructor//
			Conjugate_Gradient::Conjugate_Gradient(const Conjugate_Gradient &otherCholesky) {
			size=otherCholesky.size;
			A=otherCholesky.A;
			b=otherCholesky.b;
			x_0=otherCholesky.x_0;
			tol=otherCholesky.tol;
			}

			//Constructors
			Conjugate_Gradient::Conjugate_Gradient(double *m_A, double m_b, int m_size) {
			size=m_size;
			A=m_A;
			b=m_b;
			x_0=Initialize(size);
			tol=1e-12;
			}

			Conjugate_Gradient::Conjugate_Gradient(double *m_A, double m_b, double *m_x_0, int m_size) {
			size=m_size;
			A=m_A;
			b=m_b;
			x_0=m_x_0;
			tol=1e-12;
			}

			Conjugate_Gradient::Conjugate_Gradient(double *m_A, double m_b, double *m_x_0, double m_tol, int m_size) {
			size=m_size;
			A=m_A;
			b=m_b;
			x_0=m_x_0;
			tol=m_tol;
			}

			Conjugate_Gradient::Conjugate_Gradient(double *m_A, double m_b, double m_tol, int m_size) {
			size=m_size;
			A=m_A;
			b=m_b;
			x_0=Initialize(size);
			tol=m_tol;
			}


			double* Conjugate_Gradient::Solve() {
			int counter=0;
			double* x_1=x_0;
			double *r = Substract(b, Multiply(A, x_0, size, size, size), size);
			double *p_old = r;
			while(Multiply(r,r,size)>tol) {
			double r_old = Multiply(r, r, size);
			double alpha = r_old / Multiply(p_old, Multiply(A, p_old, size, size, size), size);
			x_1 = Add(x_1, Multiply(p_old, alpha, size), size);
			double *r_1 = Substract(r, Multiply(Multiply(A, p_old, size, size, size), alpha, size), size);
			double beta = Multiply(r_1, r_1, size) / r_old;
			double *p_new = Add(r_1, Multiply(p_old, beta, size), size);
			r=r_1;p_old=p_new;counter++;
			}
			n_iter=counter;
			return x_1;
			}