Page MenuHomec4science
Files F102733957

cg.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Feb 23, 15:30

cg.cpp
View Options

#include "cg.h"
using namespace std;
double cg::calcAlpha(double numerator,vector<double> Ap,vector<double> p){
vector<double> p_tr=p;
double denominator = operation::dotProduct(p_tr,Ap);
return numerator/denominator;
}
double cg::calcBeta(double denominator,vector<double> r,vector<double> r_new){
vector<double> r_tr=r;
vector<double> r_tr_new=r_new;
double numerator=operation::dotProduct(r_tr_new,r_new);
return numerator/denominator;
}
vector<double> cg::conjugateGradient(vector<double> A,vector<double> b,double accuracy){
int vectorSize=b.size();
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
int numProc;
MPI_Comm_size(MPI_COMM_WORLD, &numProc);
int rows=b.size();
vector<double> x(vectorSize);
vector<double> r;
r=b;
vector<double> p=r;
//keeps track of iterations
int i=0;
double pBuffer[vectorSize];
for(int i=0;i<p.size();i++){
pBuffer[i]=p[i];
}
int done=0;
int sc[numProc];
operation::sizeCounts(rows,numProc,sc);
int disp[numProc];
operation::displacement(numProc,sc,disp);
double A_arr[rows*rows];
for(int i=0;i<rows*rows;i++){
A_arr[i]=A[i];
}
while(!done){
for(int i=0;i<p.size();i++){
p[i]=pBuffer[i];
}
vector<double> Ap=operation::matrMultiply(A_arr,p,sc,disp);
if(rank==0){
vector<double> r_tr=r;
//r transposed multiplied by r
double r_tr_r=operation::dotProduct(r_tr,r);
double alpha=calcAlpha(r_tr_r,Ap,p);
//p scaled by alpha
vector<double> alpha_p=p;
for(int i=0;i<alpha_p.size();i++){
alpha_p[i]=alpha_p[i]*alpha;
}
x=operation::vecAdd(x,alpha_p);
vector<double> alpha_Ap=Ap;
for(int i=0;i<alpha_Ap.size();i++){
alpha_Ap[i]=alpha_Ap[i]*alpha;
}
vector<double> r_new=operation::vecSubtract(r,alpha_Ap);
double beta = calcBeta(r_tr_r,r,r_new);
vector<double> Bp=p;
for(int i=0;i<Bp.size();i++){
Bp[i]=Bp[i]*beta;
}
p=operation::vecAdd(r_new,Bp);
r=r_new;
for(int i=0;i<p.size();i++){
pBuffer[i]=p[i];
}
done = operation::vector_norm(r)<accuracy;
}
MPI_Bcast(&done,1, MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(pBuffer,vectorSize, MPI_DOUBLE,0,MPI_COMM_WORLD);
i=i+1;
}
if(rank==0){
cout<<"\nIterations: "<<i<<"\n";
cout<<"Final Error: "<<operation::vector_norm(r)<<"\n";
}
return x;
}

Event Timeline

c4science · Help