Advertisement
manac68974

Untitled

Dec 4th, 2019
125
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 26.10 KB | None | 0 0
  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 == M && current_node_global_offset_y == 0){
  43. return 6;
  44. }
  45.  
  46. //right top corner
  47. if (current_node_global_offset_x == M && 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 == M){
  61. return 3;
  62. }
  63.  
  64. //bottom side
  65. if (current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 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 <= M - 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 <= M - 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. if (neighbours["left"]){
  208. left_coeff *= -k(x_left + (glob_x - 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
  209. left_neighbour = Get(glob_x - 1, glob_y);
  210. }
  211.  
  212. if (neighbours["right"]){
  213. right_coeff *= -k(x_left + (glob_x + 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
  214. right_neighbour = Get(glob_x + 1, glob_y);
  215. }
  216.  
  217. if (neighbours["bottom"]){
  218. bottom_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y - 0.5) * h_y) / pow(h_y, 2);
  219. bottom_neighbour = Get(glob_x, glob_y - 1);
  220. }
  221.  
  222. if (neighbours["top"]){
  223. top_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y + 0.5) * h_y) / pow(h_y, 2);
  224. top_neighbour = Get(glob_x, glob_y + 1);
  225. }
  226.  
  227. this_coeff = q(x_left + glob_x * h_x, y_bottom + glob_y * h_y) - left_coeff - right_coeff - bottom_coeff - top_coeff;
  228. this_node = Get(glob_x, glob_y);
  229.  
  230. result = left_coeff * left_neighbour +
  231. right_coeff * right_neighbour +
  232. bottom_coeff * bottom_neighbour +
  233. top_coeff * top_neighbour +
  234. this_coeff * this_node;
  235.  
  236. return result;
  237. }
  238.  
  239. double GetNodeFromB(int current_node_global_offset_x, int current_node_global_offset_y) {
  240.  
  241. int glob_x = current_node_global_offset_x;
  242. int glob_y = current_node_global_offset_y;
  243.  
  244. double result = 0.0;
  245.  
  246. switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
  247. case 2:
  248. //left bottom corner
  249. result = F(x_left, y_bottom) + 2.0 / h_x * psi_L(y_bottom) + 2.0 / h_y * psi_B(x_left);
  250. break;
  251. case 4:
  252. //left top corner
  253. result = F(x_left, y_top) + 2.0 / h_x * psi_L(y_top) + 2.0 / h_y * psi_T(x_left);
  254. break;
  255. case 6:
  256. //right bottom corner
  257. result = F(x_right, y_bottom) + 2.0 / h_x * psi_R(y_bottom) + 2.0 / h_y * psi_B(x_right);
  258. break;
  259. case 8:
  260. //right top corner
  261. result = F(x_right, y_top) + 2.0 / h_x * psi_R(y_top) + 2.0 / h_y * psi_T(x_right);
  262. break;
  263. case 1:
  264. //left side
  265. result = F(x_left, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_L(y_bottom + glob_y * h_y);
  266. break;
  267. case 3:
  268. //right side
  269. result = F(x_right, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_R(y_bottom + glob_y * h_y);
  270. break;
  271. case 5:
  272. //bottom side
  273. result = F(x_left + glob_x * h_x, y_bottom) + 2.0 / h_y * psi_B(x_left + glob_x * h_x);
  274. break;
  275. case 7:
  276. //top side
  277. result = F(x_left + glob_x * h_x, y_top) + 2.0 / h_y * psi_T(x_left + glob_x * h_x);
  278. break;
  279. case 0:
  280. //internal
  281. result = F(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
  282. break;
  283. default:
  284. cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y <<endl;
  285.  
  286. }
  287.  
  288. return result;
  289.  
  290. }
  291.  
  292. double GetNodeFromExact(int current_node_global_offset_x, int current_node_global_offset_y) {
  293.  
  294. int glob_x = current_node_global_offset_x;
  295. int glob_y = current_node_global_offset_y;
  296.  
  297. return u(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
  298.  
  299. }
  300.  
  301. void ComputeMatrixR(){
  302. /*if(my_rank == 0)
  303. cout << "[INFO]: Computation of matrix r started"<<endl;
  304. */
  305.  
  306. vector<double> r_tmp_matrix (cur_block_size_x*cur_block_size_y, 0.0);
  307.  
  308. for(int i = 0; i < cur_block_size_x; ++i){
  309. for(int j = 0; j < cur_block_size_y; ++j){
  310.  
  311. int current_node_global_offset_x = GetGlobalX(i),
  312. current_node_global_offset_y = GetGlobalY(j);
  313.  
  314. int glob_x = current_node_global_offset_x,
  315. glob_y = current_node_global_offset_y;
  316.  
  317. r_tmp_matrix [ j + cur_block_size_y*i ] = ComputeMagicInnerProductA_iw(glob_x,glob_y) - GetNodeFromB(glob_x, glob_y);
  318.  
  319. }
  320. }
  321.  
  322. MPI_Barrier(MPI_COMM_WORLD);
  323. /*
  324. for(int i = 0; i < cur_block_size_x; ++i){
  325. for(int j = 0; j < cur_block_size_y; ++j){
  326.  
  327. old_internal_data[ j + cur_block_size_y*i ] = internal_data[ j + cur_block_size_y*i];
  328. internal_data[ j + cur_block_size_y*i ] = r_tmp_matrix[ j + cur_block_size_y*i];
  329.  
  330. }
  331. }*/
  332.  
  333. old_internal_data = internal_data;
  334. internal_data = r_tmp_matrix;
  335.  
  336. SyncMPI();
  337.  
  338. }
  339.  
  340. double ComputeTauAndStopCase(bool &should_i_stop){
  341.  
  342. double local_Ar_r_inner_product_sum = 0.0;
  343. double local_Ar_Ar_inner_product_sum = 0.0;
  344. double global_Ar_r_inner_product_sum = 0.0;
  345. double global_Ar_Ar_inner_product_sum = 0.0;
  346. double local_r_norm = 0.0;
  347. double global_r_norm = 0.0;
  348.  
  349.  
  350. for(int i = 0; i < cur_block_size_x; ++i){
  351. for(int j = 0; j < cur_block_size_y; ++j){
  352. double rho = 1.0;
  353.  
  354. int current_node_global_offset_x = GetGlobalX(i),
  355. current_node_global_offset_y = GetGlobalY(j);
  356.  
  357. int glob_x = current_node_global_offset_x,
  358. glob_y = current_node_global_offset_y;
  359.  
  360. double tmp_Ar_i_j = ComputeMagicInnerProductA_iw(glob_x, glob_y);
  361.  
  362. switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
  363. case 2:
  364. case 4:
  365. case 6:
  366. case 8:
  367. //angle
  368. rho = 0.25;
  369. break;
  370. case 1:
  371. case 3:
  372. case 5:
  373. case 7:
  374. //side
  375. rho = 0.5;
  376. break;
  377. case 0:
  378. //internal
  379. rho = 1.0;
  380. break;
  381. default:
  382. cout << "[ERROR]: Bad global coords compute tau. Global:" << glob_x << " " << glob_y << endl;
  383. }
  384.  
  385. double tmp_cur_node_value = Get(glob_x, glob_y);
  386.  
  387. local_Ar_r_inner_product_sum += rho * tmp_Ar_i_j * tmp_cur_node_value * h_x*h_y;
  388. local_Ar_Ar_inner_product_sum += rho * pow (tmp_Ar_i_j, 2) * h_x*h_y;
  389. local_r_norm += rho * pow(tmp_cur_node_value, 2) * h_x*h_y;
  390.  
  391. }
  392. }
  393.  
  394. //cout << "[DEBUG]: Local"<< local_Ar_r_inner_product_sum << endl;
  395.  
  396. MPI_Allreduce(&local_Ar_r_inner_product_sum, &global_Ar_r_inner_product_sum, 1, MPI_DOUBLE, MPI_SUM,
  397. comm);
  398.  
  399. //cout << "[DEBUG]: "<< global_Ar_r_inner_product_sum << endl;
  400.  
  401. MPI_Allreduce(&local_Ar_Ar_inner_product_sum, &global_Ar_Ar_inner_product_sum, 1, MPI_DOUBLE, MPI_SUM,
  402. MPI_COMM_WORLD);
  403.  
  404. //cout << "[DEBUG]: "<< global_Ar_Ar_inner_product_sum << endl;
  405.  
  406. double global_tau = global_Ar_r_inner_product_sum/ global_Ar_Ar_inner_product_sum;
  407.  
  408. MPI_Allreduce(&local_r_norm, &global_r_norm, 1, MPI_DOUBLE, MPI_SUM,
  409. MPI_COMM_WORLD);
  410.  
  411. double eps = 1e-06;
  412.  
  413. if (global_r_norm < 0){
  414. cout << "[ERROR]: bad global r norm" << endl;
  415. }
  416.  
  417. double tmp_norm = fabs(global_tau)*sqrt(global_r_norm);
  418.  
  419. if (my_rank == 0)
  420. cout << "[DEBUG]: solution norm "<< tmp_norm << endl;
  421.  
  422. if (tmp_norm <= eps){
  423. should_i_stop = true;
  424. }else{
  425. should_i_stop = false;
  426. }
  427.  
  428. return global_tau;
  429.  
  430. }
  431.  
  432.  
  433. void ComputeNewW(double tau){
  434.  
  435. for(int i = 0; i < cur_block_size_x; ++i){
  436. for(int j = 0; j < cur_block_size_y; ++j){
  437. 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 ];
  438. //old_internal_data[ j + cur_block_size_y*i ] = 0.0;
  439. }
  440. }
  441. }
  442.  
  443. int GetGlobalX(int i){
  444. return cur_block_global_offset_x + i;
  445. }
  446.  
  447. int GetGlobalY(int j){
  448. return cur_block_global_offset_y + j;
  449. }
  450.  
  451. public:
  452. 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){
  453.  
  454. M = inpM;
  455. N = inpN;
  456.  
  457. Gx = inpGx;
  458. Gy = inpGy;
  459.  
  460. x_left = inpx_left;
  461. x_right = inpx_right;
  462. y_bottom = inpy_bottom;
  463. y_top = inpy_top;
  464.  
  465. h_x = double((x_right - x_left)) / M;
  466. h_y = double((y_top - y_bottom)) / N;
  467.  
  468. my_rank = inprank;
  469. comm = inpcomm;
  470.  
  471. cur_block_global_coors_x = inp_cur_block_global_coors_x;
  472. cur_block_global_coors_y = inp_cur_block_global_coors_y;
  473.  
  474. global_block_size_x = (M + 1) / Gx;
  475. global_block_size_y = (N + 1) / Gy;
  476.  
  477. cur_block_size_x = global_block_size_x;
  478. cur_block_size_y = global_block_size_y;
  479.  
  480. cur_block_global_offset_x = global_block_size_x * cur_block_global_coors_x;
  481. cur_block_global_offset_y = global_block_size_y * cur_block_global_coors_y;
  482.  
  483. if (cur_block_global_offset_x + global_block_size_x > M + 1){
  484. cur_block_size_x = (M + 1) % Gx;
  485. }
  486.  
  487. if (cur_block_global_offset_y + global_block_size_y > N + 1){
  488. cur_block_size_y = (N + 1) % Gy;
  489. }
  490.  
  491. internal_data.resize(cur_block_size_x * cur_block_size_y);
  492. old_internal_data.resize(cur_block_size_x * cur_block_size_y);
  493.  
  494.  
  495. //OX
  496. external_data[0].resize(cur_block_size_y);
  497. external_data[1].resize(cur_block_size_y);
  498.  
  499. //OY
  500. external_data[2].resize(cur_block_size_x);
  501. external_data[3].resize(cur_block_size_x);
  502.  
  503. }
  504.  
  505. double Get(int i, int j) {
  506. return GetLocalIndex(i - cur_block_global_offset_x, j - cur_block_global_offset_y);
  507. }
  508.  
  509. void Set(int i, int j, double v) {
  510. return SetLocalIndex(i - cur_block_global_offset_x, j - cur_block_global_offset_y, v);
  511. }
  512.  
  513. void SyncMPI(){
  514.  
  515. /*if(my_rank == 0)
  516. cout << "[INFO]: Sync started"<< endl;
  517. */
  518.  
  519. MPI_Barrier(MPI_COMM_WORLD);
  520.  
  521. //left and right sides
  522. for(int j = 0; j < cur_block_size_y; ++j){
  523.  
  524. external_data[ 0 ][ j ] = GetLocalIndex(0,j);//internal_data[ j ];
  525. external_data[ 1 ][ j ] = GetLocalIndex(cur_block_size_x - 1,j);//internal_data[ j + cur_block_size_y * (cur_block_size_x - 1) ];
  526.  
  527. }
  528.  
  529. //bottom and top sides
  530. for(int i = 0; i < cur_block_size_x; ++i){
  531.  
  532. external_data[ 2 ][ i ] = GetLocalIndex(i,0);//internal_data[ cur_block_size_y*i ];
  533. external_data[ 3 ][ i ] = GetLocalIndex(i,cur_block_size_y - 1); //internal_data[ (cur_block_size_y - 1) + cur_block_size_y*i ];
  534.  
  535. }
  536.  
  537. int my_coords[2];
  538. int targets_ranks[4];
  539.  
  540. MPI_Cart_coords(comm, my_rank, 2, my_coords);
  541.  
  542. int neighbour_offsets[ 4 ][ 2 ] = {
  543. { -1, 0 },{ 1, 0 },
  544. { 0, -1 },{ 0, 1 }
  545. };
  546.  
  547. for(int i = 0; i < 4; i++){
  548.  
  549. int target_coords[2];
  550.  
  551. target_coords[ 0 ] = my_coords[ 0 ] + neighbour_offsets[ i ][ 0 ];
  552. target_coords[ 1 ] = my_coords[ 1 ] + neighbour_offsets[ i ][ 1 ];
  553.  
  554. if (target_coords[0] >= 0 && target_coords[0] < Gx && target_coords[1] >= 0 && target_coords[1] < Gy){
  555.  
  556. MPI_Cart_rank(comm, target_coords, &targets_ranks[ i ]);
  557.  
  558. }
  559. else{
  560. targets_ranks[i] = -1;
  561. }
  562. }
  563.  
  564. //Now we have rank for all targets
  565.  
  566. for(int axis = 0; axis < 2; axis++){
  567.  
  568. int parity_bit = (my_coords[ axis ]) % 2;
  569.  
  570. //if parity_bit == 0 then
  571. // zuerst mit links, dann mit rechts tauschen
  572. //elif parity_bit == 1:
  573. // zuerst mit rechts, dann mit links tauschen
  574.  
  575. for(int tmp = 0; tmp < 2; tmp++){
  576. parity_bit = 1 - parity_bit;
  577.  
  578. //target id in external_data und targets_ranks
  579. int target_idx = 2 * axis + (1 - parity_bit);
  580.  
  581. if (targets_ranks[target_idx] != -1){
  582.  
  583. // вычисляем теги отправки и приема
  584. // в них зашиты номер ноды, ось, направление
  585. int send_tag = 100000 + my_rank * 100 + axis * 10 + parity_bit;
  586. int recv_tag = 100000 + targets_ranks[ target_idx ] * 100 + axis * 10 + (1-parity_bit);
  587.  
  588. MPI_Status tmp_status;
  589. // если отправка не на себя, то отправляем
  590. if(my_rank != targets_ranks[ target_idx ]){
  591.  
  592. MPI_Sendrecv_replace(&external_data[ target_idx ][ 0 ], external_data[ target_idx ].size(),
  593. MPI_DOUBLE, targets_ranks[ target_idx ], send_tag, targets_ranks[ target_idx ], recv_tag,
  594. comm, &tmp_status);
  595.  
  596. }
  597. }
  598. }
  599. }
  600.  
  601. MPI_Barrier(MPI_COMM_WORLD);
  602.  
  603. }
  604.  
  605. void DoIteration(bool &should_i_stop){
  606.  
  607. ComputeMatrixR();
  608.  
  609. //Now R Matrix is in internal_data
  610.  
  611. double tau = ComputeTauAndStopCase(should_i_stop);
  612.  
  613. //We have in our block jetzt:
  614. //in internal_data: R Matrix
  615. //in old_internal_data: w aus letzte iteration
  616. //and we have tau
  617. //jetzt koennen wir naechste w finden
  618.  
  619. ComputeNewW(tau);
  620.  
  621. SyncMPI();
  622.  
  623. }
  624.  
  625. double GetLocalIndex(int i, int j){
  626. //internal data
  627. if ((j >= 0) && (j < cur_block_size_y) && (i >= 0) && (i < cur_block_size_x)){
  628. return internal_data[ j + cur_block_size_y*i ];
  629. }
  630.  
  631. //external data
  632. //OX
  633. if((j >= 0) && (j < cur_block_size_y)){
  634.  
  635. if (i == -1)
  636. return external_data[ 0 ][ j ];
  637.  
  638. if (i == cur_block_size_x)
  639. return external_data[ 1 ][ j ];
  640.  
  641. }
  642.  
  643. //OY
  644. if((i >= 0) && (i < cur_block_size_x)){
  645.  
  646. if (j == -1)
  647. return external_data[ 2 ][ i ];
  648. if (j == cur_block_size_y)
  649. return external_data[ 3 ][ i ];
  650.  
  651. }
  652.  
  653. cout << "[ERROR]: bad local index" << endl;
  654.  
  655. return nan("");
  656. }
  657.  
  658. void SetLocalIndex(int i, int j, double v){
  659. if ((j >= 0) && (j < cur_block_size_y) && (i >= 0) && (i < cur_block_size_x)){
  660. internal_data[ j + cur_block_size_y*i ] = v;
  661. }else{
  662. cout << "[ERROR]: trying to set data outside the local area" << endl;
  663. }
  664. }
  665.  
  666. double CompareWithExact() {
  667.  
  668. double local_diff_norm = 0.0;
  669. double global_diff_norm = 0.0;
  670. /*
  671. if (my_rank == 0)
  672. cout << "[INFO]: Starting computing compare with exact" << endl;
  673. */
  674. for (int i = 0; i < cur_block_size_x; ++i){
  675. for (int j = 0; j < cur_block_size_y; ++j){
  676.  
  677.  
  678. double rho = 1.0;
  679.  
  680. int current_node_global_offset_x = GetGlobalX(i),
  681. current_node_global_offset_y = GetGlobalY(j);
  682.  
  683. int glob_x = current_node_global_offset_x,
  684. glob_y = current_node_global_offset_y;
  685.  
  686. switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
  687. case 2:
  688. case 4:
  689. case 6:
  690. case 8:
  691. //angle
  692. rho = 0.25;
  693. break;
  694. case 1:
  695. case 3:
  696. case 5:
  697. case 7:
  698. //side
  699. rho = 0.5;
  700. break;
  701. case 0:
  702. //internal
  703. rho = 1.0;
  704. break;
  705. default:
  706. cout << "[ERROR]: Bad global coords compute exact. Local:" << i << " " << j << ". Global:" << glob_x << " " << glob_y <<endl;
  707. }
  708.  
  709. double tmp_elem = Get(glob_x, glob_y) - GetNodeFromExact(glob_x, glob_y);
  710.  
  711. //local_diff_norm = max( fabs(tmp_elem), local_diff_norm);
  712.  
  713. local_diff_norm += rho * pow(tmp_elem, 2) * h_x * h_y;
  714.  
  715. }
  716. }
  717.  
  718. //cout << "[INFO]: local max diff in " << cur_block_global_offset_x << " " << cur_block_global_offset_y << " " << local_diff_norm << endl;
  719.  
  720. MPI_Allreduce(&local_diff_norm, &global_diff_norm, 1, MPI_DOUBLE, MPI_SUM,
  721. MPI_COMM_WORLD);
  722.  
  723. global_diff_norm = sqrt(global_diff_norm);
  724.  
  725. return global_diff_norm;
  726. }
  727. };
  728.  
  729.  
  730.  
  731.  
  732. int main(int argc, char* argv[]){
  733.  
  734. const double x_left = -2, x_right = 3;
  735. const double y_bottom = -1, y_top = 4;
  736.  
  737. double time_start, time_stop;
  738. double loc_t1, loc_t2;
  739. int N, Gx, Gy;
  740. int dim[2], period[2], reorder;
  741.  
  742. MPI_Comm comm;
  743.  
  744. //N - global grid size
  745. N = atoi(argv[1]);
  746.  
  747. //Gx
  748. Gx = dim[0] = atoi(argv[2]);
  749.  
  750. //Gy
  751. Gy = dim[1] = atoi(argv[3]);
  752.  
  753. period[0]=0;
  754. period[1]=0;
  755.  
  756.  
  757.  
  758.  
  759. MPI_Init(&argc, &argv);
  760.  
  761. time_start = MPI_Wtime();
  762.  
  763. int world_size;
  764.  
  765. int my_rank;
  766.  
  767. MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  768. MPI_Comm_size(MPI_COMM_WORLD, &world_size);
  769.  
  770.  
  771. if(my_rank == 0){
  772. cout << "[INFO]: N = " << N << endl;
  773. cout << "[INFO]: Gx = " << dim[0] << endl;
  774. cout << "[INFO]: Gy = " << dim[1] << endl;
  775. }
  776.  
  777. if(argc <= 3){
  778. if(my_rank == 0)
  779. cout << "[ERROR]: Usage: mpieval <N> <Gx> <Gy>" << endl;
  780.  
  781. MPI_Abort(MPI_COMM_WORLD, 1);
  782. }
  783.  
  784. if(Gx * Gy != world_size){
  785. if(my_rank == 0)
  786. cout << "[ERROR]: mpi world size is not equal to "<< Gx << "*" << Gy << endl;
  787.  
  788. MPI_Abort(MPI_COMM_WORLD, 1);
  789. }
  790.  
  791. MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, 1, &comm);
  792.  
  793. if(my_rank == 0)
  794. cout << "[INFO]: Cart created"<<endl;
  795.  
  796. MPI_Comm_rank(comm, &my_rank);
  797.  
  798. int my_coords[2];
  799.  
  800. MPI_Cart_coords(comm, my_rank, 2, my_coords);
  801.  
  802.  
  803.  
  804. 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);
  805.  
  806.  
  807.  
  808. if ( my_rank == 0 ){
  809. cout << "[INFO]: Running" << endl;
  810. }
  811.  
  812. int iteration_num = 0;
  813.  
  814. bool should_i_stop = false;
  815.  
  816. while ( should_i_stop != true ){
  817. if ( (my_rank == 0) && (iteration_num % 10000 == 0) ){
  818. cout << "[INFO]: Iteration " << iteration_num << endl;
  819. }
  820.  
  821. w_func.DoIteration(should_i_stop);
  822.  
  823. iteration_num++;
  824. }
  825.  
  826.  
  827. loc_t1 = MPI_Wtime();
  828.  
  829. MPI_Barrier(MPI_COMM_WORLD);
  830.  
  831. loc_t2 = MPI_Wtime();
  832. /*if(my_rank == 0)
  833. cout << "[INFO]: Barrier jumping time: " << loc_t2 - loc_t1 << endl;
  834. */
  835. double comparing = w_func.CompareWithExact();
  836.  
  837. if (my_rank == 0)
  838. cout << "[INFO]: Diff with exact solution: " << comparing << endl;
  839.  
  840.  
  841.  
  842. time_stop = MPI_Wtime();
  843. if( my_rank == 0 )
  844. cout << "Finished!" << endl
  845. << "Total iterations: " << iteration_num << endl
  846. << "Elapsed time: " << (time_stop - time_start) << endl
  847. << "Time per iteration: " << (time_stop - time_start) / double(iteration_num) << endl;
  848.  
  849. MPI_Finalize();
  850.  
  851. return 0;
  852. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement