Homec4science
Diffusion conjugate gradient project (master) Browse

R6737/Admahl's law/Admahl_law.cpp
eb88ccadd1a9master

Last Change
Raw File

Admahl_law.cpp
View Options

/* This code is written by Mehrdad Kiani 2018*/
/* Calculating the alpha in Admahl's law*/
/* Comments are exist in the serial conjugate gradient code*/
/* PHPC_2018 */
#include<iostream>
#include<cmath>
#include <chrono>
using namespace std;
void initialization (double* b, double** a,double* x, int n)
{
for(int i=0;i<n;i++) {
b[i] = 1.0;
x[i] = 0.0;
for(int j=0;j<n;j++) {
a[i][j] = 0.0;
if(j==(i+1)){
a[i][j] = -1.0;
}
if(j==(i-1)){
a[i][j] = -1.0;
}
if(j==(i+2)){
a[i][j] = -1.0;
}
if(j==(i-2)){
a[i][j] = -1.0;
}
if(i==j){
a[i][j] = (n+1)*1.0;
}
}
}
}
void precalculation (double* h, double* x, double** a, double* r, double* b, double* d,int n)
{
for(int i=0;i<n;i++) {
h[i]=0.;
for(int j=0;j<n;j++) {
h[i] += a[i][j]*x[j];
}
r[i] = b[i]-h[i];
d[i] = r[i];
}
}
void printResults(double* x)
{
cout << "\nThe Answer is = " << " ";
for(int i=0;i<10;i++) {
cout << x[i] <<" ";
}
cout<<"\n";
}
double CalAlpha(double* dd, double* d, double** a, double* r,double alpha, int n)
{
double rr = 0.0, ddd = 0.0;
for(int j=0;j<n;j++) {
dd[j] = 0.0;
for(int i=0;i<n;i++) {
dd[j] += d[i]*a[i][j];
}
rr += r[j] * r[j];
ddd += dd[j] * d[j];
}
alpha = rr / ddd;
return alpha;
}
void secondfun(double* x, double* d, double** a, double* r, double* r_old, double* sum, double alpha, int n)
{
for(int i=0;i<n;i++) {
x[i] = x[i]+alpha*d[i];
r_old[i] = r[i];
sum[i] = 0.0;
for(int j=0;j<n;j++) {
sum[i] += alpha*a[i][j]*d[j];
}
r[i] = r_old[i] - sum[i];
}
}
int n;
int main()
{
cout<<"How size you disere to be your problem ? "<<endl;
cin>> n;
chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
double rrr = 0.0, roldr = 0.0, alpha = 0.0, beta = 1.0, beta_save=0 ;
int i,j,m = 0, m_max=5;
double accum = 0.0, norm = 10.0, tol = 1.0E-8;
double *b = new double[n];
double *x = new double[n];
double *dd = new double[n];
double *d = new double[n];
double *r = new double[n];
double *r_old = new double[n];
double *sum = new double[n];
double *h = new double[n];
double *residual = new double[n];
double **a = new double *[n];
for (int i=0; i< n; i++) {
a[i] = new double[n];
}
chrono::high_resolution_clock::time_point parallel_1_begin = chrono::high_resolution_clock::now();
initialization(b, a, x, n);
precalculation(h, x, a, r, b, d, n);
while(m <= m_max) {
alpha = CalAlpha(dd,d, a,r,alpha, n);
cout<<"\n ALPHA: "<<alpha<<endl;
secondfun(x,d,a, r,r_old,sum, alpha, n);
rrr = 0.0; roldr = 0.0;
for(i=0;i<n;i++) {
rrr += r[i]*r[i];
roldr += r_old[i]*r_old[i];
}
if(beta> 1.e-7){
beta = rrr / roldr;
beta_save = beta;
cout<<"beta : "<<beta<<endl;
}
if(beta<= 1.e-7){
beta = beta_save;
cout<<"beta : "<<beta<<endl;
}
for(i=0;i<n;i++) {
d[i] = r[i]+beta*d[i];
}
m++;
}
chrono::high_resolution_clock::time_point parallel_1_end = chrono::high_resolution_clock::now();
chrono::duration<double> time_span_parallel_1 = chrono::duration_cast<chrono::duration<double>>(parallel_1_end - parallel_1_begin);
for(int i =0; i<n; i++) {
delete[] a[i];
}
delete[] a; delete[] h; delete[] r;
delete[] d; delete[] b; delete[] r_old;
delete[] sum; delete[] dd; delete[] residual;
// printResults(x);
delete[] x;
double time_parallel=time_span_parallel_1.count();
chrono::high_resolution_clock::time_point t2 = chrono::high_resolution_clock::now();
chrono::duration<double> time_span = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "\n"<< " Total loop took: " << time_span.count()*1000000 << " usecs.\n";
cout << "\n"<< "Alpha according to Adhaml's law: " << (time_span.count()-time_parallel)/(time_span.count())<<"\n";
return 0;
}
c4science · Help