fmm_tree.txx 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686
  1. /**
  2. * \file fmm_tree.txx
  3. * \author Dhairya Malhotra, dhairya.malhotra@gmail.com
  4. * \date 12-11-2010
  5. * \brief This file contains the implementation of the class FMM_Tree.
  6. */
  7. #include <omp.h>
  8. #include <sstream>
  9. #include <iomanip>
  10. #include <cassert>
  11. #include <mpi_node.hpp>
  12. #include <fmm_node.hpp>
  13. #include <mem_utils.hpp>
  14. #include <mortonid.hpp>
  15. #include <profile.hpp>
  16. #include <vector.hpp>
  17. namespace pvfmm{
  18. template <class FMM_Mat_t>
  19. void FMM_Tree<FMM_Mat_t>::Initialize(typename FMM_Node_t::NodeData* init_data) {
  20. Profile::Tic("InitTree",this->Comm(),true);{
  21. //Build octree from points.
  22. MPI_Tree<FMM_Node_t>::Initialize(init_data);
  23. Profile::Tic("InitFMMData",this->Comm(),true,5);
  24. { //Initialize FMM data.
  25. std::vector<FMM_Node_t*>& nodes=this->GetNodeList();
  26. #pragma omp parallel for
  27. for(size_t i=0;i<nodes.size();i++){
  28. if(nodes[i]->FMMData()==NULL) nodes[i]->FMMData()=new typename FMM_Mat_t::FMMData;
  29. }
  30. }
  31. Profile::Toc();
  32. }Profile::Toc();
  33. }
  34. template <class FMM_Mat_t>
  35. void FMM_Tree<FMM_Mat_t>::InitFMM_Tree(bool refine, BoundaryType bndry_) {
  36. Profile::Tic("InitFMM_Tree",this->Comm(),true);{
  37. interac_list.Initialize(this->Dim());
  38. bndry=bndry_;
  39. if(refine){
  40. //RefineTree
  41. Profile::Tic("RefineTree",this->Comm(),true,2);
  42. this->RefineTree();
  43. Profile::Toc();
  44. }
  45. //2:1 Balancing
  46. Profile::Tic("2:1Balance",this->Comm(),true,2);
  47. this->Balance21(bndry);
  48. Profile::Toc();
  49. //Redistribute nodes.
  50. Profile::Tic("Redistribute",this->Comm(),true,3);
  51. this->RedistNodes();
  52. Profile::Toc();
  53. }Profile::Toc();
  54. }
  55. template <class FMM_Mat_t>
  56. void FMM_Tree<FMM_Mat_t>::SetupFMM(FMM_Mat_t* fmm_mat_) {
  57. Profile::Tic("SetupFMM",this->Comm(),true);{
  58. bool device=true;
  59. #ifdef __INTEL_OFFLOAD
  60. Profile::Tic("InitLocks",this->Comm(),false,3);
  61. MIC_Lock::init();
  62. Profile::Toc();
  63. #endif
  64. int omp_p=omp_get_max_threads();
  65. fmm_mat=fmm_mat_;
  66. //Construct LET
  67. Profile::Tic("ConstructLET",this->Comm(),false,2);
  68. this->ConstructLET(bndry);
  69. Profile::Toc();
  70. //Set Colleagues (Needed to build U, V, W and X lists.)
  71. Profile::Tic("SetColleagues",this->Comm(),false,3);
  72. this->SetColleagues(bndry);
  73. Profile::Toc();
  74. Profile::Tic("BuildLists",this->Comm(),false,3);
  75. BuildInteracLists();
  76. Profile::Toc();
  77. Profile::Tic("CollectNodeData",this->Comm(),false,3);
  78. //Build node list.
  79. FMM_Node_t* n=dynamic_cast<FMM_Node_t*>(this->PostorderFirst());
  80. std::vector<FMM_Node_t*> all_nodes;
  81. while(n!=NULL){
  82. all_nodes.push_back(n);
  83. n=static_cast<FMM_Node_t*>(this->PostorderNxt(n));
  84. }
  85. //Collect node data into continuous array.
  86. std::vector<Vector<FMM_Node_t*> > node_lists; // TODO: Remove this parameter, not really needed
  87. fmm_mat->CollectNodeData(all_nodes, node_data_buff, node_lists);
  88. Profile::Toc();
  89. setup_data.clear();
  90. precomp_lst.clear();
  91. setup_data.resize(8*MAX_DEPTH);
  92. precomp_lst.resize(8);
  93. Profile::Tic("UListSetup",this->Comm(),false,3);
  94. for(size_t i=0;i<MAX_DEPTH;i++){
  95. setup_data[i+MAX_DEPTH*0].precomp_data=&precomp_lst[0];
  96. fmm_mat->U_ListSetup(setup_data[i+MAX_DEPTH*0],node_data_buff,node_lists,fmm_mat->Homogen()?(i==0?-1:MAX_DEPTH+1):i, device);
  97. }
  98. Profile::Toc();
  99. Profile::Tic("WListSetup",this->Comm(),false,3);
  100. for(size_t i=0;i<MAX_DEPTH;i++){
  101. setup_data[i+MAX_DEPTH*1].precomp_data=&precomp_lst[1];
  102. fmm_mat->W_ListSetup(setup_data[i+MAX_DEPTH*1],node_data_buff,node_lists,fmm_mat->Homogen()?(i==0?-1:MAX_DEPTH+1):i, device);
  103. }
  104. Profile::Toc();
  105. Profile::Tic("XListSetup",this->Comm(),false,3);
  106. for(size_t i=0;i<MAX_DEPTH;i++){
  107. setup_data[i+MAX_DEPTH*2].precomp_data=&precomp_lst[2];
  108. fmm_mat->X_ListSetup(setup_data[i+MAX_DEPTH*2],node_data_buff,node_lists,fmm_mat->Homogen()?(i==0?-1:MAX_DEPTH+1):i, device);
  109. }
  110. Profile::Toc();
  111. Profile::Tic("VListSetup",this->Comm(),false,3);
  112. for(size_t i=0;i<MAX_DEPTH;i++){
  113. setup_data[i+MAX_DEPTH*3].precomp_data=&precomp_lst[3];
  114. fmm_mat->V_ListSetup(setup_data[i+MAX_DEPTH*3],node_data_buff,node_lists,fmm_mat->Homogen()?(i==0?-1:MAX_DEPTH+1):i, /*device*/ false);
  115. }
  116. Profile::Toc();
  117. Profile::Tic("D2DSetup",this->Comm(),false,3);
  118. for(size_t i=0;i<MAX_DEPTH;i++){
  119. setup_data[i+MAX_DEPTH*4].precomp_data=&precomp_lst[4];
  120. fmm_mat->Down2DownSetup(setup_data[i+MAX_DEPTH*4],node_data_buff,node_lists,i, /*device*/ false);
  121. }
  122. Profile::Toc();
  123. Profile::Tic("D2TSetup",this->Comm(),false,3);
  124. for(size_t i=0;i<MAX_DEPTH;i++){
  125. setup_data[i+MAX_DEPTH*5].precomp_data=&precomp_lst[5];
  126. fmm_mat->Down2TargetSetup(setup_data[i+MAX_DEPTH*5],node_data_buff,node_lists,fmm_mat->Homogen()?(i==0?-1:MAX_DEPTH+1):i, /*device*/ false);
  127. }
  128. Profile::Toc();
  129. Profile::Tic("S2USetup",this->Comm(),false,3);
  130. for(size_t i=0;i<MAX_DEPTH;i++){
  131. setup_data[i+MAX_DEPTH*6].precomp_data=&precomp_lst[6];
  132. fmm_mat->Source2UpSetup(setup_data[i+MAX_DEPTH*6],node_data_buff,node_lists,fmm_mat->Homogen()?(i==0?-1:MAX_DEPTH+1):i, /*device*/ false);
  133. }
  134. Profile::Toc();
  135. Profile::Tic("U2USetup",this->Comm(),false,3);
  136. for(size_t i=0;i<MAX_DEPTH;i++){
  137. setup_data[i+MAX_DEPTH*7].precomp_data=&precomp_lst[7];
  138. fmm_mat->Up2UpSetup(setup_data[i+MAX_DEPTH*7],node_data_buff,node_lists,i, /*device*/ false);
  139. }
  140. Profile::Toc();
  141. #ifdef __INTEL_OFFLOAD
  142. int wait_lock_idx=-1;
  143. wait_lock_idx=MIC_Lock::curr_lock();
  144. #pragma offload target(mic:0)
  145. {MIC_Lock::wait_lock(wait_lock_idx);}
  146. #endif
  147. }Profile::Toc();
  148. }
  149. template <class FMM_Mat_t>
  150. void FMM_Tree<FMM_Mat_t>::ClearFMMData() {
  151. Profile::Tic("ClearFMMData",this->Comm(),true);{
  152. bool device=true;
  153. if(setup_data[0+MAX_DEPTH*1]. input_data!=NULL) setup_data[0+MAX_DEPTH*1]. input_data->SetZero();
  154. if(setup_data[0+MAX_DEPTH*2].output_data!=NULL) setup_data[0+MAX_DEPTH*2].output_data->SetZero();
  155. if(setup_data[0+MAX_DEPTH*0].output_data!=NULL) setup_data[0+MAX_DEPTH*0].output_data->SetZero();
  156. if(device){ // Host2Device
  157. if(setup_data[0+MAX_DEPTH*1]. input_data!=NULL) setup_data[0+MAX_DEPTH*1]. input_data->AllocDevice(true);
  158. if(setup_data[0+MAX_DEPTH*2].output_data!=NULL) setup_data[0+MAX_DEPTH*2].output_data->AllocDevice(true);
  159. if(setup_data[0+MAX_DEPTH*0].output_data!=NULL) setup_data[0+MAX_DEPTH*0].output_data->AllocDevice(true);
  160. #ifdef __INTEL_OFFLOAD
  161. if(!fmm_mat->Homogen()){ // Wait
  162. int wait_lock_idx=-1;
  163. wait_lock_idx=MIC_Lock::curr_lock();
  164. #pragma offload target(mic:0)
  165. {MIC_Lock::wait_lock(wait_lock_idx);}
  166. }
  167. MIC_Lock::init();
  168. #endif
  169. }
  170. }Profile::Toc();
  171. }
  172. template <class FMM_Mat_t>
  173. void FMM_Tree<FMM_Mat_t>::RunFMM() {
  174. Profile::Tic("RunFMM",this->Comm(),true);{
  175. //Upward Pass
  176. Profile::Tic("UpwardPass",this->Comm(),false,2);
  177. UpwardPass();
  178. Profile::Toc();
  179. //Multipole Reduce Broadcast.
  180. Profile::Tic("ReduceBcast",this->Comm(),true,2);
  181. MultipoleReduceBcast();
  182. Profile::Toc();
  183. //Local 2:1 Balancing.
  184. //This can cause load imbalance, always use global 2:1 balance instead.
  185. //Profile::Tic("2:1Balance(local)",this->Comm(),false,3);
  186. //this->Balance21_local(bndry);
  187. //UpwardPass(true);
  188. //Profile::Toc();
  189. //Downward Pass
  190. Profile::Tic("DownwardPass",this->Comm(),true,2);
  191. DownwardPass();
  192. Profile::Toc();
  193. }Profile::Toc();
  194. }
  195. template <class FMM_Mat_t>
  196. void FMM_Tree<FMM_Mat_t>::UpwardPass() {
  197. bool device=true;
  198. //Upward Pass (initialize all leaf nodes)
  199. Profile::Tic("S2U",this->Comm(),false,5);
  200. for(int i=0; i<(fmm_mat->Homogen()?1:MAX_DEPTH); i++){ // Source2Up
  201. if(!fmm_mat->Homogen()) fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*6],/*device*/ false);
  202. fmm_mat->Source2Up(setup_data[i+MAX_DEPTH*6]);
  203. }
  204. Profile::Toc();
  205. //Upward Pass (level by level)
  206. Profile::Tic("U2U",this->Comm(),false,5);
  207. for(int i=MAX_DEPTH-1; i>=0; i--){ // Up2Up
  208. if(!fmm_mat->Homogen()) fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*7],/*device*/ false);
  209. fmm_mat->Up2Up(setup_data[i+MAX_DEPTH*7]);
  210. }
  211. Profile::Toc();
  212. }
  213. template <class FMM_Mat_t>
  214. void FMM_Tree<FMM_Mat_t>::BuildInteracLists() {
  215. std::vector<FMM_Node_t*>& n_list=this->GetNodeList();
  216. // Build interaction lists.
  217. int omp_p=omp_get_max_threads();
  218. {
  219. size_t k=n_list.size();
  220. #pragma omp parallel for
  221. for(int j=0;j<omp_p;j++){
  222. size_t a=(k*j)/omp_p;
  223. size_t b=(k*(j+1))/omp_p;
  224. for(size_t i=a;i<b;i++){
  225. FMM_Node_t* n=n_list[i];
  226. n->interac_list.resize(Type_Count);
  227. n->interac_list[S2U_Type]=interac_list.BuildList(n,S2U_Type);
  228. n->interac_list[U2U_Type]=interac_list.BuildList(n,U2U_Type);
  229. n->interac_list[D2D_Type]=interac_list.BuildList(n,D2D_Type);
  230. n->interac_list[D2T_Type]=interac_list.BuildList(n,D2T_Type);
  231. n->interac_list[U0_Type]=interac_list.BuildList(n,U0_Type);
  232. n->interac_list[U1_Type]=interac_list.BuildList(n,U1_Type);
  233. n->interac_list[U2_Type]=interac_list.BuildList(n,U2_Type);
  234. n->interac_list[V_Type]=interac_list.BuildList(n,V_Type);
  235. n->interac_list[V1_Type]=interac_list.BuildList(n,V1_Type);
  236. n->interac_list[W_Type]=interac_list.BuildList(n,W_Type);
  237. n->interac_list[X_Type]=interac_list.BuildList(n,X_Type);
  238. }
  239. }
  240. }
  241. }
  242. template <class FMM_Mat_t>
  243. void FMM_Tree<FMM_Mat_t>::MultipoleReduceBcast() {
  244. int num_p,rank;
  245. MPI_Comm_size(*this->Comm(),&num_p);
  246. MPI_Comm_rank(*this->Comm(),&rank );
  247. if(num_p==1) return;
  248. Profile::Tic("Reduce",this->Comm(),true,3);
  249. std::vector<MortonId> mins=this->GetMins();
  250. size_t bit_mask=1;
  251. size_t max_child=(1UL<<this->Dim());
  252. //Initialize initial send nodes.
  253. std::vector<FMM_Node_t*> send_nodes[2];
  254. //Initialize send_node[0]
  255. FMM_Node_t* tmp_node=static_cast<FMM_Node_t*>(this->RootNode());
  256. assert(!tmp_node->IsGhost());
  257. while(!tmp_node->IsLeaf()){
  258. FMM_Node_t* tmp_node_=NULL;
  259. for(size_t i=0;i<max_child;i++){
  260. tmp_node_=static_cast<FMM_Node_t*>(tmp_node->Child(i));
  261. if(tmp_node_!=NULL) if(!tmp_node_->IsGhost()) break;
  262. }
  263. tmp_node=tmp_node_; assert(tmp_node!=NULL);
  264. }
  265. int n[2];
  266. n[0]=tmp_node->Depth()+1;
  267. send_nodes[0].resize(n[0]);
  268. send_nodes[0][n[0]-1]=tmp_node;
  269. for(int i=n[0]-1;i>0;i--)
  270. send_nodes[0][i-1]=static_cast<FMM_Node_t*>(send_nodes[0][i]->Parent());
  271. //Initialize send_node[1]
  272. tmp_node=static_cast<FMM_Node_t*>(this->RootNode());
  273. while(!tmp_node->IsLeaf()){
  274. FMM_Node_t* tmp_node_=NULL;
  275. for(int i=max_child-1;i>=0;i--){
  276. tmp_node_=static_cast<FMM_Node_t*>(tmp_node->Child(i));
  277. if(tmp_node_!=NULL) if(!tmp_node_->IsGhost()) break;
  278. }
  279. tmp_node=tmp_node_; assert(tmp_node!=NULL);
  280. }
  281. n[1]=tmp_node->Depth()+1;
  282. send_nodes[1].resize(n[1]);
  283. send_nodes[1][n[1]-1]=tmp_node;
  284. for(int i=n[1]-1;i>0;i--)
  285. send_nodes[1][i-1]=static_cast<FMM_Node_t*>(send_nodes[1][i]->Parent());
  286. //Hypercube reduction.
  287. while(bit_mask<(size_t)num_p){
  288. size_t partner=rank^bit_mask; //Partner process id
  289. int merge_indx=((bit_mask & rank)==0?1:0);
  290. bit_mask=bit_mask<<1;
  291. //if(rank >= num_p - (num_p % bit_mask)) break;
  292. //Initialize send data.
  293. size_t s_node_cnt[2]={send_nodes[0].size(),send_nodes[1].size()};
  294. int send_size=2*sizeof(size_t)+(s_node_cnt[0]+s_node_cnt[1])*sizeof(MortonId);
  295. std::vector<PackedData> send_data(s_node_cnt[0]+s_node_cnt[1]);
  296. size_t s_iter=0;
  297. for(int i=0;i<2;i++)
  298. for(size_t j=0;j<s_node_cnt[i];j++){
  299. assert(send_nodes[i][j]!=NULL);
  300. send_data[s_iter]=send_nodes[i][j]->PackMultipole();
  301. send_size+=send_data[s_iter].length+sizeof(size_t);
  302. s_iter++;
  303. }
  304. char* send_buff=new char[send_size];
  305. char* buff_iter=send_buff;
  306. ((size_t*)buff_iter)[0]=s_node_cnt[0];
  307. ((size_t*)buff_iter)[1]=s_node_cnt[1];
  308. buff_iter+=2*sizeof(size_t);
  309. s_iter=0;
  310. for(int i=0;i<2;i++)
  311. for(size_t j=0;j<s_node_cnt[i];j++){
  312. ((MortonId*)buff_iter)[0]=send_nodes[i][j]->GetMortonId();
  313. buff_iter+=sizeof(MortonId);
  314. ((size_t*)buff_iter)[0]=send_data[s_iter].length;
  315. buff_iter+=sizeof(size_t);
  316. mem::memcopy((void*)buff_iter,send_data[s_iter].data,send_data[s_iter].length);
  317. buff_iter+=send_data[s_iter].length;
  318. s_iter++;
  319. }
  320. //Exchange send and recv sizes
  321. int recv_size=0;
  322. MPI_Status status;
  323. char* recv_buff=NULL;
  324. if(partner<(size_t)num_p){
  325. MPI_Sendrecv(&send_size, 1, MPI_INT, partner, 0, &recv_size, 1, MPI_INT, partner, 0, *this->Comm(), &status);
  326. recv_buff=new char[recv_size];
  327. MPI_Sendrecv(send_buff, send_size, MPI_BYTE, partner, 0, recv_buff, recv_size, MPI_BYTE, partner, 0, *this->Comm(), &status);
  328. }
  329. //Need an extra broadcast for incomplete hypercubes.
  330. size_t p0_start=num_p - (num_p % (bit_mask ));
  331. size_t p0_end =num_p - (num_p % (bit_mask>>1));
  332. if(((size_t)rank >= p0_start) && ((size_t)num_p>p0_end) && ((size_t)rank < p0_end) ){
  333. size_t bit_mask0=1;
  334. size_t num_p0=p0_end-p0_start;
  335. while( bit_mask0 < num_p0 ){
  336. if( (bit_mask0<<1) > (num_p - p0_end) ){
  337. size_t partner0=rank^bit_mask0;
  338. if( rank-p0_start < bit_mask0 ){
  339. //Send
  340. MPI_Send(&recv_size, 1, MPI_INT , partner0, 0, *this->Comm());
  341. MPI_Send( recv_buff, recv_size, MPI_BYTE, partner0, 0, *this->Comm());
  342. }else if( rank-p0_start < (bit_mask0<<1) ){
  343. //Receive
  344. if(recv_size>0) delete[] recv_buff;
  345. MPI_Recv(&recv_size, 1, MPI_INT , partner0, 0, *this->Comm(), &status);
  346. recv_buff=new char[recv_size];
  347. MPI_Recv( recv_buff, recv_size, MPI_BYTE, partner0, 0, *this->Comm(), &status);
  348. }
  349. }
  350. bit_mask0=bit_mask0<<1;
  351. }
  352. }
  353. //Construct nodes from received data.
  354. if(recv_size>0){
  355. buff_iter=recv_buff;
  356. size_t r_node_cnt[2]={((size_t*)buff_iter)[0],((size_t*)buff_iter)[1]};
  357. buff_iter+=2*sizeof(size_t);
  358. std::vector<MortonId> r_mid[2];
  359. r_mid[0].resize(r_node_cnt[0]);
  360. r_mid[1].resize(r_node_cnt[1]);
  361. std::vector<FMM_Node_t*> recv_nodes[2];
  362. recv_nodes[0].resize(r_node_cnt[0]);
  363. recv_nodes[1].resize(r_node_cnt[1]);
  364. std::vector<PackedData> recv_data[2];
  365. recv_data[0].resize(r_node_cnt[0]);
  366. recv_data[1].resize(r_node_cnt[1]);
  367. for(int i=0;i<2;i++)
  368. for(size_t j=0;j<r_node_cnt[i];j++){
  369. r_mid[i][j]=((MortonId*)buff_iter)[0];
  370. buff_iter+=sizeof(MortonId);
  371. recv_data[i][j].length=((size_t*)buff_iter)[0];
  372. buff_iter+=sizeof(size_t);
  373. recv_data[i][j].data=(void*)buff_iter;
  374. buff_iter+=recv_data[i][j].length;
  375. }
  376. // Add multipole expansion to existing nodes.
  377. for(size_t i=0;i<r_node_cnt[1-merge_indx];i++){
  378. if(i<send_nodes[merge_indx].size()){
  379. if(r_mid[1-merge_indx][i]==send_nodes[merge_indx][i]->GetMortonId()){
  380. send_nodes[merge_indx][i]->AddMultipole(recv_data[1-merge_indx][i]);
  381. }else break;
  382. }else break;
  383. }
  384. bool new_branch=false;
  385. for(size_t i=0;i<r_node_cnt[merge_indx];i++){
  386. if(i<send_nodes[merge_indx].size() && !new_branch){
  387. if(r_mid[merge_indx][i]==send_nodes[merge_indx][i]->GetMortonId()){
  388. recv_nodes[merge_indx][i]=send_nodes[merge_indx][i];
  389. }else{
  390. new_branch=true;
  391. size_t n_=(i<(size_t)n[merge_indx]?n[merge_indx]:i);
  392. for(size_t j=n_;j<send_nodes[merge_indx].size();j++)
  393. delete send_nodes[merge_indx][j];
  394. if(i<(size_t)n[merge_indx]) n[merge_indx]=i;
  395. }
  396. }
  397. if(i>=send_nodes[merge_indx].size() || new_branch){
  398. recv_nodes[merge_indx][i]=static_cast<FMM_Node_t*>(this->NewNode());
  399. recv_nodes[merge_indx][i]->SetCoord(r_mid[merge_indx][i]);
  400. recv_nodes[merge_indx][i]->InitMultipole(recv_data[merge_indx][i]);
  401. }
  402. }
  403. send_nodes[merge_indx]=recv_nodes[merge_indx];
  404. }
  405. delete[] send_buff;
  406. delete[] recv_buff;
  407. }
  408. for(int i=0;i<2;i++)
  409. for(size_t j=n[i];j<send_nodes[i].size();j++)
  410. delete send_nodes[i][j];
  411. Profile::Toc();
  412. //Now Broadcast nodes to build LET.
  413. Profile::Tic("Broadcast",this->Comm(),true,4);
  414. this->ConstructLET(bndry);
  415. Profile::Toc();
  416. }
  417. template <class FMM_Mat_t>
  418. void FMM_Tree<FMM_Mat_t>::DownwardPass() {
  419. bool device=true;
  420. Profile::Tic("Setup",this->Comm(),true,3);
  421. std::vector<FMM_Node_t*> leaf_nodes;
  422. int max_depth=0;
  423. { // Build leaf node list
  424. int max_depth_loc=0;
  425. std::vector<FMM_Node_t*>& nodes=this->GetNodeList();
  426. for(size_t i=0;i<nodes.size();i++){
  427. FMM_Node_t* n=nodes[i];
  428. if(!n->IsGhost() && n->IsLeaf()) leaf_nodes.push_back(n);
  429. if(n->Depth()>max_depth_loc) max_depth_loc=n->Depth();
  430. }
  431. MPI_Allreduce(&max_depth_loc, &max_depth, 1, MPI_INT, MPI_MAX, *this->Comm());
  432. }
  433. Profile::Toc();
  434. #ifdef __INTEL_OFFLOAD
  435. if(device){ // Host2Device:Src
  436. Profile::Tic("Host2Device:Src",this->Comm(),false,5);
  437. if(setup_data[0+MAX_DEPTH*2]. coord_data!=NULL) setup_data[0+MAX_DEPTH*2]. coord_data->AllocDevice(true);
  438. if(setup_data[0+MAX_DEPTH*2]. input_data!=NULL) setup_data[0+MAX_DEPTH*2]. input_data->AllocDevice(true);
  439. Profile::Toc();
  440. }
  441. #endif
  442. if(bndry==Periodic){ //Add contribution from periodic infinite tiling.
  443. Profile::Tic("BoundaryCondition",this->Comm(),false,5);
  444. fmm_mat->PeriodicBC(dynamic_cast<FMM_Node_t*>(this->RootNode()));
  445. Profile::Toc();
  446. }
  447. for(size_t i=0; i<=(fmm_mat->Homogen()?0:max_depth); i++){ // U,V,W,X-lists
  448. if(!fmm_mat->Homogen()){ // Precomp
  449. std::stringstream level_str;
  450. level_str<<"Level-"<<std::setfill('0')<<std::setw(2)<<i<<"\0";
  451. Profile::Tic(level_str.str().c_str(),this->Comm(),false,5);
  452. Profile::Tic("Precomp",this->Comm(),false,5);
  453. {// Precomp U
  454. Profile::Tic("Precomp-U",this->Comm(),false,10);
  455. fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*0],device);
  456. Profile::Toc();
  457. }
  458. {// Precomp W
  459. Profile::Tic("Precomp-W",this->Comm(),false,10);
  460. fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*1],device);
  461. Profile::Toc();
  462. }
  463. {// Precomp X
  464. Profile::Tic("Precomp-X",this->Comm(),false,10);
  465. fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*2],device);
  466. Profile::Toc();
  467. }
  468. {// Precomp V
  469. Profile::Tic("Precomp-V",this->Comm(),false,10);
  470. fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*3], /*device*/ false);
  471. Profile::Toc();
  472. }
  473. Profile::Toc();
  474. }
  475. {// X-List
  476. Profile::Tic("X-List",this->Comm(),false,5);
  477. fmm_mat->X_List(setup_data[i+MAX_DEPTH*2], device);
  478. Profile::Toc();
  479. }
  480. #ifdef __INTEL_OFFLOAD
  481. if(i==0 && device){ // Host2Device:Mult
  482. Profile::Tic("Host2Device:Mult",this->Comm(),false,5);
  483. if(setup_data[0+MAX_DEPTH*1]. input_data!=NULL) setup_data[0+MAX_DEPTH*1]. input_data->AllocDevice(true);
  484. Profile::Toc();
  485. }
  486. if(device) if(i==(fmm_mat->Homogen()?0:max_depth)){ // Device2Host: LocalExp
  487. Profile::Tic("Device2Host:LocExp",this->Comm(),false,5);
  488. if(setup_data[0+MAX_DEPTH*2].output_data!=NULL){
  489. Matrix<Real_t>& output_data=*setup_data[0+MAX_DEPTH*2].output_data;
  490. assert(fmm_mat->dev_buffer.Dim()>=output_data.Dim(0)*output_data.Dim(1));
  491. output_data.Device2Host((Real_t*)&fmm_mat->dev_buffer[0]);
  492. }
  493. Profile::Toc();
  494. }
  495. #endif
  496. {// W-List
  497. Profile::Tic("W-List",this->Comm(),false,5);
  498. fmm_mat->W_List(setup_data[i+MAX_DEPTH*1], device);
  499. Profile::Toc();
  500. }
  501. {// U-List
  502. Profile::Tic("U-List",this->Comm(),false,5);
  503. fmm_mat->U_List(setup_data[i+MAX_DEPTH*0], device);
  504. Profile::Toc();
  505. }
  506. {// V-List
  507. Profile::Tic("V-List",this->Comm(),false,5);
  508. fmm_mat->V_List(setup_data[i+MAX_DEPTH*3], /*device*/ false);
  509. Profile::Toc();
  510. }
  511. if(!fmm_mat->Homogen()){ // Wait
  512. #ifdef __INTEL_OFFLOAD
  513. int wait_lock_idx=-1;
  514. if(device) wait_lock_idx=MIC_Lock::curr_lock();
  515. #pragma offload if(device) target(mic:0)
  516. {if(device) MIC_Lock::wait_lock(wait_lock_idx);}
  517. #endif
  518. Profile::Toc();
  519. }
  520. }
  521. #ifdef __INTEL_OFFLOAD
  522. Profile::Tic("D2H_Wait:LocExp",this->Comm(),false,5);
  523. if(device) if(setup_data[0+MAX_DEPTH*2].output_data!=NULL){
  524. Real_t* dev_ptr=(Real_t*)&fmm_mat->dev_buffer[0];
  525. Matrix<Real_t>& output_data=*setup_data[0+MAX_DEPTH*2].output_data;
  526. size_t n=output_data.Dim(0)*output_data.Dim(1);
  527. Real_t* host_ptr=output_data[0];
  528. output_data.Device2HostWait();
  529. #pragma omp parallel for
  530. for(size_t i=0;i<n;i++){
  531. host_ptr[i]+=dev_ptr[i];
  532. }
  533. }
  534. Profile::Toc();
  535. Profile::Tic("Device2Host:Trg",this->Comm(),false,5);
  536. if(device) if(setup_data[0+MAX_DEPTH*0].output_data!=NULL){ // Device2Host: Target
  537. Matrix<Real_t>& output_data=*setup_data[0+MAX_DEPTH*0].output_data;
  538. assert(fmm_mat->dev_buffer.Dim()>=output_data.Dim(0)*output_data.Dim(1));
  539. output_data.Device2Host((Real_t*)&fmm_mat->dev_buffer[0]);
  540. }
  541. Profile::Toc();
  542. #endif
  543. Profile::Tic("D2D",this->Comm(),false,5);
  544. for(size_t i=0; i<=max_depth; i++){ // Down2Down
  545. if(!fmm_mat->Homogen()) fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*4],/*device*/ false);
  546. fmm_mat->Down2Down(setup_data[i+MAX_DEPTH*4]);
  547. }
  548. Profile::Toc();
  549. Profile::Tic("D2T",this->Comm(),false,5);
  550. for(int i=0; i<=(fmm_mat->Homogen()?0:max_depth); i++){ // Down2Target
  551. if(!fmm_mat->Homogen()) fmm_mat->SetupPrecomp(setup_data[i+MAX_DEPTH*5],/*device*/ false);
  552. fmm_mat->Down2Target(setup_data[i+MAX_DEPTH*5]);
  553. }
  554. Profile::Toc();
  555. #ifdef __INTEL_OFFLOAD
  556. Profile::Tic("D2H_Wait:Trg",this->Comm(),false,5);
  557. if(device) if(setup_data[0+MAX_DEPTH*0].output_data!=NULL){
  558. Real_t* dev_ptr=(Real_t*)&fmm_mat->dev_buffer[0];
  559. Matrix<Real_t>& output_data=*setup_data[0+MAX_DEPTH*0].output_data;
  560. size_t n=output_data.Dim(0)*output_data.Dim(1);
  561. Real_t* host_ptr=output_data[0];
  562. output_data.Device2HostWait();
  563. #pragma omp parallel for
  564. for(size_t i=0;i<n;i++){
  565. host_ptr[i]+=dev_ptr[i];
  566. }
  567. }
  568. Profile::Toc();
  569. #endif
  570. Profile::Tic("PostProc",this->Comm(),false,5);
  571. fmm_mat->PostProcessing(leaf_nodes);
  572. Profile::Toc();
  573. }
  574. template <class FMM_Mat_t>
  575. void FMM_Tree<FMM_Mat_t>::Copy_FMMOutput() {
  576. std::vector<FMM_Node_t*>& all_nodes=this->GetNodeList();
  577. int omp_p=omp_get_max_threads();
  578. // Copy output to the tree.
  579. {
  580. size_t k=all_nodes.size();
  581. #pragma omp parallel for
  582. for(int j=0;j<omp_p;j++){
  583. size_t a=(k*j)/omp_p;
  584. size_t b=(k*(j+1))/omp_p;
  585. fmm_mat->CopyOutput(&(all_nodes[a]),b-a);
  586. }
  587. }
  588. }
  589. }//end namespace