SHARE
TWEET

Untitled

manac68974 Dec 3rd, 2019 90 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <mpi.h>
  2.  
  3. #include <iomanip>
  4. #include <iostream>
  5. #include <vector>
  6. #include <map>
  7. #include <cmath>
  8. #include <cstdlib>
  9. #include <algorithm>
  10.  
  11. using namespace std;
  12.  
  13. class MPIComputations{
  14.     int M, N;
  15.     int cur_block_size_x, cur_block_size_y, cur_block_global_offset_x, cur_block_global_offset_y;
  16.     int Gx, Gy;
  17.     int cur_block_global_coors_x, cur_block_global_coors_y;
  18.     double h_x, h_y;
  19.     int global_block_size_x, global_block_size_y;
  20.     int my_rank;
  21.     MPI_Comm comm;
  22.     double x_left, x_right, y_bottom, y_top;
  23.  
  24.     vector<double> internal_data;
  25.     vector<double> old_internal_data;
  26.     vector<double> external_data[4];
  27.  
  28.     int IsNodeInternalCornerOrSide(int current_node_global_offset_x, int current_node_global_offset_y){
  29.  
  30.         //corners
  31.         //left bottom corner
  32.         if (current_node_global_offset_x == 0 && current_node_global_offset_y == 0){
  33.             return 2;
  34.         }
  35.  
  36.         //left top corner
  37.         if (current_node_global_offset_x == 0 && current_node_global_offset_y == N){
  38.             return 4;
  39.         }
  40.  
  41.         //right bottom corner
  42.         if (current_node_global_offset_x == N && current_node_global_offset_y == 0){
  43.             return 6;
  44.         }
  45.  
  46.         //right top corner
  47.         if (current_node_global_offset_x == N && current_node_global_offset_y == N){
  48.             return 8;
  49.         }
  50.  
  51.         //sides
  52.         //left side
  53.         if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
  54.             current_node_global_offset_x == 0){
  55.             return 1;
  56.         }
  57.  
  58.         //right side
  59.         if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
  60.             current_node_global_offset_x == N){
  61.             return 3;
  62.         }
  63.  
  64.         //bottom side
  65.         if (current_node_global_offset_x >= 1 && current_node_global_offset_x <= N - 1 &&
  66.             current_node_global_offset_y == 0){
  67.             return 5;
  68.         }
  69.  
  70.         //top side
  71.         if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= N - 1 &&
  72.             current_node_global_offset_y == N)){
  73.             return 7;
  74.         }
  75.  
  76.         //internal
  77.         if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= N - 1) &&
  78.             (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1)){
  79.             return 0;
  80.         }
  81.  
  82.         return -1;
  83.     }
  84.  
  85.     double k(double x, double y) {
  86.         return 1 + pow(x + y, 2);
  87.     }
  88.  
  89.     double q(double x, double y) {
  90.         return 1;
  91.     }
  92.  
  93.     double u(double x, double y) {
  94.         return 2.0 / (1 + pow(x, 2) + pow(y, 2));
  95.     }
  96.  
  97.     // psi_R(x, y) = k(A2, y) * du/dx(A2, y)
  98.     double psi_R(double y) {
  99.         return (-12) * (pow((y + 3), 2) + 1) / pow((pow(y, 2) + 10), 2);
  100.     }
  101.  
  102.     // psi_L(x, y) = -k(A1, y) * du/dx(A1, y)
  103.     double psi_L(double y) {
  104.         return (-8) * (pow((y - 2), 2) + 1) / pow((pow(y, 2) + 5), 2);
  105.     }
  106.  
  107.     // psi_T(x, y) = k(x, B2) * du/dy(x, B2)
  108.     double psi_T(double x) {
  109.         return (-16) * (pow((x + 4), 2) + 1) / pow((pow(x, 2) + 17), 2);
  110.     }
  111.  
  112.     // psi_B(x, y) = -k(x, B1) * du/dy(x, B1)
  113.     double psi_B(double x) {
  114.         return (-4) * (pow((x - 1), 2) + 1) / pow((pow(x, 2) + 2), 2);
  115.     }
  116.  
  117.     // right-part function of Poisson equation
  118.     double F(double x, double y) {
  119.         return 2 * (pow(x,4) + pow(y,4) + 2 * (pow(x,2) + 3) * pow(y,2) + 6 * pow(x,2) + 16*x*y + 5)
  120.         / pow((1 + pow(x, 2) + pow(y, 2)), 3);
  121.     }
  122.  
  123.     //inner_product(A[i], internal_data)
  124.     double ComputeMagicInnerProductA_iw (int current_node_global_offset_x, int current_node_global_offset_y){
  125.  
  126.         int glob_x = current_node_global_offset_x;
  127.         int glob_y = current_node_global_offset_y;
  128.  
  129.         double result = 0.0;
  130.  
  131.         map <string,bool> neighbours = {
  132.             {"left", true},
  133.             {"right", true},
  134.             {"bottom", true},
  135.             {"top", true}
  136.         };
  137.  
  138.         double left_neighbour = 0.0, right_neighbour = 0.0, bottom_neighbour = 0.0, top_neighbour = 0.0, this_node = 0.0;
  139.         double left_coeff = 1.0, right_coeff = 1.0, bottom_coeff = 1.0, top_coeff = 1.0, this_coeff = 1.0;
  140.  
  141.         switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
  142.             case 2:
  143.                 //left bottom corner
  144.                 neighbours["left"] = false;
  145.                 neighbours["bottom"] = false;
  146.                 break;
  147.             case 4:
  148.                 //left top corner
  149.                 neighbours["left"] = false;
  150.                 neighbours["top"] = false;
  151.                 break;
  152.             case 6:
  153.                 //right bottom corner
  154.                 neighbours["right"] = false;
  155.                 neighbours["bottom"] = false;
  156.                 break;
  157.             case 8:
  158.                 //right top corner
  159.                 neighbours["right"] = false;
  160.                 neighbours["top"] = false;
  161.                 break;
  162.             case 1:
  163.                 //left side
  164.                 neighbours["left"] = false;
  165.                 break;
  166.             case 3:
  167.                 //right side
  168.                 neighbours["right"] = false;
  169.                 break;
  170.             case 5:
  171.                 //bottom side
  172.                 neighbours["bottom"] = false;
  173.                 break;
  174.             case 7:
  175.                 //top side
  176.                 neighbours["top"] = false;
  177.                 break;
  178.             case 0:
  179.                 //internal
  180.                 break;
  181.             default:
  182.                 cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y<<endl;
  183.         }
  184.  
  185.         if (!neighbours["left"]){
  186.             right_coeff = 2.0;
  187.             //left_coeff = 0.0;
  188.         }
  189.  
  190.         if (!neighbours["right"]){
  191.             left_coeff = 2.0;
  192.             //right_coeff = 0.0;
  193.         }
  194.  
  195.         if (!neighbours["bottom"]){
  196.             top_coeff = 2.0;
  197.             //bottom_coeff = 0.0;
  198.         }
  199.  
  200.         if (!neighbours["top"]){
  201.             bottom_coeff = 2.0;
  202.             //top_coeff = 0.0;
  203.         }
  204.  
  205.  
  206.  
  207.  
  208.  
  209.         if (neighbours["left"]){
  210.             left_coeff *= -k(x_left + (glob_x - 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
  211.             left_neighbour = Get(glob_x - 1, glob_y);
  212.         }
  213.  
  214.         if (neighbours["right"]){
  215.             right_coeff *= -k(x_left + (glob_x + 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
  216.             right_neighbour = Get(glob_x + 1, glob_y);
  217.         }
  218.  
  219.         if (neighbours["bottom"]){
  220.             bottom_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y - 0.5) * h_y) / pow(h_y, 2);
  221.             bottom_neighbour = Get(glob_x, glob_y - 1);
  222.         }
  223.  
  224.         if (neighbours["top"]){
  225.             top_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y + 0.5) * h_y) / pow(h_y, 2);
  226.             top_neighbour = Get(glob_x, glob_y + 1);
  227.         }
  228.  
  229.         this_coeff = q(x_left + glob_x * h_x, y_bottom + glob_y * h_y) - left_neighbour - right_neighbour - bottom_neighbour - top_neighbour;
  230.         this_node = Get(glob_x, glob_y);
  231.  
  232.         result = left_coeff * left_neighbour +
  233.                  right_coeff * right_neighbour +
  234.                  bottom_coeff * bottom_neighbour +
  235.                  top_coeff * top_neighbour +
  236.                  this_coeff * this_node;
  237.  
  238.         return result;
  239.     }
  240.  
  241.     double GetNodeFromB(int current_node_global_offset_x, int current_node_global_offset_y) {
  242.  
  243.         int glob_x = current_node_global_offset_x;
  244.         int glob_y = current_node_global_offset_y;
  245.  
  246.         double result = 0.0;
  247.  
  248.         switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
  249.             case 2:
  250.                 //left bottom corner
  251.                 result = F(x_left, y_bottom) + 2.0 / h_x * psi_L(y_bottom) + 2.0 / h_y * psi_B(x_left);
  252.                 break;
  253.             case 4:
  254.                 //left top corner
  255.                 result = F(x_left, y_top) + 2.0 / h_x * psi_L(y_top) + 2.0 / h_y * psi_T(x_left);
  256.                 break;
  257.             case 6:
  258.                 //right bottom corner
  259.                 result = F(x_right, y_bottom) + 2.0 / h_x * psi_R(y_bottom) + 2.0 / h_y * psi_B(x_right);
  260.                 break;
  261.             case 8:
  262.                 //right top corner
  263.                 result = F(x_right, y_top) + 2.0 / h_x * psi_R(y_top) + 2.0 / h_y * psi_T(x_right);
  264.                 break;
  265.             case 1:
  266.                 //left side
  267.                 result = F(x_left, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_L(y_bottom + glob_y * h_y);
  268.                 break;
  269.             case 3:
  270.                 //right side
  271.                 result = F(x_right, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_R(y_bottom + glob_y * h_y);
  272.                 break;
  273.             case 5:
  274.                 //bottom side
  275.                 result = F(x_left + glob_x * h_x, y_bottom) + 2.0 / h_y * psi_B(x_left + glob_x * h_x);
  276.                 break;
  277.             case 7:
  278.                 //top side
  279.                 result = F(x_left + glob_x * h_x, y_top) + 2.0 / h_y * psi_T(x_left + glob_x * h_x);
  280.                 break;
  281.             case 0:
  282.                 //internal
  283.                 result = F(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
  284.                 break;
  285.             default:
  286.                 cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y <<endl;
  287.  
  288.         }
  289.  
  290.         return result;
  291.  
  292.     }
  293.  
  294.     double GetNodeFromExact(int current_node_global_offset_x, int current_node_global_offset_y) {
  295.  
  296.         int glob_x = current_node_global_offset_x;
  297.         int glob_y = current_node_global_offset_y;
  298.  
  299.         return u(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
  300.  
  301.     }
  302.  
  303.  
  304.  
  305.     void ComputeMatrixR(){
  306.         if(my_rank == 0)
  307.             cout << "[INFO]: Computation of matrix r started"<<endl;
  308.  
  309.         double loc_t1, loc_t2;
  310.  
  311.         vector<double> r_tmp_matrix (cur_block_size_x*cur_block_size_y, 0.0);
  312.  
  313.         loc_t1 = MPI_Wtime();
  314.  
  315.         for(int i = 0; i < cur_block_size_x; ++i){
  316.             for(int j = 0; j < cur_block_size_y; ++j){
  317.  
  318.                 int current_node_global_offset_x = GetGlobalX(i),
  319.                     current_node_global_offset_y = GetGlobalY(j);
  320.  
  321.                 int glob_x = current_node_global_offset_x,
  322.                     glob_y = current_node_global_offset_y;
  323.  
  324.                 r_tmp_matrix [ j + cur_block_size_y*i ] = ComputeMagicInnerProductA_iw(glob_x,glob_y) - GetNodeFromB(glob_x, glob_y);
  325.  
  326.             }
  327.         }
  328.  
  329.  
  330.         loc_t2 = MPI_Wtime();
  331.  
  332.         if(my_rank == 0)
  333.             cout << "[INFO]: Computing r matrix time: " << loc_t2 - loc_t1 << endl;
  334.  
  335.         loc_t1 = MPI_Wtime();
  336.  
  337.         MPI_Barrier(MPI_COMM_WORLD);
  338.  
  339.         loc_t2 = MPI_Wtime();
  340.         /*if(my_rank == 0)
  341.             cout << "[INFO]: Barrier jumping time: " << loc_t2 - loc_t1 << endl;
  342.         */
  343.  
  344.         for(int i = 0; i < cur_block_size_x; ++i){
  345.             for(int j = 0; j < cur_block_size_y; ++j){
  346.  
  347.                 old_internal_data[ j + cur_block_size_y*i ] = internal_data[ j + cur_block_size_y*i];
  348.                 internal_data[j + cur_block_size_y*i ] = r_tmp_matrix[j + cur_block_size_y*i];
  349.  
  350.             }
  351.         }
  352.  
  353.         SyncMPI();
  354.  
  355.     }
  356.  
  357.     double ComputeTauAndStopCase(bool &should_i_stop){
  358.  
  359.         double local_Ar_r_inner_product_sum = 0.0;
  360.         double local_Ar_Ar_inner_product_sum = 0.0;
  361.         double global_Ar_r_inner_product_sum = 0.0;
  362.         double global_Ar_Ar_inner_product_sum = 0.0;
  363.         double local_r_norm = 0.0;
  364.         double global_r_norm = 0.0;
  365.  
  366.  
  367.         for(int i = 0; i < cur_block_size_x; ++i){
  368.             for(int j = 0; j < cur_block_size_y; ++j){
  369.                 double rho = 1.0;
  370.                
  371.                 int current_node_global_offset_x = GetGlobalX(i),
  372.                     current_node_global_offset_y = GetGlobalY(j);
  373.  
  374.                 int glob_x = current_node_global_offset_x,
  375.                     glob_y = current_node_global_offset_y;
  376.  
  377.                 double tmp_Ar_i_j = ComputeMagicInnerProductA_iw(glob_x, glob_y);
  378.  
  379.                 switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
  380.                     case 2:
  381.                     case 4:
  382.                     case 6:
  383.                     case 8:
  384.                         //angle
  385.                         rho = 0.25;
  386.                         break;
  387.                     case 1:
  388.                     case 3:
  389.                     case 5:
  390.                     case 7:
  391.                         //side
  392.                         rho = 0.5;
  393.                         break;
  394.                     case 0:
  395.                         //internal
  396.                         rho = 1.0;
  397.                         break;
  398.                     default:
  399.                         cout << "[ERROR]: Bad global coords compute tau. Global:" << glob_x << " " << glob_y << endl;
  400.                 }
  401.  
  402.                 double tmp_cur_node_value = Get(glob_x, glob_y);
  403.  
  404.                 local_Ar_r_inner_product_sum += rho * tmp_Ar_i_j * tmp_cur_node_value * h_x*h_y;
  405.                 local_Ar_Ar_inner_product_sum += rho * pow (tmp_Ar_i_j, 2) * h_x*h_y;
  406.                 local_r_norm += rho * pow(tmp_cur_node_value, 2) * h_x*h_y;
  407.  
  408.             }
  409.         }
  410.  
  411.         //cout << "[DEBUG]: Local"<< local_Ar_r_inner_product_sum << endl;
  412.  
  413.         MPI_Allreduce(&local_Ar_r_inner_product_sum, &global_Ar_r_inner_product_sum, 1, MPI_DOUBLE, MPI_SUM,
  414.             comm);
  415.  
  416.  
  417.         //cout << "[DEBUG]: "<< global_Ar_r_inner_product_sum << endl;
  418.  
  419.         MPI_Allreduce(&local_Ar_Ar_inner_product_sum, &global_Ar_Ar_inner_product_sum, 1, MPI_DOUBLE, MPI_SUM,
  420.             MPI_COMM_WORLD);
  421.  
  422.         //cout << "[DEBUG]: "<< global_Ar_Ar_inner_product_sum << endl;
  423.  
  424.         double global_tau = global_Ar_r_inner_product_sum/ global_Ar_Ar_inner_product_sum;
  425.  
  426.         MPI_Allreduce(&local_r_norm, &global_r_norm, 1, MPI_DOUBLE, MPI_SUM,
  427.             MPI_COMM_WORLD);
  428.  
  429.         double eps = 1e-06;
  430.  
  431.         if (global_r_norm < 0){
  432.             cout << "[ERROR]: bad global r norm" << endl;
  433.         }
  434.  
  435.         double tmp_norm = fabs(global_tau)*sqrt(global_r_norm);
  436.  
  437.         if (my_rank == 0)
  438.             cout << "[DEBUG]: solution norm "<< tmp_norm << endl;
  439.  
  440.         if (tmp_norm <= eps){
  441.             should_i_stop = true;
  442.         }else{
  443.             should_i_stop = false;
  444.         }
  445.  
  446.  
  447.         return global_tau;
  448.  
  449.     }
  450.  
  451.  
  452.     void ComputeNewW(double tau){
  453.  
  454.         for(int i = 0; i < cur_block_size_x; ++i){
  455.             for(int j = 0; j < cur_block_size_y; ++j){
  456.                 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 ];
  457.                 old_internal_data[j + cur_block_size_y*i ] = 0.0;
  458.             }
  459.         }
  460.     }
  461.  
  462.     int GetGlobalX(int i){
  463.         return cur_block_global_offset_x + i;
  464.     }
  465.  
  466.     int GetGlobalY(int j){
  467.         return cur_block_global_offset_y + j;
  468.     }
  469. public:
  470.     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){
  471.  
  472.         M = inpM;
  473.         N = inpN;
  474.  
  475.         Gx = inpGx;
  476.         Gy = inpGy;
  477.  
  478.         x_left = inpx_left;
  479.         x_right = inpx_right;
  480.         y_bottom = inpy_bottom;
  481.         y_top = inpy_top;
  482.  
  483.         h_x = double((x_right - x_left)) / N;
  484.         h_y = double((y_top - y_bottom)) / N;
  485.  
  486.         my_rank = inprank;
  487.         comm = inpcomm;
  488.  
  489.         cur_block_global_coors_x = inp_cur_block_global_coors_x;
  490.         cur_block_global_coors_y = inp_cur_block_global_coors_y;
  491.  
  492.         global_block_size_x = (N + 1) / Gx;
  493.         global_block_size_y = (N + 1) / Gy;
  494.  
  495.         cur_block_size_x = global_block_size_x;
  496.         cur_block_size_y = global_block_size_y;
  497.        
  498.         cur_block_global_offset_x = global_block_size_x * cur_block_global_coors_x;
  499.         cur_block_global_offset_y = global_block_size_y * cur_block_global_coors_y;
  500.  
  501.         if (cur_block_global_offset_x + global_block_size_x > N + 1){
  502.             cur_block_size_x = (N + 1) % Gx;
  503.         }
  504.  
  505.         if (cur_block_global_offset_y + global_block_size_y > N + 1){
  506.             cur_block_size_y = (N + 1) % Gy;
  507.         }
  508.  
  509.         //cout << "[INFO]: For block" << my_rank << " size " << cur_block_size_x << " " << cur_block_size_y << "offset"<<cur_block_global_offset_x<<" " <<cur_block_global_offset_y <<endl;
  510.  
  511.  
  512.         //vector< vector<double> > tmp_internal_data(cur_block_size_x, vector<double>(cur_block_size_y, 0.0));
  513.         //internal_data = tmp_internal_data;
  514.  
  515.         //vector< vector<double> > tmp_old_internal_data(cur_block_size_x, vector<double>(cur_block_size_y, 0.0));
  516.         //old_internal_data = tmp_old_internal_data;
  517.  
  518.         internal_data.resize(cur_block_size_x * cur_block_size_y);
  519.         old_internal_data.resize(cur_block_size_x * cur_block_size_y);
  520.  
  521.  
  522.         //OX
  523.         external_data[0].resize(cur_block_size_y);
  524.         external_data[1].resize(cur_block_size_y);
  525.  
  526.         //OY
  527.         external_data[2].resize(cur_block_size_x);
  528.         external_data[3].resize(cur_block_size_x);
  529.  
  530.     }
  531.  
  532.     double Get(int i, int j) {
  533.         return GetLocalIndex(i - cur_block_global_offset_x, j - cur_block_global_offset_y);
  534.     }
  535.  
  536.     void Set(int i, int j, double v) {
  537.         return SetLocalIndex(i - cur_block_global_offset_x, j - cur_block_global_offset_y, v);
  538.     }
  539.  
  540.     void SyncMPI(){
  541.  
  542.         double loc_t1, loc_t2;
  543.         /*if(my_rank == 0)
  544.             cout << "[INFO]: Sync started"<< endl;
  545.         */
  546.         //left and right sides
  547.  
  548.  
  549.         loc_t1 = MPI_Wtime();
  550.  
  551.         MPI_Barrier(MPI_COMM_WORLD);
  552.  
  553.         loc_t2 = MPI_Wtime();
  554.  
  555.        
  556.         /*
  557.         if(my_rank == 0)
  558.             cout << "[INFO]: Barrier jumping time: " << loc_t2 - loc_t1 << endl;
  559.         */
  560.  
  561.         for(int j = 0; j < cur_block_size_y; ++j){
  562.  
  563.             external_data[ 0 ][ j ] = GetLocalIndex(0,j);//internal_data[ j ];
  564.  
  565.             external_data[ 1 ][ j ] = GetLocalIndex(cur_block_size_x - 1,j);//internal_data[ j + cur_block_size_y * (cur_block_size_x - 1) ];
  566.  
  567.         }
  568.  
  569.         //bottom and top sides
  570.         for(int i = 0; i < cur_block_size_x; ++i){
  571.  
  572.             external_data[ 2 ][ i ] = GetLocalIndex(i,0);//internal_data[ cur_block_size_y*i ];
  573.  
  574.             external_data[ 3 ][ i ] = GetLocalIndex(i,cur_block_size_y - 1); //internal_data[ (cur_block_size_y - 1) + cur_block_size_y*i ];
  575.  
  576.         }
  577.  
  578.         int my_coords[2];
  579.         int targets_ranks[4];
  580.        
  581.         MPI_Cart_coords(comm, my_rank, 2, my_coords);
  582.  
  583.         int neighbour_offsets[ 4 ][ 2 ] = {
  584.             { -1, 0 },{ 1, 0 },
  585.             { 0, -1 },{ 0, 1 }
  586.         };
  587.  
  588.         for(int i = 0; i < 4; i++){
  589.  
  590.             int target_coords[2];
  591.  
  592.             target_coords[ 0 ] = my_coords[ 0 ] + neighbour_offsets[ i ][ 0 ];
  593.             target_coords[ 1 ] = my_coords[ 1 ] + neighbour_offsets[ i ][ 1 ];
  594.            
  595.             if (target_coords[0] >= 0 && target_coords[0] < Gx && target_coords[1] >= 0 && target_coords[1] < Gy){
  596.                 //cout << "[INFO]: Try to get target"<< target_coords[0] << " " << target_coords[1] << endl;
  597.                 MPI_Cart_rank(comm, target_coords, &targets_ranks[ i ]);
  598.                 //cout << "[INFO]: Target"<< target_coords[0] << " " << target_coords[1] << "getted!" << endl;
  599.                 //cout << "[INFO]: Target ranks: " << targets_ranks[i]<<endl;
  600.             }
  601.             else{
  602.                 targets_ranks[i] = -1;
  603.             }
  604.  
  605.         }
  606.  
  607.         //Now we have rank for all targets
  608.         /*
  609.         if(my_rank == 0)
  610.             cout << "[INFO]: Targets getted"<< endl;
  611.         */
  612.         for(int axis = 0; axis < 2; axis++){
  613.            
  614.             int parity_bit = (my_coords[ axis ]) % 2;
  615.  
  616.             //if parity_bit == 0 then
  617.             //  zuerst mit links, dann mit rechts tauschen
  618.             //elif parity_bit == 1:
  619.             //  zuerst mit rechts, dann mit links tauschen
  620.  
  621.             for(int tmp = 0; tmp < 2; tmp++){
  622.                 parity_bit = 1 - parity_bit;
  623.  
  624.                 //target id in external_data und targets_ranks
  625.                 int target_idx = 2 * axis + (1 - parity_bit);
  626.  
  627.                 if (targets_ranks[target_idx] != -1){
  628.  
  629.                     // вычисляем теги отправки и приема
  630.                     // в них зашиты номер ноды, ось, направление
  631.                     int send_tag = 100000 + my_rank * 100 + axis * 10 + parity_bit;
  632.                     int recv_tag = 100000 + targets_ranks[ target_idx ] * 100 + axis * 10 + (1-parity_bit);
  633.                    
  634.                     MPI_Status tmp_status;
  635.                     // если отправка не на себя, то отправляем
  636.                     if(my_rank != targets_ranks[ target_idx ]){
  637.  
  638.                         MPI_Sendrecv_replace(&external_data[ target_idx ][ 0 ], external_data[ target_idx ].size(),
  639.                             MPI_DOUBLE, targets_ranks[ target_idx ], send_tag, targets_ranks[ target_idx ], recv_tag,
  640.                             comm, &tmp_status);
  641.  
  642.                     }
  643.                 }
  644.             }            
  645.         }
  646.         /*
  647.         if(my_rank == 0)
  648.             cout << "[INFO]: Sync finished"<< endl;
  649.         */
  650.  
  651.  
  652.         loc_t1 = MPI_Wtime();
  653.  
  654.         MPI_Barrier(MPI_COMM_WORLD);
  655.  
  656.         loc_t2 = MPI_Wtime();
  657.  
  658.        
  659.         /*
  660.         if(my_rank == 0)
  661.             cout << "[INFO]: Barrier jumping time: " << loc_t2 - loc_t1 << endl;
  662.         */
  663.  
  664.     }
  665.  
  666.  
  667.  
  668.     void DoIteration(bool &should_i_stop){
  669.  
  670.         double loc_t1, loc_t2;
  671.  
  672.         double metrics = 0;
  673.  
  674.         ComputeMatrixR();
  675.  
  676.         //Now R Matrix is in internal_data
  677.  
  678.         double tau = ComputeTauAndStopCase(should_i_stop);
  679.  
  680.         //We have in our block jetzt:
  681.             //in internal_data:  R Matrix
  682.             //in old_internal_data: w aus letzte iteration
  683.         //and we have tau
  684.         //jetzt koennen wir naechste w finden
  685.  
  686.         ComputeNewW(tau);
  687.  
  688.  
  689.         SyncMPI();
  690.  
  691.     }
  692.  
  693.     double GetLocalIndex(int i, int j){
  694.         //internal data
  695.         if ((j >= 0) && (j < cur_block_size_y) && (i >= 0) && (i < cur_block_size_x)){
  696.             return internal_data[ j + cur_block_size_y*i ];
  697.         }
  698.  
  699.         //external data
  700.         //OX
  701.         if((j >= 0) && (j < cur_block_size_y)){
  702.  
  703.             if (i == -1)
  704.                 return external_data[ 0 ][ j ];
  705.  
  706.             if (i == cur_block_size_x)
  707.                 return external_data[ 1 ][ j ];
  708.  
  709.         }
  710.  
  711.         //OY
  712.         if((i >= 0) && (i < cur_block_size_x)){
  713.  
  714.             if (j == -1)
  715.                 return external_data[ 2 ][ i ];
  716.             if (j == cur_block_size_y)
  717.                 return external_data[ 3 ][ i ];
  718.  
  719.         }
  720.  
  721.         return nan("");
  722.     }
  723.  
  724.    
  725.     void SetLocalIndex(int i, int j, double v){
  726.        
  727.         if ((j >= 0) && (j < cur_block_size_y) && (i >= 0) && (i < cur_block_size_x)){
  728.             internal_data[ j + cur_block_size_y*i ] = v;
  729.         }else{
  730.             cout << "[ERROR]: trying to set data outside the local area" << endl;
  731.         }
  732.  
  733.     }
  734.  
  735.  
  736.  
  737.     double CompareWithExact() {
  738.  
  739.         double local_diff_norm = 0.0;
  740.         double global_diff_norm = 0.0;
  741.         /*
  742.         if (my_rank == 0)
  743.             cout << "[INFO]: Starting computing compare with exact" << endl;
  744.         */
  745.         for (int i = 0; i < cur_block_size_x; ++i){
  746.             for (int j = 0; j < cur_block_size_y; ++j){
  747.                
  748.  
  749.                 double rho = 1.0;
  750.                
  751.                 int current_node_global_offset_x = GetGlobalX(i),
  752.                     current_node_global_offset_y = GetGlobalY(j);
  753.  
  754.                 int glob_x = current_node_global_offset_x,
  755.                     glob_y = current_node_global_offset_y;
  756.  
  757.                 switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
  758.                     case 2:
  759.                     case 4:
  760.                     case 6:
  761.                     case 8:
  762.                         //angle
  763.                         rho = 0.25;
  764.                         break;
  765.                     case 1:
  766.                     case 3:
  767.                     case 5:
  768.                     case 7:
  769.                         //side
  770.                         rho = 0.5;
  771.                         break;
  772.                     case 0:
  773.                         //internal
  774.                         rho = 1.0;
  775.                         break;
  776.                     default:
  777.                         cout << "[ERROR]: Bad global coords compute exact. Local:" << i << " " << j << ". Global:" << glob_x << " " << glob_y <<endl;
  778.                 }
  779.  
  780.                 double tmp_elem = Get(glob_x, glob_y) - GetNodeFromExact(glob_x, glob_y);
  781.  
  782.                 //local_diff_norm = max( fabs(tmp_elem), local_diff_norm);
  783.  
  784.                 local_diff_norm += rho * pow(tmp_elem, 2) * h_x * h_y;
  785.  
  786.             }
  787.         }
  788.  
  789.         //cout << "[INFO]: local max diff in " << cur_block_global_offset_x << " " << cur_block_global_offset_y << " " << local_diff_norm << endl;
  790.  
  791.         MPI_Allreduce(&local_diff_norm, &global_diff_norm, 1, MPI_DOUBLE, MPI_SUM,
  792.             MPI_COMM_WORLD);
  793.  
  794.         global_diff_norm = sqrt(global_diff_norm);
  795.  
  796.         return global_diff_norm;
  797.     }
  798. };
  799.  
  800.  
  801.  
  802.  
  803. int main(int argc, char* argv[]){
  804.  
  805.     const double x_left = -2, x_right = 3;
  806.     const double y_bottom = -1, y_top = 4;
  807.  
  808.     double time_start, time_stop;
  809.     double loc_t1, loc_t2;
  810.     int N, Gx, Gy;
  811.     int dim[2], period[2], reorder;
  812.  
  813.     MPI_Comm comm;
  814.  
  815.     //N - global grid size
  816.     N = atoi(argv[1]);
  817.  
  818.     //Gx
  819.     Gx = dim[0] = atoi(argv[2]);
  820.  
  821.     //Gy
  822.     Gy = dim[1] = atoi(argv[3]);
  823.  
  824.     period[0]=0;
  825.     period[1]=0;
  826.  
  827.    
  828.  
  829.  
  830.     MPI_Init(&argc, &argv);
  831.  
  832.     time_start = MPI_Wtime();
  833.  
  834.     int world_size;
  835.  
  836.     int my_rank;
  837.  
  838.     MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  839.     MPI_Comm_size(MPI_COMM_WORLD, &world_size);    
  840.  
  841.  
  842.     if(my_rank == 0){
  843.         cout << "[INFO]: N = " << N << endl;
  844.         cout << "[INFO]: Gx = " << dim[0] << endl;
  845.         cout << "[INFO]: Gy = " << dim[1] << endl;
  846.     }
  847.  
  848.     if(argc <= 3){
  849.         if(my_rank == 0)
  850.             cout << "[ERROR]: Usage: mpieval <N> <Gx> <Gy>" << endl;
  851.  
  852.         MPI_Abort(MPI_COMM_WORLD, 1);
  853.     }
  854.  
  855.     if(Gx * Gy != world_size){
  856.         if(my_rank == 0)
  857.             cout << "[ERROR]: mpi world size is not equal to "<< Gx << "*" << Gy << endl;
  858.  
  859.         MPI_Abort(MPI_COMM_WORLD, 1);
  860.     }
  861.  
  862.     MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, 1, &comm);
  863.  
  864.     if(my_rank == 0)
  865.             cout << "[INFO]: Cart created"<<endl;
  866.  
  867.     MPI_Comm_rank(comm, &my_rank);
  868.  
  869.     int my_coords[2];
  870.  
  871.     MPI_Cart_coords(comm, my_rank, 2, my_coords);
  872.  
  873.  
  874.  
  875.     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);
  876.    
  877.  
  878.  
  879.     if ( my_rank == 0 ){
  880.         cout << "[INFO]: Running" << endl;
  881.     }
  882.  
  883.     int iteration_num = 0;
  884.  
  885.     bool should_i_stop = false;
  886.  
  887.     while ( should_i_stop != true ){
  888.         if ( (my_rank == 0) && (iteration_num % 10000 == 0) ){
  889.             cout << "[INFO]: Iteration " << iteration_num << endl;
  890.         }
  891.  
  892.         w_func.DoIteration(should_i_stop);
  893.        
  894.         iteration_num++;
  895.     }
  896.  
  897.  
  898.     loc_t1 = MPI_Wtime();
  899.  
  900.     MPI_Barrier(MPI_COMM_WORLD);
  901.  
  902.     loc_t2 = MPI_Wtime();
  903.     /*if(my_rank == 0)
  904.         cout << "[INFO]: Barrier jumping time: " << loc_t2 - loc_t1 << endl;
  905.     */
  906.     double comparing = w_func.CompareWithExact();
  907.  
  908.     if (my_rank == 0)
  909.         cout << "[INFO]: Diff with exact solution: " << comparing << endl;
  910.  
  911.  
  912.  
  913.     time_stop = MPI_Wtime();
  914.     if( my_rank == 0 )
  915.         cout << "Finished!" << endl
  916.             << "Total iterations: " << iteration_num << endl
  917.             << "Elapsed time: " << (time_stop - time_start) << endl
  918.             << "Time per iteration: " << (time_stop - time_start) / double(iteration_num) << endl;
  919.  
  920.     MPI_Finalize();
  921.  
  922.     return 0;
  923. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top