Page MenuHomec4science
Files F120910258

CG_Serial.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, Jul 7, 19:22

CG_Serial.cpp
View Options

//
// Created by shernand on 29/05/19.
//
/**
* @file CG_Serial.cpp
* @author Sergio Hernandez
* This file is part of the Conjugate Gradient Project
*
* This class implements the serial version of the Conjugate Gradient
*
*/
#include <ctime>
#include <fstream>
#include "CG_Serial.hpp"
/**
* Complete constructor of the CG solver.
* @param matrixA std::vector<double> NxN symmetric definite positive
* @param vectorB Nx1 std::vector<double> right side vector
* @param tol double Tolerance for our stopping criterion
* @param maxIter int Maximum number of iterations of the stopping criterion
*/
CG_Serial::CG_Serial(std::vector<double> const &matrixA, std::vector<double> const &vectorB, double const &tol,
int const &maxIter) {
setMatrixA(matrixA);
setVectorB(vectorB);
setTol(tol);
setMaxIterations(maxIter);
}
/**
* Computes the CG solver
* @param x_0 Nx1 std::vecotr<double> the initial guess for the solution
* @return std::tuple<std::vector<double> x, int n_iterations> The tuple <solution, (Nx1 double) and how many iterations
* it took.
*/
void CG_Serial::computeCG(std::vector<double> &x_0) {
std::vector<double> b = getVectorB();
std::vector<double> mA = getMatrixA();
int max_iter = getMaxIterations();
double tol = getTol();
std::vector<double> x_k = x_0;
std::vector<double> r_k = matrixVector(mA,x_k);
r_k = vectorScalarMul(r_k,-1.0);
r_k = vectorVectorSum(b,r_k);
std::vector<double> d_k = r_k;
double norm_res = vectorVectorDot(r_k,r_k);
int k = 0;
while ((norm_res>tol)&&(k<max_iter)){
std::vector<double> A_dk = matrixVector(mA,d_k);
double dk_A_dk = vectorVectorDot(d_k,A_dk);
double alpha_k = norm_res/dk_A_dk;
x_k = vectorVectorSum(x_k,vectorScalarMul(d_k,alpha_k));
A_dk = vectorScalarMul(A_dk,-alpha_k);
r_k= vectorVectorSum(r_k,A_dk);
double norm_r_k_1 = vectorVectorDot(r_k,r_k);
double beta_k_1 = norm_r_k_1/norm_res;
d_k = vectorVectorSum(r_k,vectorScalarMul(d_k,beta_k_1));
norm_res=norm_r_k_1;
k+=1;
}
setSolution(x_k);
setNIterations(k);
}
/**
* Saves the solution vector, and the number of iterations into a file.
* @param vector
* @param n_iterations
* @param filename
*/
void CG_Serial::saveVector(std::string &filename){
std::ofstream myfile;
std::vector<double> vector = getSolution();
int n_iterations = getNIterations();
int n = vector.size();
//Opens the file
myfile.open (filename);
myfile << n_iterations<<"\n";
for(int i=0;i<n;i++){
myfile << vector[i]<<"\n";
}
myfile.close(); //Closes the file
}
/**
* The operation scalar*vector
* @param pVector1 Nx1 std::vector<double>
* @param pScalar double
* @return mul_vector Nx1 std::vector<double> mult element wise of the scalar and the vector
*/
std::vector<double> CG_Serial::vectorScalarMul(const std::vector<double> &pVector1, const double pScalar) {
std::vector<double> mul_vector = pVector1;
for(int i=0;i<pVector1.size();i++){
mul_vector[i]= pVector1[i]*pScalar;
}
return mul_vector;
}
/**
* The operation vector+vector
* @param pVector1 Nx1 std::vector<double>
* @param pVector2 Nx1 std::vector<double>
* @return sum_vector Nx1 std::vector<double> sum element wise of the two vectors
*/
std::vector<double> CG_Serial::vectorVectorSum(const std::vector<double> &pVector1, const std::vector<double> pVector2) {
std::vector<double> sum_vector = pVector1;
for(int i=0;i<pVector1.size();i++){
sum_vector[i]= pVector1[i]+pVector2[i];
}
return sum_vector;
}
/**
* The operation vector*vector
* @param pVector1 Nx1 std::vector<double>
* @param pVector2 Nx1 std::vector<double>
* @return double the resulting dot product
*/
double CG_Serial::vectorVectorDot(const std::vector<double> &pVector1, const std::vector<double> pVector2) {
double dot_product = 0;
for(int i=0;i<pVector1.size();i++){
dot_product+= pVector1[i]*pVector2[i];
}
return dot_product;
}
/**
* The operation Matrix*Vector. This is the most expensive operation.
* @param pMatrix NxN matrix std::vector<double>
* @param pVector Nx1 vector std::vector<double>
* @return prod_vecotr Nx1 std::vector<double>
*/
std::vector<double> CG_Serial::matrixVector(const std::vector<double> &pMatrix, const std::vector<double> pVector) {
std::vector<double> prod_vector = pVector;
for(int i=0;i<pVector.size();i++){
double sum_term = 0.0;
for(int j=0;j<pVector.size();j++){
sum_term+= pMatrix[i+pVector.size()*j]*pVector[j];
}
prod_vector[i]=sum_term;
}
return prod_vector;
}
//-------------------Getters-Setters----------------------
const std::vector<double> &CG_Serial::getMatrixA() const {
return matrixA;
}
void CG_Serial::setMatrixA(const std::vector<double> &matrixA) {
CG_Serial::matrixA = matrixA;
}
const std::vector<double> &CG_Serial::getVectorB() const {
return vectorB;
}
void CG_Serial::setVectorB(const std::vector<double> &vectorB) {
CG_Serial::vectorB = vectorB;
}
double CG_Serial::getTol() const {
return tol;
}
void CG_Serial::setTol(double tol) {
CG_Serial::tol = tol;
}
int CG_Serial::getMaxIterations() const {
return max_iterations;
}
void CG_Serial::setMaxIterations(int maxIterations) {
max_iterations = maxIterations;
}
int CG_Serial::getNIterations() const {
return n_iterations;
}
void CG_Serial::setNIterations(int nIterations) {
n_iterations = nIterations;
}
const std::vector<double> &CG_Serial::getSolution() const {
return vectorSolution;
}
void CG_Serial::setSolution(const std::vector<double> &vectorSolution) {
CG_Serial::vectorSolution = vectorSolution;
}

Event Timeline

c4science · Help