Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <mpi.h>
- #include <omp.h>
- #include <iomanip>
- #include <iostream>
- #include <vector>
- #include <map>
- #include <cmath>
- #include <cstdlib>
- #include <algorithm>
- #include <fstream>
- #include<sstream>
- using namespace std;
- class MPIComputations{
- int M, N;
- int cur_block_size_x, cur_block_size_y, cur_block_global_offset_x, cur_block_global_offset_y;
- int Gx, Gy;
- int cur_block_global_coors_x, cur_block_global_coors_y;
- double h_x, h_y;
- int global_block_size_x, global_block_size_y;
- int my_rank;
- double current_norm;
- MPI_Comm comm;
- double x_left, x_right, y_bottom, y_top;
- ofstream my_file;
- vector<double> internal_data;
- vector<double> old_internal_data;
- vector<double> external_data[4];
- int IsNodeInternalCornerOrSide(int current_node_global_offset_x, int current_node_global_offset_y){
- //corners
- //left bottom corner
- if (current_node_global_offset_x == 0 && current_node_global_offset_y == 0){
- return 2;
- }
- //left top corner
- if (current_node_global_offset_x == 0 && current_node_global_offset_y == N){
- return 4;
- }
- //right bottom corner
- if (current_node_global_offset_x == M && current_node_global_offset_y == 0){
- return 6;
- }
- //right top corner
- if (current_node_global_offset_x == M && current_node_global_offset_y == N){
- return 8;
- }
- //sides
- //left side
- if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
- current_node_global_offset_x == 0){
- return 1;
- }
- //right side
- if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
- current_node_global_offset_x == M){
- return 3;
- }
- //bottom side
- if (current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1 &&
- current_node_global_offset_y == 0){
- return 5;
- }
- //top side
- if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1 &&
- current_node_global_offset_y == N)){
- return 7;
- }
- //internal
- if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1) &&
- (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1)){
- return 0;
- }
- return -1;
- }
- double k(double x, double y) {
- return 1 + pow(x + y, 2);
- }
- double q(double x, double y) {
- return 1;
- }
- double u(double x, double y) {
- return 2.0 / (1 + pow(x, 2) + pow(y, 2));
- }
- // psi_R(x, y) = k(A2, y) * du/dx(A2, y)
- double psi_R(double y) {
- return (-12) * (pow((y + 3), 2) + 1) / pow((pow(y, 2) + 10), 2);
- }
- // psi_L(x, y) = -k(A1, y) * du/dx(A1, y)
- double psi_L(double y) {
- return (-8) * (pow((y - 2), 2) + 1) / pow((pow(y, 2) + 5), 2);
- }
- // psi_T(x, y) = k(x, B2) * du/dy(x, B2)
- double psi_T(double x) {
- return (-16) * (pow((x + 4), 2) + 1) / pow((pow(x, 2) + 17), 2);
- }
- // psi_B(x, y) = -k(x, B1) * du/dy(x, B1)
- double psi_B(double x) {
- return (-4) * (pow((x - 1), 2) + 1) / pow((pow(x, 2) + 2), 2);
- }
- // right-part function of Poisson equation
- double F(double x, double y) {
- return 2 * (pow(x,4) + pow(y,4) + 2 * (pow(x,2) + 3) * pow(y,2) + 6 * pow(x,2) + 16*x*y + 5)
- / pow((1 + pow(x, 2) + pow(y, 2)), 3);
- }
- //inner_product(A[i], internal_data)
- double ComputeMagicInnerProductA_iw (int current_node_global_offset_x, int current_node_global_offset_y){
- int glob_x = current_node_global_offset_x;
- int glob_y = current_node_global_offset_y;
- double result = 0.0;
- map <string,bool> neighbours = {
- {"left", true},
- {"right", true},
- {"bottom", true},
- {"top", true}
- };
- double left_neighbour = 0.0, right_neighbour = 0.0, bottom_neighbour = 0.0, top_neighbour = 0.0, this_node = 0.0;
- double left_coeff = 1.0, right_coeff = 1.0, bottom_coeff = 1.0, top_coeff = 1.0, this_coeff = 1.0;
- switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
- case 2:
- //left bottom corner
- neighbours["left"] = false;
- neighbours["bottom"] = false;
- break;
- case 4:
- //left top corner
- neighbours["left"] = false;
- neighbours["top"] = false;
- break;
- case 6:
- //right bottom corner
- neighbours["right"] = false;
- neighbours["bottom"] = false;
- break;
- case 8:
- //right top corner
- neighbours["right"] = false;
- neighbours["top"] = false;
- break;
- case 1:
- //left side
- neighbours["left"] = false;
- break;
- case 3:
- //right side
- neighbours["right"] = false;
- break;
- case 5:
- //bottom side
- neighbours["bottom"] = false;
- break;
- case 7:
- //top side
- neighbours["top"] = false;
- break;
- case 0:
- //internal
- break;
- default:
- cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y<<endl;
- }
- if (!neighbours["left"]){
- right_coeff = 2.0;
- left_coeff = 0.0;
- }
- if (!neighbours["right"]){
- left_coeff = 2.0;
- right_coeff = 0.0;
- }
- if (!neighbours["bottom"]){
- top_coeff = 2.0;
- bottom_coeff = 0.0;
- }
- if (!neighbours["top"]){
- bottom_coeff = 2.0;
- top_coeff = 0.0;
- }
- if (neighbours["left"]){
- left_coeff *= -k(x_left + (glob_x - 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
- left_neighbour = Get(glob_x - 1, glob_y);
- }
- if (neighbours["right"]){
- right_coeff *= -k(x_left + (glob_x + 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
- right_neighbour = Get(glob_x + 1, glob_y);
- }
- if (neighbours["bottom"]){
- bottom_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y - 0.5) * h_y) / pow(h_y, 2);
- bottom_neighbour = Get(glob_x, glob_y - 1);
- }
- if (neighbours["top"]){
- top_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y + 0.5) * h_y) / pow(h_y, 2);
- top_neighbour = Get(glob_x, glob_y + 1);
- }
- this_coeff = q(x_left + glob_x * h_x, y_bottom + glob_y * h_y) - left_coeff - right_coeff - bottom_coeff - top_coeff;
- this_node = Get(glob_x, glob_y);
- result = left_coeff * left_neighbour +
- right_coeff * right_neighbour +
- bottom_coeff * bottom_neighbour +
- top_coeff * top_neighbour +
- this_coeff * this_node;
- return result;
- }
- double GetNodeFromB(int current_node_global_offset_x, int current_node_global_offset_y) {
- int glob_x = current_node_global_offset_x;
- int glob_y = current_node_global_offset_y;
- double result = 0.0;
- switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
- case 2:
- //left bottom corner
- result = F(x_left, y_bottom) + 2.0 / h_x * psi_L(y_bottom) + 2.0 / h_y * psi_B(x_left);
- break;
- case 4:
- //left top corner
- result = F(x_left, y_top) + 2.0 / h_x * psi_L(y_top) + 2.0 / h_y * psi_T(x_left);
- break;
- case 6:
- //right bottom corner
- result = F(x_right, y_bottom) + 2.0 / h_x * psi_R(y_bottom) + 2.0 / h_y * psi_B(x_right);
- break;
- case 8:
- //right top corner
- result = F(x_right, y_top) + 2.0 / h_x * psi_R(y_top) + 2.0 / h_y * psi_T(x_right);
- break;
- case 1:
- //left side
- result = F(x_left, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_L(y_bottom + glob_y * h_y);
- break;
- case 3:
- //right side
- result = F(x_right, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_R(y_bottom + glob_y * h_y);
- break;
- case 5:
- //bottom side
- result = F(x_left + glob_x * h_x, y_bottom) + 2.0 / h_y * psi_B(x_left + glob_x * h_x);
- break;
- case 7:
- //top side
- result = F(x_left + glob_x * h_x, y_top) + 2.0 / h_y * psi_T(x_left + glob_x * h_x);
- break;
- case 0:
- //internal
- result = F(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
- break;
- default:
- cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y <<endl;
- }
- return result;
- }
- double GetNodeFromExact(int current_node_global_offset_x, int current_node_global_offset_y) {
- int glob_x = current_node_global_offset_x;
- int glob_y = current_node_global_offset_y;
- return u(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
- }
- void ComputeMatrixR(){
- /*if(my_rank == 0)
- cout << "[INFO]: Computation of matrix r started"<<endl;
- */
- vector<double> r_tmp_matrix (cur_block_size_x*cur_block_size_y, 0.0);
- //#pragma omp parallel for
- for(int i = 0; i < cur_block_size_x; ++i){
- for(int j = 0; j < cur_block_size_y; ++j){
- int current_node_global_offset_x = GetGlobalX(i),
- current_node_global_offset_y = GetGlobalY(j);
- int glob_x = current_node_global_offset_x,
- glob_y = current_node_global_offset_y;
- r_tmp_matrix [ j + cur_block_size_y*i ] = ComputeMagicInnerProductA_iw(glob_x,glob_y) - GetNodeFromB(glob_x, glob_y);
- }
- }
- //#pragma omp parallel for
- for(int i = 0; i < cur_block_size_x; ++i){
- for(int j = 0; j < cur_block_size_y; ++j){
- old_internal_data[ j + cur_block_size_y*i ] = internal_data[ j + cur_block_size_y*i];
- internal_data[ j + cur_block_size_y*i ] = r_tmp_matrix[ j + cur_block_size_y*i];
- }
- }
- //old_internal_data = internal_data;
- //internal_data = r_tmp_matrix;
- SyncMPI();
- }
- double ComputeTauAndStopCase(bool &should_i_stop){
- double local_Ar_r_inner_product_sum = 0.0;
- double local_Ar_Ar_inner_product_sum = 0.0;
- double global_Ar_r_inner_product_sum = 0.0;
- double global_Ar_Ar_inner_product_sum = 0.0;
- double local_r_norm = 0.0;
- double global_r_norm = 0.0;
- //#pragma omp parallel for
- for(int i = 0; i < cur_block_size_x; ++i){
- for(int j = 0; j < cur_block_size_y; ++j){
- double rho = 1.0;
- int current_node_global_offset_x = GetGlobalX(i),
- current_node_global_offset_y = GetGlobalY(j);
- int glob_x = current_node_global_offset_x,
- glob_y = current_node_global_offset_y;
- double tmp_Ar_i_j = ComputeMagicInnerProductA_iw(glob_x, glob_y);
- switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
- case 2:
- case 4:
- case 6:
- case 8:
- //angle
- rho = 0.25;
- break;
- case 1:
- case 3:
- case 5:
- case 7:
- //side
- rho = 0.5;
- break;
- case 0:
- //internal
- rho = 1.0;
- break;
- default:
- cout << "[ERROR]: Bad global coords compute tau. Global:" << glob_x << " " << glob_y << endl;
- }
- double tmp_cur_node_value = Get(glob_x, glob_y);
- local_Ar_r_inner_product_sum += rho * tmp_Ar_i_j * tmp_cur_node_value * h_x*h_y;
- local_Ar_Ar_inner_product_sum += rho * pow (tmp_Ar_i_j, 2) * h_x*h_y;
- local_r_norm += rho * pow(tmp_cur_node_value, 2) * h_x*h_y;
- }
- }
- //cout << "[DEBUG]: Local"<< local_Ar_r_inner_product_sum << endl;
- MPI_Allreduce(&local_Ar_r_inner_product_sum, &global_Ar_r_inner_product_sum, 1, MPI_DOUBLE, MPI_SUM,
- comm);
- //cout << "[DEBUG]: "<< global_Ar_r_inner_product_sum << endl;
- MPI_Allreduce(&local_Ar_Ar_inner_product_sum, &global_Ar_Ar_inner_product_sum, 1, MPI_DOUBLE, MPI_SUM,
- MPI_COMM_WORLD);
- //cout << "[DEBUG]: "<< global_Ar_Ar_inner_product_sum << endl;
- double global_tau = global_Ar_r_inner_product_sum/ global_Ar_Ar_inner_product_sum;
- MPI_Allreduce(&local_r_norm, &global_r_norm, 1, MPI_DOUBLE, MPI_SUM,
- MPI_COMM_WORLD);
- double eps = 1e-06;
- if (global_r_norm < 0){
- cout << "[ERROR]: bad global r norm" << endl;
- }
- current_norm = fabs(global_tau)*sqrt(global_r_norm);
- //if (my_rank == 0)
- // cout << "[DEBUG]: solution norm "<< tmp_norm << endl;
- if (current_norm <= eps){
- if (my_rank == 0)
- cout << "[INFO]: r norm is " << sqrt(global_r_norm) << endl;
- should_i_stop = true;
- }else{
- should_i_stop = false;
- }
- return global_tau;
- }
- void ComputeNewW(double tau){
- //#pragma omp parallel for
- for(int i = 0; i < cur_block_size_x; ++i){
- for(int j = 0; j < cur_block_size_y; ++j){
- internal_data[ j + cur_block_size_y*i ] = old_internal_data[ j + cur_block_size_y*i ] - tau * internal_data[ j + cur_block_size_y*i ];
- //old_internal_data[ j + cur_block_size_y*i ] = 0.0;
- }
- }
- }
- int GetGlobalX(int i){
- return cur_block_global_offset_x + i;
- }
- int GetGlobalY(int j){
- return cur_block_global_offset_y + j;
- }
- public:
- MPIComputations(int inpM, int inpN, int inpGx, int inpGy, int inpx_left, int inpx_right, int inpy_bottom, int inpy_top, int inp_cur_block_global_coors_x, int inp_cur_block_global_coors_y, int inprank, MPI_Comm inpcomm){
- M = inpM;
- N = inpN;
- Gx = inpGx;
- Gy = inpGy;
- x_left = inpx_left;
- x_right = inpx_right;
- y_bottom = inpy_bottom;
- y_top = inpy_top;
- h_x = double((x_right - x_left)) / M;
- h_y = double((y_top - y_bottom)) / N;
- my_rank = inprank;
- comm = inpcomm;
- ofstream myfile;
- stringstream ss;
- ss<< "./file" << inp_cur_block_global_coors_x << "_" << inp_cur_block_global_coors_y << ".res";
- string filename;
- ss >> filename;
- my_file.open (filename);
- cur_block_global_coors_x = inp_cur_block_global_coors_x;
- cur_block_global_coors_y = inp_cur_block_global_coors_y;
- global_block_size_x = (M + 1) / Gx;
- global_block_size_y = (N + 1) / Gy;
- cur_block_size_x = global_block_size_x;
- cur_block_size_y = global_block_size_y;
- cur_block_global_offset_x = global_block_size_x * cur_block_global_coors_x;
- cur_block_global_offset_y = global_block_size_y * cur_block_global_coors_y;
- if ((cur_block_global_offset_x + 2 * global_block_size_x) > (M + 1)){
- cout << "I am "<< cur_block_global_coors_x<< " " << cur_block_global_coors_y << ". My cur xblock size is too big" << endl;
- cur_block_size_x += (M + 1) % Gx;
- }
- if ((cur_block_global_offset_y + 2 * global_block_size_y) > (N + 1)){
- cout << "I am "<< cur_block_global_coors_x<< " " << cur_block_global_coors_y << ". My cur yblock size is too big" << endl;
- cur_block_size_y += (N + 1) % Gy;
- }
- cout << "I am "<< cur_block_global_coors_x<< " " << cur_block_global_coors_y << ". My cur block size is " << cur_block_size_x << " " << cur_block_size_y << endl;
- internal_data.resize(cur_block_size_x * cur_block_size_y);
- old_internal_data.resize(cur_block_size_x * cur_block_size_y);
- //OX
- external_data[0].resize(cur_block_size_y);
- external_data[1].resize(cur_block_size_y);
- //OY
- external_data[2].resize(cur_block_size_x);
- external_data[3].resize(cur_block_size_x);
- }
- double Get(int i, int j) {
- return GetLocalIndex(i - cur_block_global_offset_x, j - cur_block_global_offset_y);
- }
- void Set(int i, int j, double v) {
- return SetLocalIndex(i - cur_block_global_offset_x, j - cur_block_global_offset_y, v);
- }
- void SyncMPI(){
- /*if(my_rank == 0)
- cout << "[INFO]: Sync started"<< endl;
- */
- MPI_Barrier(MPI_COMM_WORLD);
- //left and right sides
- //#pragma omp parallel for
- for(int j = 0; j < cur_block_size_y; ++j){
- external_data[ 0 ][ j ] = GetLocalIndex(0, j);//internal_data[ j ];
- external_data[ 1 ][ j ] = GetLocalIndex(cur_block_size_x - 1, j);//internal_data[ j + cur_block_size_y * (cur_block_size_x - 1) ];
- }
- //bottom and top sides
- //#pragma omp parallel for
- for(int i = 0; i < cur_block_size_x; ++i){
- external_data[ 2 ][ i ] = GetLocalIndex(i, 0);//internal_data[ cur_block_size_y*i ];
- external_data[ 3 ][ i ] = GetLocalIndex(i, cur_block_size_y - 1); //internal_data[ (cur_block_size_y - 1) + cur_block_size_y*i ];
- }
- int my_coords[2];
- int targets_ranks[4];
- MPI_Cart_coords(comm, my_rank, 2, my_coords);
- int neighbour_offsets[ 4 ][ 2 ] = {
- { -1, 0 },{ 1, 0 },
- { 0, -1 },{ 0, 1 }
- };
- //#pragma omp parallel for
- for(int i = 0; i < 4; i++){
- int target_coords[2];
- target_coords[ 0 ] = my_coords[ 0 ] + neighbour_offsets[ i ][ 0 ];
- target_coords[ 1 ] = my_coords[ 1 ] + neighbour_offsets[ i ][ 1 ];
- if (target_coords[0] >= 0 && target_coords[0] < Gx && target_coords[1] >= 0 && target_coords[1] < Gy){
- MPI_Cart_rank(comm, target_coords, &targets_ranks[ i ]);
- }
- else{
- targets_ranks[i] = -1;
- }
- }
- //Now we have rank for all targets
- for(int axis = 0; axis < 2; axis++){
- int parity_bit = (my_coords[ axis ]) % 2;
- //if parity_bit == 0 then
- // zuerst mit links, dann mit rechts tauschen
- //elif parity_bit == 1:
- // zuerst mit rechts, dann mit links tauschen
- for(int tmp = 0; tmp < 2; tmp++){
- parity_bit = 1 - parity_bit;
- //target id in external_data und targets_ranks
- int target_idx = 2 * axis + (1 - parity_bit);
- if (targets_ranks[target_idx] != -1){
- // вычисляем теги отправки и приема
- // в них зашиты номер ноды, ось, направление
- int send_tag = 100000 + my_rank * 100 + axis * 10 + parity_bit;
- int recv_tag = 100000 + targets_ranks[ target_idx ] * 100 + axis * 10 + (1-parity_bit);
- MPI_Status tmp_status;
- // если отправка не на себя, то отправляем
- if(my_rank != targets_ranks[ target_idx ]){
- MPI_Sendrecv_replace(&external_data[ target_idx ][ 0 ], external_data[ target_idx ].size(),
- MPI_DOUBLE, targets_ranks[ target_idx ], send_tag, targets_ranks[ target_idx ], recv_tag,
- comm, &tmp_status);
- }
- }
- }
- }
- MPI_Barrier(MPI_COMM_WORLD);
- }
- void DoIteration(bool &should_i_stop){
- ComputeMatrixR();
- //Now R Matrix is in internal_data
- double tau = ComputeTauAndStopCase(should_i_stop);
- //We have in our block jetzt:
- //in internal_data: R Matrix
- //in old_internal_data: w aus letzte iteration
- //and we have tau
- //jetzt koennen wir naechste w finden
- ComputeNewW(tau);
- SyncMPI();
- }
- double GetLocalIndex(int i, int j){
- //internal data
- if ((j >= 0) && (j < cur_block_size_y) && (i >= 0) && (i < cur_block_size_x)){
- return internal_data[ j + cur_block_size_y*i ];
- }
- //external data
- //OX
- if((j >= 0) && (j < cur_block_size_y)){
- if (i == -1)
- return external_data[ 0 ][ j ];
- if (i == cur_block_size_x)
- return external_data[ 1 ][ j ];
- }
- //OY
- if((i >= 0) && (i < cur_block_size_x)){
- if (j == -1)
- return external_data[ 2 ][ i ];
- if (j == cur_block_size_y)
- return external_data[ 3 ][ i ];
- }
- cout << "[ERROR]: bad local index" << endl;
- return nan("");
- }
- void SetLocalIndex(int i, int j, double v){
- if ((j >= 0) && (j < cur_block_size_y) && (i >= 0) && (i < cur_block_size_x)){
- internal_data[ j + cur_block_size_y*i ] = v;
- }else{
- cout << "[ERROR]: trying to set data outside the local area" << endl;
- }
- }
- double CompareWithExact(bool print_data_to_files) {
- double local_diff_norm = 0.0;
- double global_diff_norm = 0.0;
- /*
- if (my_rank == 0)
- cout << "[INFO]: Starting computing compare with exact" << endl;
- */
- vector <double> tmp_elements(cur_block_size_x * cur_block_size_y, 0.0);
- //#pragma omp parallel for
- for (int i = 0; i < cur_block_size_x; ++i){
- for (int j = 0; j < cur_block_size_y; ++j){
- int current_node_global_offset_x = GetGlobalX(i),
- current_node_global_offset_y = GetGlobalY(j);
- int glob_x = current_node_global_offset_x,
- glob_y = current_node_global_offset_y;
- double tmp_elem = Get(glob_x, glob_y) - GetNodeFromExact(glob_x, glob_y);
- //local_diff_norm = max( fabs(tmp_elem), local_diff_norm);
- //local_diff_norm += rho * pow(tmp_elem, 2) * h_x * h_y;
- tmp_elements [j + cur_block_size_y*i] = tmp_elem;
- }
- }
- //cout << "[INFO]: local max diff in " << cur_block_global_offset_x << " " << cur_block_global_offset_y << " " << local_diff_norm << endl;
- for (int i = 0; i < cur_block_size_x; ++i){
- for (int j = 0; j < cur_block_size_y; ++j){
- double rho = 1.0;
- int current_node_global_offset_x = GetGlobalX(i),
- current_node_global_offset_y = GetGlobalY(j);
- int glob_x = current_node_global_offset_x,
- glob_y = current_node_global_offset_y;
- switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
- case 2:
- case 4:
- case 6:
- case 8:
- //angle
- rho = 0.25;
- break;
- case 1:
- case 3:
- case 5:
- case 7:
- //side
- rho = 0.5;
- break;
- case 0:
- //internal
- rho = 1.0;
- break;
- default:
- cout << "[ERROR]: Bad global coords compute exact. Local:" << i << " " << j << ". Global:" << glob_x << " " << glob_y <<endl;
- }
- double tmp_elem = tmp_elements [j + cur_block_size_y*i];
- local_diff_norm = max( fabs(tmp_elem), local_diff_norm);
- //local_diff_norm += rho * pow(tmp_elem, 2) * h_x * h_y;
- }
- }
- MPI_Allreduce(&local_diff_norm, &global_diff_norm, 1, MPI_DOUBLE, MPI_MAX,
- MPI_COMM_WORLD);
- //global_diff_norm = sqrt(global_diff_norm);
- if (print_data_to_files){
- for (int i = 0; i < cur_block_size_x; ++i){
- for (int j = 0; j < cur_block_size_y; ++j){
- int current_node_global_offset_x = GetGlobalX(i),
- current_node_global_offset_y = GetGlobalY(j);
- int glob_x = current_node_global_offset_x,
- glob_y = current_node_global_offset_y;
- my_file << Get(glob_x, glob_y) << ":" << GetNodeFromExact(glob_x, glob_y) << "\t";
- }
- my_file << endl;
- }
- }
- return global_diff_norm;
- }
- double GetCurrentNorm(){
- return current_norm;
- }
- };
- int main(int argc, char* argv[]){
- const double x_left = -2, x_right = 3;
- const double y_bottom = -1, y_top = 4;
- double time_start, time_stop;
- int N, Gx, Gy;
- int dim[2], period[2];
- MPI_Comm comm;
- //N - global grid size
- N = atoi(argv[1]);
- //Gx
- Gx = dim[0] = atoi(argv[2]);
- //Gy
- Gy = dim[1] = atoi(argv[3]);
- bool print_data_to_files = false;
- if (argc == 4){
- print_data_to_files = true;
- }
- period[0]=0;
- period[1]=0;
- MPI_Init(&argc, &argv);
- time_start = MPI_Wtime();
- int world_size;
- int my_rank;
- MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
- MPI_Comm_size(MPI_COMM_WORLD, &world_size);
- if(my_rank == 0){
- cout << "[INFO]: N = " << N << endl;
- cout << "[INFO]: Gx = " << dim[0] << endl;
- cout << "[INFO]: Gy = " << dim[1] << endl;
- }
- if(argc <= 3){
- if(my_rank == 0)
- cout << "[ERROR]: Usage: mpieval <N> <Gx> <Gy> [<print_data_to_files>]" << endl;
- MPI_Abort(MPI_COMM_WORLD, 1);
- }
- if(Gx * Gy != world_size){
- if(my_rank == 0)
- cout << "[ERROR]: mpi world size is not equal to "<< Gx << "*" << Gy << endl;
- MPI_Abort(MPI_COMM_WORLD, 1);
- }
- MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, 1, &comm);
- /*if(my_rank == 0)
- cout << "[INFO]: Cart created"<<endl;
- */
- MPI_Comm_rank(comm, &my_rank);
- int my_coords[2];
- MPI_Cart_coords(comm, my_rank, 2, my_coords);
- class MPIComputations w_func(N, N, Gx, Gy, x_left, x_right, y_bottom, y_top, my_coords[0], my_coords[1], my_rank, comm);
- int iteration_num = 0;
- bool should_i_stop = false;
- while ( should_i_stop != true ){
- w_func.DoIteration(should_i_stop);
- if ( (my_rank == 0) && (iteration_num % 10000 == 0) ){
- cout << "[INFO]: Iteration " << iteration_num << " difference with previous is: "<< w_func.GetCurrentNorm() << endl;
- }
- iteration_num++;
- }
- MPI_Barrier(MPI_COMM_WORLD);
- //vector <double> local_elements;
- double comparing = w_func.CompareWithExact(print_data_to_files);
- if (my_rank == 0)
- cout << "[INFO]: Diff with exact solution: " << comparing << endl;
- //myfile.close();
- time_stop = MPI_Wtime();
- if( my_rank == 0 )
- cout << "Finished!" << endl
- << "Total iterations: " << iteration_num << endl
- << "Elapsed time: " << (time_stop - time_start) << endl
- << "Time per iteration: " << (time_stop - time_start) / double(iteration_num) << endl;
- MPI_Finalize();
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement