Homec4science
Diffusion PHPC_2019_Hanut (master) Browse

R8988/Parallel/Type1/MainParallelWOIO10239.cpp
10ba0eb21349master

Last Change
Raw File

MainParallelWOIO10239.cpp
View Options

//Main file for parallel code of Conjugate Gradient Method.
//Hanut Vemulapalli April 2019
#include <iostream>
#include <cmath>
#include <cassert>
#include "ArrayTools.h"
#include <fstream>
#include <mpi.h>
int main()
{
MPI_Init(NULL, NULL);
double time=MPI_Wtime();
int n =10239; //Ensure that this matches the data file sizes.
//
int world_size;
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
// Get the rank of the process
int my_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
//std::cout <<"My rank is"<<my_rank<<std::endl;
MPI_File Af;
MPI_File bf;
MPI_File x0f;
if(MPI_File_open(MPI_COMM_WORLD, "/home/vemulapa/PHPC_2019/code/Parallel/Matrix_A/A10239parallel.dat", MPI_MODE_RDONLY, MPI_INFO_NULL, &Af) != MPI_SUCCESS)
{
std::cerr<<"The file opening was not successful"<<std::endl;
};
if(MPI_File_open(MPI_COMM_WORLD, "/home/vemulapa/PHPC_2019/code/Parallel/Matrix_b/b10239parallel.dat", MPI_MODE_RDONLY, MPI_INFO_NULL, &bf)!= MPI_SUCCESS)
{
std::cerr<<"The file opening was not successful"<<std::endl;
};
if(MPI_File_open(MPI_COMM_WORLD, "/home/vemulapa/PHPC_2019/code/Parallel/Matrix_x0/x010239parallel.dat", MPI_MODE_RDONLY, MPI_INFO_NULL, &x0f)!= MPI_SUCCESS)
{
std::cerr<<"The file opening was not successful"<<std::endl;
};
MPI_Datatype etype;
MPI_Datatype file_type;
MPI_Status statusb;
etype=MPI_DOUBLE;
//Defining domain decomposition
int size_double=sizeof(etype);
int num_rows=sqrt(world_size);
int num_blocks=world_size/num_rows;
assert(world_size%num_blocks==0);
int block_size_r=n/num_rows;
int block_size=n/num_blocks;
int extra_column_size=n%num_blocks;
int extra_row_size=n%num_rows;
//Defining row and relative rank of processor in each row.
int pseudo_rank,my_row;
my_row=(int) my_rank/num_blocks; //row index starting from 0
pseudo_rank=(int) my_rank-(my_row*num_blocks);
int subarray_column_size,subarray_row_size;
if(pseudo_rank!=num_blocks-1)
{
subarray_column_size=block_size;
}
else
{
subarray_column_size=block_size+extra_column_size;
};
if(my_row != num_rows-1)
{
subarray_row_size=block_size_r;
}
else
{
subarray_row_size=block_size_r+extra_row_size;
};
//Reading A
int buff_size=subarray_row_size*subarray_column_size;
double* A=new double[buff_size];
const int startloc[2]={pseudo_rank*block_size,0};
const int n_size[2]={n,1};
const int subarray_size_c[2]={subarray_column_size,1};
if(MPI_Type_create_subarray(2,n_size, subarray_size_c, startloc, MPI_ORDER_C, MPI_DOUBLE, &file_type)!=MPI_SUCCESS)
{
std::cout<<"Subarray creation failed"<<std::endl;
}
MPI_Type_commit(&file_type);
MPI_File_set_view(Af, my_row*block_size_r*n*size_double, etype, file_type ,"/home/vemulapa/PHPC_2019/code/Parallel/Matrix_A/A10239parallel.dat", MPI_INFO_NULL );
MPI_File_read(Af,A, buff_size, MPI_DOUBLE, &statusb);
//Reading b and x0
double b[subarray_column_size],x0[subarray_column_size];
MPI_File_read_at(bf, pseudo_rank*block_size*size_double, &b, subarray_column_size, MPI_DOUBLE, &statusb);
MPI_File_read_at(x0f, pseudo_rank*block_size*size_double, &x0, subarray_column_size, MPI_DOUBLE, &statusb);
MPI_Type_free(&file_type);
MPI_File_close(&Af);
MPI_File_close(&bf);
MPI_File_close(&x0f);
//Creating new communicators with processors grouped row wise.
int row_rank,col_rank;
MPI_Comm row_comm, col_comm;
MPI_Comm_split( MPI_COMM_WORLD, my_row, my_rank, &row_comm);
MPI_Comm_rank(row_comm, &row_rank); //Same as pseudo rank :P
//Creating new communicator for first column processors.
MPI_Comm_split(MPI_COMM_WORLD, row_rank, my_rank, &col_comm);
MPI_Comm_rank(col_comm, &col_rank);
double ak =0 , bk = 0, my_nrk =0 , my_nrk1 =0 , my_dak = 0;
double rk[subarray_column_size], pk[subarray_column_size], xk[subarray_column_size],rk1[subarray_column_size],xk1[subarray_column_size],pk1[subarray_column_size];
double temp[subarray_row_size], redist_temp[subarray_column_size], full_temp[n];
int recvcounts[num_rows], displs[num_rows], sendcounts[num_blocks], displscol[num_blocks];
for(int i = 0; i < num_rows; ++i)
{
recvcounts[i]=block_size_r;
displs[i]=i*block_size_r;
}
recvcounts[num_rows-1]+=extra_row_size;
for(int j = 0; j < num_blocks; ++j)
{
sendcounts[j]=block_size;
displscol[j]=j*block_size;
}
sendcounts[num_blocks-1]+=extra_column_size;
MPI_Request request;
//rk=b-A*x0
MatVecProd(temp,A,x0,subarray_row_size,subarray_column_size);
if(!row_rank) {
MPI_Reduce(MPI_IN_PLACE, &temp, subarray_row_size, MPI_DOUBLE, MPI_SUM, 0, row_comm);
} else
{
MPI_Reduce(&temp, &temp, subarray_row_size, MPI_DOUBLE, MPI_SUM, 0, row_comm);
}
if(num_rows!=num_blocks)
{
if(!row_rank)
{
MPI_Gatherv(&temp, subarray_row_size, MPI_DOUBLE, &full_temp, recvcounts, displs, MPI_DOUBLE, 0, col_comm);
}
if(!col_rank)
{
MPI_Scatterv( &full_temp, sendcounts, displscol, MPI_DOUBLE, &redist_temp, subarray_column_size, MPI_DOUBLE, 0, row_comm);
}
}
else
{
if(!row_rank)
{
int send_rank=my_rank-(num_blocks-1)*my_row;
MPI_Isend(&temp,subarray_row_size,MPI_DOUBLE,send_rank,1,MPI_COMM_WORLD,&request);
}
if(!col_rank)
{
int recv_rank=my_rank*num_blocks;
MPI_Recv(&redist_temp,subarray_column_size,MPI_DOUBLE,recv_rank,1,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
}
}
//Send from fist row to other rows
MPI_Bcast( &redist_temp, subarray_column_size, MPI_DOUBLE, 0, col_comm);
subtract(rk,b,redist_temp,subarray_column_size);
//pk=rk
copy(pk,rk,subarray_column_size);
//xk=x0
copy(xk,x0,subarray_column_size);
double res=10.0;
int k=1;
//
while (res > 1e-15)
{
// if(!my_rank)
// {
// std::cout<<"Iteration: "<<k<<" Residual "<< res<< std::endl;
//}
//ak=rk.'rk/(pk.'*A*pk)
MatVecProd(temp, A, pk, subarray_row_size, subarray_column_size);
if(!row_rank) {
MPI_Reduce(MPI_IN_PLACE, &temp, subarray_row_size, MPI_DOUBLE, MPI_SUM, 0, row_comm);
} else
{
MPI_Reduce(&temp, &temp, subarray_row_size, MPI_DOUBLE, MPI_SUM, 0, row_comm);
}
if(num_rows!=num_blocks)
{
if(!row_rank)
{
MPI_Gatherv(&temp, subarray_row_size, MPI_DOUBLE, &full_temp, recvcounts, displs, MPI_DOUBLE, 0, col_comm);
}
if(!col_rank)
{
MPI_Scatterv( &full_temp, sendcounts, displscol, MPI_DOUBLE, &redist_temp, subarray_column_size, MPI_DOUBLE, 0, row_comm);
}
}
else
{
if(!row_rank)
{
int send_rank=my_rank-(num_blocks-1)*my_row;
MPI_Isend(&temp,subarray_row_size,MPI_DOUBLE,send_rank,1,MPI_COMM_WORLD,&request);
}
if(!col_rank)
{
int recv_rank=my_rank*num_blocks;
MPI_Recv(&redist_temp,subarray_column_size,MPI_DOUBLE,recv_rank,1,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
}
}
//Send from fist row to other rows
MPI_Bcast(&redist_temp, subarray_column_size, MPI_DOUBLE, 0, col_comm);
//Calculating ak
if (!my_row) {
my_nrk = dot(rk, rk, subarray_column_size);
my_dak = dot(pk, redist_temp, subarray_column_size);
if(!row_rank) {
MPI_Reduce(MPI_IN_PLACE, &my_nrk, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
MPI_Reduce(MPI_IN_PLACE, &my_dak, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
} else
{
MPI_Reduce(&my_nrk, &my_nrk, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
MPI_Reduce(&my_dak, &my_dak, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
}
}
if (!my_rank) {
ak = my_nrk / my_dak;
}
MPI_Bcast(&ak, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
//xk1=xk+ak*pk
scalar(ak, pk, subarray_column_size);
add(xk1, xk, pk, subarray_column_size);
//rk1=rk-ak*A*pk;
scalar(ak, redist_temp, subarray_column_size);
subtract(rk1, rk, redist_temp, subarray_column_size);
//Calculating bk
if (!my_row) {
my_nrk1 = dot(rk1, rk1, subarray_column_size);
// Insert code to send to processor 0 and add partial sums
if(!row_rank) {
MPI_Reduce(MPI_IN_PLACE, &my_nrk1, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
} else
{
MPI_Reduce(&my_nrk1, &my_nrk1, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
}
}
if (!my_rank) //Only processor 0 does this computation. Only 0 has full my_nrk1 and my_nrk
{
bk = my_nrk1 / my_nrk;
}
MPI_Bcast(&bk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
//pk1=rk1+bk*pk;
scalar(bk / ak, pk, subarray_column_size); //Since pk=ak*pk now. Losing some accuracy here because of division.
add(pk1, rk1, pk, subarray_column_size);
//rk=rk1;
copy(rk, rk1, subarray_column_size);
//pk=pk1;
copy(pk, pk1, subarray_column_size);
//xk=xk1;
copy(xk, xk1, subarray_column_size);
//Calculating residual
if (!my_row) {
res = norm(rk1, subarray_column_size);
//Code to reduce res values to processor
if(!row_rank) {
MPI_Reduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
} else
{
MPI_Reduce(&res, &res, 1, MPI_DOUBLE, MPI_SUM, 0, row_comm);
}
}
if (!my_rank) {
res = sqrt(res);
}
//Code to broadcast res values. BCast is collective
MPI_Bcast(&res, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
k += 1;
}
MPI_Comm_free(&col_comm);
MPI_Comm_free(&row_comm);
time=time-MPI_Wtime();
MPI_Finalize();
std::cout<< time <<std::endl;
return 0;
}
c4science · Help