API tools faq
paste
Advertisement
manac68974

Untitled

manac68974
Dec 3rd, 2019
143
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 27.90 KB | None | 0 0
raw download clone embed print report
  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. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!