fmm_pts.txx 132 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310331133123313331433153316331733183319332033213322332333243325332633273328332933303331333233333334333533363337333833393340334133423343334433453346334733483349335033513352335333543355335633573358335933603361
  1. /**
  2. * \file fmm_pts.txx
  3. * \author Dhairya Malhotra, dhairya.malhotra@gmail.com
  4. * \date 3-07-2011
  5. * \brief This file contains the implementation of the FMM_Pts class.
  6. */
  7. #include <mpi.h>
  8. #include <set>
  9. #include <sstream>
  10. #include <fft_wrapper.hpp>
  11. #include <mat_utils.hpp>
  12. #ifdef PVFMM_HAVE_SYS_STAT_H
  13. #include <sys/stat.h>
  14. #endif
  15. #ifdef __SSE__
  16. #include <xmmintrin.h>
  17. #endif
  18. #ifdef __SSE3__
  19. #include <pmmintrin.h>
  20. #endif
  21. #ifdef __AVX__
  22. #include <immintrin.h>
  23. #endif
  24. #ifdef __INTEL_OFFLOAD
  25. #include <immintrin.h>
  26. #endif
  27. #ifdef __INTEL_OFFLOAD0
  28. #pragma offload_attribute(push,target(mic))
  29. #endif
  30. namespace pvfmm{
  31. /**
  32. * \brief Returns the coordinates of points on the surface of a cube.
  33. * \param[in] p Number of points on an edge of the cube is (n+1)
  34. * \param[in] c Coordinates to the centre of the cube (3D array).
  35. * \param[in] alpha Scaling factor for the size of the cube.
  36. * \param[in] depth Depth of the cube in the octree.
  37. * \return Vector with coordinates of points on the surface of the cube in the
  38. * format [x0 y0 z0 x1 y1 z1 .... ].
  39. */
  40. template <class Real_t>
  41. std::vector<Real_t> surface(int p, Real_t* c, Real_t alpha, int depth){
  42. size_t n_=(6*(p-1)*(p-1)+2); //Total number of points.
  43. std::vector<Real_t> coord(n_*3);
  44. coord[0]=coord[1]=coord[2]=-1.0;
  45. size_t cnt=1;
  46. for(int i=0;i<p-1;i++)
  47. for(int j=0;j<p-1;j++){
  48. coord[cnt*3 ]=-1.0;
  49. coord[cnt*3+1]=(2.0*(i+1)-p+1)/(p-1);
  50. coord[cnt*3+2]=(2.0*j-p+1)/(p-1);
  51. cnt++;
  52. }
  53. for(int i=0;i<p-1;i++)
  54. for(int j=0;j<p-1;j++){
  55. coord[cnt*3 ]=(2.0*i-p+1)/(p-1);
  56. coord[cnt*3+1]=-1.0;
  57. coord[cnt*3+2]=(2.0*(j+1)-p+1)/(p-1);
  58. cnt++;
  59. }
  60. for(int i=0;i<p-1;i++)
  61. for(int j=0;j<p-1;j++){
  62. coord[cnt*3 ]=(2.0*(i+1)-p+1)/(p-1);
  63. coord[cnt*3+1]=(2.0*j-p+1)/(p-1);
  64. coord[cnt*3+2]=-1.0;
  65. cnt++;
  66. }
  67. for(size_t i=0;i<(n_/2)*3;i++)
  68. coord[cnt*3+i]=-coord[i];
  69. Real_t r = 0.5*pow(0.5,depth);
  70. Real_t b = alpha*r;
  71. for(size_t i=0;i<n_;i++){
  72. coord[i*3+0]=(coord[i*3+0]+1.0)*b+c[0];
  73. coord[i*3+1]=(coord[i*3+1]+1.0)*b+c[1];
  74. coord[i*3+2]=(coord[i*3+2]+1.0)*b+c[2];
  75. }
  76. return coord;
  77. }
  78. /**
  79. * \brief Returns the coordinates of points on the upward check surface of cube.
  80. * \see surface()
  81. */
  82. template <class Real_t>
  83. std::vector<Real_t> u_check_surf(int p, Real_t* c, int depth){
  84. Real_t r=0.5*pow(0.5,depth);
  85. Real_t coord[3]={c[0]-r*(RAD1-1.0),c[1]-r*(RAD1-1.0),c[2]-r*(RAD1-1.0)};
  86. return surface(p,coord,(Real_t)RAD1,depth);
  87. }
  88. /**
  89. * \brief Returns the coordinates of points on the upward equivalent surface of cube.
  90. * \see surface()
  91. */
  92. template <class Real_t>
  93. std::vector<Real_t> u_equiv_surf(int p, Real_t* c, int depth){
  94. Real_t r=0.5*pow(0.5,depth);
  95. Real_t coord[3]={c[0]-r*(RAD0-1.0),c[1]-r*(RAD0-1.0),c[2]-r*(RAD0-1.0)};
  96. return surface(p,coord,(Real_t)RAD0,depth);
  97. }
  98. /**
  99. * \brief Returns the coordinates of points on the downward check surface of cube.
  100. * \see surface()
  101. */
  102. template <class Real_t>
  103. std::vector<Real_t> d_check_surf(int p, Real_t* c, int depth){
  104. Real_t r=0.5*pow(0.5,depth);
  105. Real_t coord[3]={c[0]-r*(RAD0-1.0),c[1]-r*(RAD0-1.0),c[2]-r*(RAD0-1.0)};
  106. return surface(p,coord,(Real_t)RAD0,depth);
  107. }
  108. /**
  109. * \brief Returns the coordinates of points on the downward equivalent surface of cube.
  110. * \see surface()
  111. */
  112. template <class Real_t>
  113. std::vector<Real_t> d_equiv_surf(int p, Real_t* c, int depth){
  114. Real_t r=0.5*pow(0.5,depth);
  115. Real_t coord[3]={c[0]-r*(RAD1-1.0),c[1]-r*(RAD1-1.0),c[2]-r*(RAD1-1.0)};
  116. return surface(p,coord,(Real_t)RAD1,depth);
  117. }
  118. /**
  119. * \brief Defines the 3D grid for convolution in FFT acceleration of V-list.
  120. * \see surface()
  121. */
  122. template <class Real_t>
  123. std::vector<Real_t> conv_grid(int p, Real_t* c, int depth){
  124. Real_t r=pow(0.5,depth);
  125. Real_t a=r*RAD0;
  126. Real_t coord[3]={c[0],c[1],c[2]};
  127. int n1=p*2;
  128. int n2=(int)pow((Real_t)n1,2);
  129. int n3=(int)pow((Real_t)n1,3);
  130. std::vector<Real_t> grid(n3*3);
  131. for(int i=0;i<n1;i++)
  132. for(int j=0;j<n1;j++)
  133. for(int k=0;k<n1;k++){
  134. grid[(i+n1*j+n2*k)*3+0]=(i-p)*a/(p-1)+coord[0];
  135. grid[(i+n1*j+n2*k)*3+1]=(j-p)*a/(p-1)+coord[1];
  136. grid[(i+n1*j+n2*k)*3+2]=(k-p)*a/(p-1)+coord[2];
  137. }
  138. return grid;
  139. }
  140. #ifdef __INTEL_OFFLOAD0
  141. #pragma offload_attribute(pop)
  142. #endif
  143. template <class Real_t>
  144. void FMM_Data<Real_t>::Clear(){
  145. upward_equiv.Resize(0);
  146. }
  147. template <class Real_t>
  148. PackedData FMM_Data<Real_t>::PackMultipole(void* buff_ptr){
  149. PackedData p0; p0.data=buff_ptr;
  150. p0.length=upward_equiv.Dim()*sizeof(Real_t);
  151. if(p0.length==0) return p0;
  152. if(p0.data==NULL) p0.data=(char*)&upward_equiv[0];
  153. else mem::memcopy(p0.data,&upward_equiv[0],p0.length);
  154. return p0;
  155. }
  156. template <class Real_t>
  157. void FMM_Data<Real_t>::AddMultipole(PackedData p0){
  158. Real_t* data=(Real_t*)p0.data;
  159. size_t n=p0.length/sizeof(Real_t);
  160. assert(upward_equiv.Dim()==n);
  161. Matrix<Real_t> v0(1,n,&upward_equiv[0],false);
  162. Matrix<Real_t> v1(1,n,data,false);
  163. v0+=v1;
  164. }
  165. template <class Real_t>
  166. void FMM_Data<Real_t>::InitMultipole(PackedData p0, bool own_data){
  167. Real_t* data=(Real_t*)p0.data;
  168. size_t n=p0.length/sizeof(Real_t);
  169. if(n==0) return;
  170. if(own_data){
  171. upward_equiv=Vector<Real_t>(n, &data[0], false);
  172. }else{
  173. upward_equiv.ReInit(n, &data[0], false);
  174. }
  175. }
  176. template <class FMMNode>
  177. FMM_Pts<FMMNode>::~FMM_Pts() {
  178. if(mat!=NULL){
  179. // int rank;
  180. // MPI_Comm_rank(comm,&rank);
  181. // if(rank==0) mat->Save2File("Precomp.data");
  182. delete mat;
  183. mat=NULL;
  184. }
  185. if(vprecomp_fft_flag) FFTW_t<Real_t>::fft_destroy_plan(vprecomp_fftplan);
  186. #ifdef __INTEL_OFFLOAD0
  187. #pragma offload target(mic:0)
  188. #endif
  189. {
  190. if(vlist_fft_flag ) FFTW_t<Real_t>::fft_destroy_plan(vlist_fftplan );
  191. if(vlist_ifft_flag) FFTW_t<Real_t>::fft_destroy_plan(vlist_ifftplan);
  192. vlist_fft_flag =false;
  193. vlist_ifft_flag=false;
  194. }
  195. }
  196. template <class FMMNode>
  197. void FMM_Pts<FMMNode>::Initialize(int mult_order, const MPI_Comm& comm_, const Kernel<Real_t>* kernel_, const Kernel<Real_t>* aux_kernel_){
  198. Profile::Tic("InitFMM_Pts",&comm_,true);{
  199. multipole_order=mult_order;
  200. comm=comm_;
  201. kernel=*kernel_;
  202. aux_kernel=*(aux_kernel_?aux_kernel_:kernel_);
  203. assert(kernel.ker_dim[0]==aux_kernel.ker_dim[0]);
  204. mat=new PrecompMat<Real_t>(Homogen(), MAX_DEPTH+1);
  205. if(this->mat_fname.size()==0){
  206. std::stringstream st;
  207. st<<PVFMM_PRECOMP_DATA_PATH;
  208. if(!st.str().size()){ // look in PVFMM_DIR
  209. char* pvfmm_dir = getenv ("PVFMM_DIR");
  210. if(pvfmm_dir) st<<pvfmm_dir<<'/';
  211. }
  212. #ifndef STAT_MACROS_BROKEN
  213. if(st.str().size()){ // check if the path is a directory
  214. struct stat stat_buff;
  215. if(stat(st.str().c_str(), &stat_buff) || !S_ISDIR(stat_buff.st_mode)){
  216. std::cout<<"error: path not found: "<<st.str()<<'\n';
  217. exit(0);
  218. }
  219. }
  220. #endif
  221. st<<"Precomp_"<<kernel.ker_name.c_str()<<"_m"<<mult_order<<(typeid(Real_t)==typeid(float)?"_f":"")<<".data";
  222. this->mat_fname=st.str();
  223. }
  224. this->mat->LoadFile(mat_fname.c_str(), this->comm);
  225. interac_list.Initialize(COORD_DIM, this->mat);
  226. Profile::Tic("PrecompUC2UE",&comm,false,4);
  227. this->PrecompAll(UC2UE_Type);
  228. Profile::Toc();
  229. Profile::Tic("PrecompDC2DE",&comm,false,4);
  230. this->PrecompAll(DC2DE_Type);
  231. Profile::Toc();
  232. Profile::Tic("PrecompBC",&comm,false,4);
  233. { /*
  234. int type=BC_Type;
  235. for(int l=0;l<MAX_DEPTH;l++)
  236. for(size_t indx=0;indx<this->interac_list.ListCount((Mat_Type)type);indx++){
  237. Matrix<Real_t>& M=this->mat->Mat(l, (Mat_Type)type, indx);
  238. M.Resize(0,0);
  239. } // */
  240. }
  241. this->PrecompAll(BC_Type,0);
  242. Profile::Toc();
  243. Profile::Tic("PrecompU2U",&comm,false,4);
  244. this->PrecompAll(U2U_Type);
  245. Profile::Toc();
  246. Profile::Tic("PrecompD2D",&comm,false,4);
  247. this->PrecompAll(D2D_Type);
  248. Profile::Toc();
  249. Profile::Tic("PrecompV",&comm,false,4);
  250. this->PrecompAll(V_Type);
  251. Profile::Toc();
  252. Profile::Tic("PrecompV1",&comm,false,4);
  253. this->PrecompAll(V1_Type);
  254. Profile::Toc();
  255. }Profile::Toc();
  256. }
  257. template <class FMMNode>
  258. Permutation<typename FMMNode::Real_t>& FMM_Pts<FMMNode>::PrecompPerm(Mat_Type type, Perm_Type perm_indx){
  259. //Check if the matrix already exists.
  260. Permutation<Real_t>& P_ = mat->Perm((Mat_Type)type, perm_indx);
  261. if(P_.Dim()!=0) return P_;
  262. Matrix<size_t> swap_xy(10,9);
  263. Matrix<size_t> swap_xz(10,9);
  264. {
  265. for(int i=0;i<9;i++)
  266. for(int j=0;j<9;j++){
  267. swap_xy[i][j]=j;
  268. swap_xz[i][j]=j;
  269. }
  270. swap_xy[3][0]=1; swap_xy[3][1]=0; swap_xy[3][2]=2;
  271. swap_xz[3][0]=2; swap_xz[3][1]=1; swap_xz[3][2]=0;
  272. swap_xy[6][0]=1; swap_xy[6][1]=0; swap_xy[6][2]=2;
  273. swap_xy[6][3]=4; swap_xy[6][4]=3; swap_xy[6][5]=5;
  274. swap_xz[6][0]=2; swap_xz[6][1]=1; swap_xz[6][2]=0;
  275. swap_xz[6][3]=5; swap_xz[6][4]=4; swap_xz[6][5]=3;
  276. swap_xy[9][0]=4; swap_xy[9][1]=3; swap_xy[9][2]=5;
  277. swap_xy[9][3]=1; swap_xy[9][4]=0; swap_xy[9][5]=2;
  278. swap_xy[9][6]=7; swap_xy[9][7]=6; swap_xy[9][8]=8;
  279. swap_xz[9][0]=8; swap_xz[9][1]=7; swap_xz[9][2]=6;
  280. swap_xz[9][3]=5; swap_xz[9][4]=4; swap_xz[9][5]=3;
  281. swap_xz[9][6]=2; swap_xz[9][7]=1; swap_xz[9][8]=0;
  282. }
  283. //Compute the matrix.
  284. Permutation<Real_t> P;
  285. switch (type){
  286. case UC2UE_Type:
  287. {
  288. break;
  289. }
  290. case DC2DE_Type:
  291. {
  292. break;
  293. }
  294. case S2U_Type:
  295. {
  296. break;
  297. }
  298. case U2U_Type:
  299. {
  300. P=PrecompPerm(D2D_Type, perm_indx);
  301. break;
  302. }
  303. case D2D_Type:
  304. {
  305. Real_t eps=1e-10;
  306. int dof=kernel.ker_dim[0];
  307. size_t p_indx=perm_indx % C_Perm;
  308. Real_t c[3]={-0.5,-0.5,-0.5};
  309. std::vector<Real_t> trg_coord=d_check_surf(this->MultipoleOrder(),c,0);
  310. int n_trg=trg_coord.size()/3;
  311. P=Permutation<Real_t>(n_trg*dof);
  312. if(p_indx==ReflecX || p_indx==ReflecY || p_indx==ReflecZ){
  313. for(int i=0;i<n_trg;i++)
  314. for(int j=0;j<n_trg;j++){
  315. if(fabs(trg_coord[i*3+0]-trg_coord[j*3+0]*(p_indx==ReflecX?-1.0:1.0))<eps)
  316. if(fabs(trg_coord[i*3+1]-trg_coord[j*3+1]*(p_indx==ReflecY?-1.0:1.0))<eps)
  317. if(fabs(trg_coord[i*3+2]-trg_coord[j*3+2]*(p_indx==ReflecZ?-1.0:1.0))<eps){
  318. for(int k=0;k<dof;k++){
  319. P.perm[j*dof+k]=i*dof+k;
  320. }
  321. }
  322. }
  323. if(dof==3) //stokes_vel (and like kernels)
  324. for(int j=0;j<n_trg;j++)
  325. P.scal[j*dof+(int)p_indx]*=-1.0;
  326. }else if(p_indx==SwapXY || p_indx==SwapXZ)
  327. for(int i=0;i<n_trg;i++)
  328. for(int j=0;j<n_trg;j++){
  329. if(fabs(trg_coord[i*3+0]-trg_coord[j*3+(p_indx==SwapXY?1:2)])<eps)
  330. if(fabs(trg_coord[i*3+1]-trg_coord[j*3+(p_indx==SwapXY?0:1)])<eps)
  331. if(fabs(trg_coord[i*3+2]-trg_coord[j*3+(p_indx==SwapXY?2:0)])<eps){
  332. for(int k=0;k<dof;k++){
  333. P.perm[j*dof+k]=i*dof+(p_indx==SwapXY?swap_xy[dof][k]:swap_xz[dof][k]);
  334. }
  335. }
  336. }
  337. break;
  338. }
  339. case D2T_Type:
  340. {
  341. break;
  342. }
  343. case U0_Type:
  344. {
  345. break;
  346. }
  347. case U1_Type:
  348. {
  349. break;
  350. }
  351. case U2_Type:
  352. {
  353. break;
  354. }
  355. case V_Type:
  356. {
  357. break;
  358. }
  359. case V1_Type:
  360. {
  361. break;
  362. }
  363. case W_Type:
  364. {
  365. break;
  366. }
  367. case X_Type:
  368. {
  369. break;
  370. }
  371. case BC_Type:
  372. {
  373. break;
  374. }
  375. default:
  376. break;
  377. }
  378. //Save the matrix for future use.
  379. #pragma omp critical (PRECOMP_MATRIX_PTS)
  380. {
  381. if(P_.Dim()==0) P_=P;
  382. }
  383. return P_;
  384. }
  385. template <class FMMNode>
  386. Matrix<typename FMMNode::Real_t>& FMM_Pts<FMMNode>::Precomp(int level, Mat_Type type, size_t mat_indx){
  387. if(this->Homogen()) level=0;
  388. //Check if the matrix already exists.
  389. Matrix<Real_t>& M_ = this->mat->Mat(level, type, mat_indx);
  390. if(M_.Dim(0)!=0 && M_.Dim(1)!=0) return M_;
  391. else{ //Compute matrix from symmetry class (if possible).
  392. size_t class_indx = this->interac_list.InteracClass(type, mat_indx);
  393. if(class_indx!=mat_indx){
  394. Matrix<Real_t>& M0 = this->Precomp(level, type, class_indx);
  395. Permutation<Real_t>& Pr = this->interac_list.Perm_R(level, type, mat_indx);
  396. Permutation<Real_t>& Pc = this->interac_list.Perm_C(level, type, mat_indx);
  397. if(Pr.Dim()>0 && Pc.Dim()>0 && M0.Dim(0)>0 && M0.Dim(1)>0) return M_;
  398. }
  399. }
  400. //Compute the matrix.
  401. Matrix<Real_t> M;
  402. int* ker_dim=kernel.ker_dim;
  403. int* aux_ker_dim=aux_kernel.ker_dim;
  404. //int omp_p=omp_get_max_threads();
  405. switch (type){
  406. case UC2UE_Type:
  407. {
  408. if(MultipoleOrder()==0) break;
  409. // Coord of upward check surface
  410. Real_t c[3]={0,0,0};
  411. std::vector<Real_t> uc_coord=u_check_surf(MultipoleOrder(),c,level);
  412. size_t n_uc=uc_coord.size()/3;
  413. // Coord of upward equivalent surface
  414. std::vector<Real_t> ue_coord=u_equiv_surf(MultipoleOrder(),c,level);
  415. size_t n_ue=ue_coord.size()/3;
  416. // Evaluate potential at check surface due to equivalent surface.
  417. Matrix<Real_t> M_e2c(n_ue*aux_ker_dim[0],n_uc*aux_ker_dim[1]);
  418. Kernel<Real_t>::Eval(&ue_coord[0], n_ue,
  419. &uc_coord[0], n_uc, &(M_e2c[0][0]),
  420. aux_kernel.ker_poten, aux_ker_dim);
  421. M=M_e2c.pinv(); //check 2 equivalent
  422. break;
  423. }
  424. case DC2DE_Type:
  425. {
  426. if(MultipoleOrder()==0) break;
  427. // Coord of downward check surface
  428. Real_t c[3]={0,0,0};
  429. std::vector<Real_t> check_surf=d_check_surf(MultipoleOrder(),c,level);
  430. size_t n_ch=check_surf.size()/3;
  431. // Coord of downward equivalent surface
  432. std::vector<Real_t> equiv_surf=d_equiv_surf(MultipoleOrder(),c,level);
  433. size_t n_eq=equiv_surf.size()/3;
  434. // Evaluate potential at check surface due to equivalent surface.
  435. Matrix<Real_t> M_e2c(n_eq*aux_ker_dim[0],n_ch*aux_ker_dim[1]);
  436. Kernel<Real_t>::Eval(&equiv_surf[0], n_eq,
  437. &check_surf[0], n_ch, &(M_e2c[0][0]),
  438. aux_kernel.ker_poten,aux_ker_dim);
  439. M=M_e2c.pinv(); //check 2 equivalent
  440. break;
  441. }
  442. case S2U_Type:
  443. {
  444. break;
  445. }
  446. case U2U_Type:
  447. {
  448. if(MultipoleOrder()==0) break;
  449. // Coord of upward check surface
  450. Real_t c[3]={0,0,0};
  451. std::vector<Real_t> check_surf=u_check_surf(MultipoleOrder(),c,level);
  452. size_t n_uc=check_surf.size()/3;
  453. // Coord of child's upward equivalent surface
  454. Real_t s=pow(0.5,(level+2));
  455. int* coord=interac_list.RelativeCoord(type,mat_indx);
  456. Real_t child_coord[3]={(coord[0]+1)*s,(coord[1]+1)*s,(coord[2]+1)*s};
  457. std::vector<Real_t> equiv_surf=u_equiv_surf(MultipoleOrder(),child_coord,level+1);
  458. size_t n_ue=equiv_surf.size()/3;
  459. // Evaluate potential at check surface due to equivalent surface.
  460. Matrix<Real_t> M_ce2c(n_ue*aux_ker_dim[0],n_uc*aux_ker_dim[1]);
  461. Kernel<Real_t>::Eval(&equiv_surf[0], n_ue,
  462. &check_surf[0], n_uc, &(M_ce2c[0][0]),
  463. aux_kernel.ker_poten, aux_ker_dim);
  464. Matrix<Real_t>& M_c2e = Precomp(level, UC2UE_Type, 0);
  465. M=M_ce2c*M_c2e;
  466. break;
  467. }
  468. case D2D_Type:
  469. {
  470. if(MultipoleOrder()==0) break;
  471. // Coord of downward check surface
  472. Real_t s=pow(0.5,level+1);
  473. int* coord=interac_list.RelativeCoord(type,mat_indx);
  474. Real_t c[3]={(coord[0]+1)*s,(coord[1]+1)*s,(coord[2]+1)*s};
  475. std::vector<Real_t> check_surf=d_check_surf(MultipoleOrder(),c,level);
  476. size_t n_dc=check_surf.size()/3;
  477. // Coord of parent's downward equivalent surface
  478. Real_t parent_coord[3]={0,0,0};
  479. std::vector<Real_t> equiv_surf=d_equiv_surf(MultipoleOrder(),parent_coord,level-1);
  480. size_t n_de=equiv_surf.size()/3;
  481. // Evaluate potential at check surface due to equivalent surface.
  482. Matrix<Real_t> M_pe2c(n_de*aux_ker_dim[0],n_dc*aux_ker_dim[1]);
  483. Kernel<Real_t>::Eval(&equiv_surf[0], n_de,
  484. &check_surf[0], n_dc, &(M_pe2c[0][0]),
  485. aux_kernel.ker_poten,aux_ker_dim);
  486. Matrix<Real_t>& M_c2e=Precomp(level,DC2DE_Type,0);
  487. M=M_pe2c*M_c2e;
  488. break;
  489. }
  490. case D2T_Type:
  491. {
  492. if(MultipoleOrder()==0) break;
  493. std::vector<Real_t>& rel_trg_coord=mat->RelativeTrgCoord();
  494. // Coord of target points
  495. Real_t r=pow(0.5,level);
  496. size_t n_trg=rel_trg_coord.size()/3;
  497. std::vector<Real_t> trg_coord(n_trg*3);
  498. for(size_t i=0;i<n_trg*COORD_DIM;i++) trg_coord[i]=rel_trg_coord[i]*r;
  499. // Coord of downward equivalent surface
  500. Real_t c[3]={0,0,0};
  501. std::vector<Real_t> equiv_surf=d_equiv_surf(MultipoleOrder(),c,level);
  502. size_t n_eq=equiv_surf.size()/3;
  503. // Evaluate potential at target points due to equivalent surface.
  504. {
  505. M .Resize(n_eq*ker_dim [0], n_trg*ker_dim [1]);
  506. kernel.BuildMatrix(&equiv_surf[0], n_eq, &trg_coord[0], n_trg, &(M [0][0]));
  507. }
  508. break;
  509. }
  510. case U0_Type:
  511. {
  512. break;
  513. }
  514. case U1_Type:
  515. {
  516. break;
  517. }
  518. case U2_Type:
  519. {
  520. break;
  521. }
  522. case V_Type:
  523. {
  524. if(MultipoleOrder()==0) break;
  525. int n1=MultipoleOrder()*2;
  526. int n3 =n1*n1*n1;
  527. int n3_=n1*n1*(n1/2+1);
  528. //Compute the matrix.
  529. Real_t s=pow(0.5,level);
  530. int* coord2=interac_list.RelativeCoord(type,mat_indx);
  531. Real_t coord_diff[3]={coord2[0]*s,coord2[1]*s,coord2[2]*s};
  532. //Evaluate potential.
  533. std::vector<Real_t> r_trg(COORD_DIM,0.0);
  534. std::vector<Real_t> conv_poten(n3*aux_ker_dim[0]*aux_ker_dim[1]);
  535. std::vector<Real_t> conv_coord=conv_grid(MultipoleOrder(),coord_diff,level);
  536. Kernel<Real_t>::Eval(&conv_coord[0],n3,&r_trg[0],1,&conv_poten[0],aux_kernel.ker_poten,aux_ker_dim);
  537. //Rearrange data.
  538. Matrix<Real_t> M_conv(n3,aux_ker_dim[0]*aux_ker_dim[1],&conv_poten[0],false);
  539. M_conv=M_conv.Transpose();
  540. //Compute FFTW plan.
  541. int nnn[3]={n1,n1,n1};
  542. void *fftw_in, *fftw_out;
  543. fftw_in = fftw_malloc( n3 *aux_ker_dim[0]*aux_ker_dim[1]*sizeof(Real_t));
  544. fftw_out = fftw_malloc(2*n3_*aux_ker_dim[0]*aux_ker_dim[1]*sizeof(Real_t));
  545. #pragma omp critical (FFTW_PLAN)
  546. {
  547. if (!vprecomp_fft_flag){
  548. vprecomp_fftplan = FFTW_t<Real_t>::fft_plan_many_dft_r2c(COORD_DIM, nnn, aux_ker_dim[0]*aux_ker_dim[1],
  549. (Real_t*)fftw_in, NULL, 1, n3, (typename FFTW_t<Real_t>::cplx*) fftw_out, NULL, 1, n3_, FFTW_ESTIMATE);
  550. vprecomp_fft_flag=true;
  551. }
  552. }
  553. //Compute FFT.
  554. mem::memcopy(fftw_in, &conv_poten[0], n3*aux_ker_dim[0]*aux_ker_dim[1]*sizeof(Real_t));
  555. FFTW_t<Real_t>::fft_execute_dft_r2c(vprecomp_fftplan, (Real_t*)fftw_in, (typename FFTW_t<Real_t>::cplx*)(fftw_out));
  556. Matrix<Real_t> M_(2*n3_*aux_ker_dim[0]*aux_ker_dim[1],1,(Real_t*)fftw_out,false);
  557. M=M_;
  558. //Free memory.
  559. fftw_free(fftw_in);
  560. fftw_free(fftw_out);
  561. break;
  562. }
  563. case V1_Type:
  564. {
  565. if(MultipoleOrder()==0) break;
  566. size_t mat_cnt =interac_list.ListCount( V_Type);
  567. for(size_t k=0;k<mat_cnt;k++) Precomp(level, V_Type, k);
  568. const size_t chld_cnt=1UL<<COORD_DIM;
  569. size_t n1=MultipoleOrder()*2;
  570. size_t M_dim=n1*n1*(n1/2+1);
  571. size_t n3=n1*n1*n1;
  572. Vector<Real_t> zero_vec(M_dim*aux_ker_dim[0]*aux_ker_dim[1]*2);
  573. zero_vec.SetZero();
  574. Vector<Real_t*> M_ptr(chld_cnt*chld_cnt);
  575. for(size_t i=0;i<chld_cnt*chld_cnt;i++) M_ptr[i]=&zero_vec[0];
  576. int* rel_coord_=interac_list.RelativeCoord(V1_Type, mat_indx);
  577. for(int j1=0;j1<chld_cnt;j1++)
  578. for(int j2=0;j2<chld_cnt;j2++){
  579. int rel_coord[3]={rel_coord_[0]*2-(j1/1)%2+(j2/1)%2,
  580. rel_coord_[1]*2-(j1/2)%2+(j2/2)%2,
  581. rel_coord_[2]*2-(j1/4)%2+(j2/4)%2};
  582. for(size_t k=0;k<mat_cnt;k++){
  583. int* ref_coord=interac_list.RelativeCoord(V_Type, k);
  584. if(ref_coord[0]==rel_coord[0] &&
  585. ref_coord[1]==rel_coord[1] &&
  586. ref_coord[2]==rel_coord[2]){
  587. Matrix<Real_t>& M = this->mat->Mat(level, V_Type, k);
  588. M_ptr[j2*chld_cnt+j1]=&M[0][0];
  589. break;
  590. }
  591. }
  592. }
  593. // Build matrix aux_ker_dim0 x aux_ker_dim1 x M_dim x 8 x 8
  594. M.Resize(aux_ker_dim[0]*aux_ker_dim[1]*M_dim, 2*chld_cnt*chld_cnt);
  595. for(int j=0;j<aux_ker_dim[0]*aux_ker_dim[1]*M_dim;j++){
  596. for(size_t k=0;k<chld_cnt*chld_cnt;k++){
  597. M[j][k*2+0]=M_ptr[k][j*2+0]/n3;
  598. M[j][k*2+1]=M_ptr[k][j*2+1]/n3;
  599. }
  600. }
  601. break;
  602. }
  603. case W_Type:
  604. {
  605. if(MultipoleOrder()==0) break;
  606. std::vector<Real_t>& rel_trg_coord=mat->RelativeTrgCoord();
  607. // Coord of target points
  608. Real_t s=pow(0.5,level);
  609. size_t n_trg=rel_trg_coord.size()/3;
  610. std::vector<Real_t> trg_coord(n_trg*3);
  611. for(size_t j=0;j<n_trg*COORD_DIM;j++) trg_coord[j]=rel_trg_coord[j]*s;
  612. // Coord of downward equivalent surface
  613. int* coord2=interac_list.RelativeCoord(type,mat_indx);
  614. Real_t c[3]={(coord2[0]+1)*s*0.25,(coord2[1]+1)*s*0.25,(coord2[2]+1)*s*0.25};
  615. std::vector<Real_t> equiv_surf=u_equiv_surf(MultipoleOrder(),c,level+1);
  616. size_t n_eq=equiv_surf.size()/3;
  617. // Evaluate potential at target points due to equivalent surface.
  618. {
  619. M .Resize(n_eq*ker_dim [0],n_trg*ker_dim [1]);
  620. kernel.BuildMatrix(&equiv_surf[0], n_eq, &trg_coord[0], n_trg, &(M [0][0]));
  621. }
  622. break;
  623. }
  624. case X_Type:
  625. {
  626. break;
  627. }
  628. case BC_Type:
  629. {
  630. if(MultipoleOrder()==0) break;
  631. size_t mat_cnt_m2m=interac_list.ListCount(U2U_Type);
  632. size_t n_surf=(6*(MultipoleOrder()-1)*(MultipoleOrder()-1)+2); //Total number of points.
  633. if((M.Dim(0)!=n_surf*aux_ker_dim[0] || M.Dim(1)!=n_surf*aux_ker_dim[1]) && level==0){
  634. Matrix<Real_t> M_m2m[BC_LEVELS+1];
  635. Matrix<Real_t> M_m2l[BC_LEVELS+1];
  636. Matrix<Real_t> M_l2l[BC_LEVELS+1];
  637. Matrix<Real_t> M_zero_avg(n_surf*aux_ker_dim[0],n_surf*aux_ker_dim[0]);
  638. { // Set average multipole charge to zero. (improves stability for large BC_LEVELS)
  639. M_zero_avg.SetZero();
  640. for(size_t i=0;i<n_surf*aux_ker_dim[0];i++)
  641. M_zero_avg[i][i]+=1;
  642. for(size_t i=0;i<n_surf;i++)
  643. for(size_t j=0;j<n_surf;j++)
  644. for(size_t k=0;k<aux_ker_dim[0];k++)
  645. M_zero_avg[i*aux_ker_dim[0]+k][j*aux_ker_dim[0]+k]-=1.0/n_surf;
  646. }
  647. for(int level=0; level>-BC_LEVELS; level--){
  648. M_l2l[-level] = this->Precomp(level, D2D_Type, 0);
  649. if(M_l2l[-level].Dim(0)==0 || M_l2l[-level].Dim(1)==0){
  650. Matrix<Real_t>& M0 = interac_list.ClassMat(level, D2D_Type, 0);
  651. Permutation<Real_t>& Pr = this->interac_list.Perm_R(level, D2D_Type, 0);
  652. Permutation<Real_t>& Pc = this->interac_list.Perm_C(level, D2D_Type, 0);
  653. M_l2l[-level] = Pr*M0*Pc;
  654. }
  655. M_m2m[-level] = M_zero_avg*this->Precomp(level, U2U_Type, 0);
  656. for(size_t mat_indx=1; mat_indx<mat_cnt_m2m; mat_indx++)
  657. M_m2m[-level] += M_zero_avg*this->Precomp(level, U2U_Type, mat_indx);
  658. }
  659. for(int level=-BC_LEVELS;level<=0;level++){
  660. if(!Homogen() || level==-BC_LEVELS){
  661. Real_t s=(1UL<<(-level));
  662. Real_t ue_coord[3]={0,0,0};
  663. Real_t dc_coord[3]={0,0,0};
  664. std::vector<Real_t> src_coord=u_equiv_surf(MultipoleOrder(), ue_coord, level);
  665. std::vector<Real_t> trg_coord=d_check_surf(MultipoleOrder(), dc_coord, level);
  666. Matrix<Real_t> M_ue2dc(n_surf*aux_ker_dim[0], n_surf*aux_ker_dim[1]);
  667. M_ue2dc.SetZero();
  668. for(int x0=-2;x0<4;x0++)
  669. for(int x1=-2;x1<4;x1++)
  670. for(int x2=-2;x2<4;x2++)
  671. if(abs(x0)>1 || abs(x1)>1 || abs(x2)>1){
  672. ue_coord[0]=x0*s; ue_coord[1]=x1*s; ue_coord[2]=x2*s;
  673. std::vector<Real_t> src_coord=u_equiv_surf(MultipoleOrder(), ue_coord, level);
  674. Matrix<Real_t> M_tmp(n_surf*aux_ker_dim[0], n_surf*aux_ker_dim[1]);
  675. Kernel<Real_t>::Eval(&src_coord[0], n_surf,
  676. &trg_coord[0], n_surf, &(M_tmp[0][0]),
  677. aux_kernel.ker_poten, aux_ker_dim);
  678. M_ue2dc+=M_tmp;
  679. }
  680. // Shift by constant.
  681. Real_t scale_adj=(Homogen()?pow(0.5, level*aux_kernel.poten_scale):1);
  682. for(size_t i=0;i<M_ue2dc.Dim(0);i++){
  683. std::vector<Real_t> avg(aux_ker_dim[1],0);
  684. for(size_t j=0; j<M_ue2dc.Dim(1); j+=aux_ker_dim[1])
  685. for(int k=0; k<aux_ker_dim[1]; k++) avg[k]+=M_ue2dc[i][j+k];
  686. for(int k=0; k<aux_ker_dim[1]; k++) avg[k]/=n_surf;
  687. for(size_t j=0; j<M_ue2dc.Dim(1); j+=aux_ker_dim[1])
  688. for(int k=0; k<aux_ker_dim[1]; k++) M_ue2dc[i][j+k]=(M_ue2dc[i][j+k]-avg[k])*scale_adj;
  689. }
  690. Matrix<Real_t>& M_dc2de = Precomp(level, DC2DE_Type, 0);
  691. M_m2l[-level]=M_ue2dc*M_dc2de;
  692. }else M_m2l[-level]=M_m2l[1-level];
  693. if(level==-BC_LEVELS) M = M_m2l[-level];
  694. else M = M_m2l[-level] + M_m2m[-level]*M*M_l2l[-level];
  695. { // Shift by constant. (improves stability for large BC_LEVELS)
  696. Matrix<Real_t> M_de2dc(n_surf*aux_ker_dim[0], n_surf*aux_ker_dim[1]);
  697. { // M_de2dc TODO: For homogeneous kernels, compute only once.
  698. // Coord of downward check surface
  699. Real_t c[3]={0,0,0};
  700. std::vector<Real_t> check_surf=d_check_surf(MultipoleOrder(),c,0);
  701. size_t n_ch=check_surf.size()/3;
  702. // Coord of downward equivalent surface
  703. std::vector<Real_t> equiv_surf=d_equiv_surf(MultipoleOrder(),c,0);
  704. size_t n_eq=equiv_surf.size()/3;
  705. // Evaluate potential at check surface due to equivalent surface.
  706. Kernel<Real_t>::Eval(&equiv_surf[0], n_eq,
  707. &check_surf[0], n_ch, &(M_de2dc[0][0]),
  708. aux_kernel.ker_poten,aux_ker_dim);
  709. }
  710. Matrix<Real_t> M_ue2dc=M*M_de2dc;
  711. for(size_t i=0;i<M_ue2dc.Dim(0);i++){
  712. std::vector<Real_t> avg(aux_ker_dim[1],0);
  713. for(size_t j=0; j<M_ue2dc.Dim(1); j+=aux_ker_dim[1])
  714. for(int k=0; k<aux_ker_dim[1]; k++) avg[k]+=M_ue2dc[i][j+k];
  715. for(int k=0; k<aux_ker_dim[1]; k++) avg[k]/=n_surf;
  716. for(size_t j=0; j<M_ue2dc.Dim(1); j+=aux_ker_dim[1])
  717. for(int k=0; k<aux_ker_dim[1]; k++) M_ue2dc[i][j+k]-=avg[k];
  718. }
  719. Matrix<Real_t>& M_dc2de = Precomp(level, DC2DE_Type, 0);
  720. M=M_ue2dc*M_dc2de;
  721. }
  722. }
  723. { // ax+by+cz+d correction.
  724. std::vector<Real_t> corner_pts;
  725. corner_pts.push_back(0); corner_pts.push_back(0); corner_pts.push_back(0);
  726. corner_pts.push_back(1); corner_pts.push_back(0); corner_pts.push_back(0);
  727. corner_pts.push_back(0); corner_pts.push_back(1); corner_pts.push_back(0);
  728. corner_pts.push_back(0); corner_pts.push_back(0); corner_pts.push_back(1);
  729. size_t n_corner=corner_pts.size()/3;
  730. // Coord of downward equivalent surface
  731. Real_t c[3]={0,0,0};
  732. std::vector<Real_t> up_equiv_surf=u_equiv_surf(MultipoleOrder(),c,0);
  733. std::vector<Real_t> dn_equiv_surf=d_equiv_surf(MultipoleOrder(),c,0);
  734. std::vector<Real_t> dn_check_surf=d_check_surf(MultipoleOrder(),c,0);
  735. Matrix<Real_t> M_err;
  736. { // Evaluate potential at corner due to upward and dnward equivalent surface.
  737. { // Error from local expansion.
  738. Matrix<Real_t> M_e2pt(n_surf*aux_ker_dim[0],n_corner*aux_ker_dim[1]);
  739. Kernel<Real_t>::Eval(&dn_equiv_surf[0], n_surf,
  740. &corner_pts[0], n_corner, &(M_e2pt[0][0]),
  741. aux_kernel.ker_poten,aux_ker_dim);
  742. M_err=M*M_e2pt;
  743. }
  744. for(size_t k=0;k<4;k++){ // Error from colleagues of root.
  745. for(int j0=-1;j0<=1;j0++)
  746. for(int j1=-1;j1<=1;j1++)
  747. for(int j2=-1;j2<=1;j2++){
  748. Real_t pt_coord[3]={corner_pts[k*COORD_DIM+0]-j0,
  749. corner_pts[k*COORD_DIM+1]-j1,
  750. corner_pts[k*COORD_DIM+2]-j2};
  751. if(fabs(pt_coord[0]-0.5)>1.0 || fabs(pt_coord[1]-0.5)>1.0 || fabs(pt_coord[2]-0.5)>1.0){
  752. Matrix<Real_t> M_e2pt(n_surf*aux_ker_dim[0],aux_ker_dim[1]);
  753. Kernel<Real_t>::Eval(&up_equiv_surf[0], n_surf,
  754. &pt_coord[0], 1, &(M_e2pt[0][0]),
  755. aux_kernel.ker_poten,aux_ker_dim);
  756. for(size_t i=0;i<M_e2pt.Dim(0);i++)
  757. for(size_t j=0;j<M_e2pt.Dim(1);j++)
  758. M_err[i][k*aux_ker_dim[1]+j]+=M_e2pt[i][j];
  759. }
  760. }
  761. }
  762. }
  763. Matrix<Real_t> M_grad(M_err.Dim(0),n_surf*aux_ker_dim[1]);
  764. for(size_t i=0;i<M_err.Dim(0);i++)
  765. for(size_t k=0;k<aux_ker_dim[1];k++)
  766. for(size_t j=0;j<n_surf;j++){
  767. M_grad[i][j*aux_ker_dim[1]+k]=(M_err[i][0*aux_ker_dim[1]+k] )*1.0 +
  768. (M_err[i][1*aux_ker_dim[1]+k]-M_err[i][0*aux_ker_dim[1]+k])*dn_check_surf[j*COORD_DIM+0]+
  769. (M_err[i][2*aux_ker_dim[1]+k]-M_err[i][0*aux_ker_dim[1]+k])*dn_check_surf[j*COORD_DIM+1]+
  770. (M_err[i][3*aux_ker_dim[1]+k]-M_err[i][0*aux_ker_dim[1]+k])*dn_check_surf[j*COORD_DIM+2];
  771. }
  772. Matrix<Real_t>& M_dc2de = Precomp(0, DC2DE_Type, 0);
  773. M-=M_grad*M_dc2de;
  774. }
  775. { // Free memory
  776. Mat_Type type=D2D_Type;
  777. for(int l=-BC_LEVELS;l<0;l++)
  778. for(size_t indx=0;indx<this->interac_list.ListCount(type);indx++){
  779. Matrix<Real_t>& M=this->mat->Mat(l, type, indx);
  780. M.Resize(0,0);
  781. }
  782. type=U2U_Type;
  783. for(int l=-BC_LEVELS;l<0;l++)
  784. for(size_t indx=0;indx<this->interac_list.ListCount(type);indx++){
  785. Matrix<Real_t>& M=this->mat->Mat(l, type, indx);
  786. M.Resize(0,0);
  787. }
  788. type=DC2DE_Type;
  789. for(int l=-BC_LEVELS;l<0;l++)
  790. for(size_t indx=0;indx<this->interac_list.ListCount(type);indx++){
  791. Matrix<Real_t>& M=this->mat->Mat(l, type, indx);
  792. M.Resize(0,0);
  793. }
  794. type=UC2UE_Type;
  795. for(int l=-BC_LEVELS;l<0;l++)
  796. for(size_t indx=0;indx<this->interac_list.ListCount(type);indx++){
  797. Matrix<Real_t>& M=this->mat->Mat(l, type, indx);
  798. M.Resize(0,0);
  799. }
  800. }
  801. }
  802. break;
  803. }
  804. default:
  805. break;
  806. }
  807. //Save the matrix for future use.
  808. #pragma omp critical (PRECOMP_MATRIX_PTS)
  809. if(M_.Dim(0)==0 && M_.Dim(1)==0){
  810. M_=M;
  811. /*
  812. M_.Resize(M.Dim(0),M.Dim(1));
  813. int dof=aux_ker_dim[0]*aux_ker_dim[1];
  814. for(int j=0;j<dof;j++){
  815. size_t a=(M.Dim(0)*M.Dim(1)* j )/dof;
  816. size_t b=(M.Dim(0)*M.Dim(1)*(j+1))/dof;
  817. #pragma omp parallel for // NUMA
  818. for(int tid=0;tid<omp_p;tid++){
  819. size_t a_=a+((b-a)* tid )/omp_p;
  820. size_t b_=a+((b-a)*(tid+1))/omp_p;
  821. mem::memcopy(&M_[0][a_], &M[0][a_], (b_-a_)*sizeof(Real_t));
  822. }
  823. }
  824. */
  825. }
  826. return M_;
  827. }
  828. template <class FMMNode>
  829. void FMM_Pts<FMMNode>::PrecompAll(Mat_Type type, int level){
  830. int depth=(Homogen()?1:MAX_DEPTH);
  831. if(level==-1){
  832. for(int l=0;l<depth;l++){
  833. std::stringstream level_str;
  834. level_str<<"level="<<l;
  835. PrecompAll(type, l);
  836. }
  837. return;
  838. }
  839. //Compute basic permutations.
  840. for(size_t i=0;i<Perm_Count;i++)
  841. this->PrecompPerm(type, (Perm_Type) i);
  842. {
  843. //Allocate matrices.
  844. size_t mat_cnt=interac_list.ListCount((Mat_Type)type);
  845. mat->Mat(level, (Mat_Type)type, mat_cnt-1);
  846. { // Compute InteracClass matrices.
  847. std::vector<size_t> indx_lst;
  848. for(size_t i=0; i<mat_cnt; i++){
  849. if(interac_list.InteracClass((Mat_Type)type,i)==i)
  850. indx_lst.push_back(i);
  851. }
  852. //Compute Transformations.
  853. //#pragma omp parallel for //lets use fine grained parallelism
  854. for(size_t i=0; i<indx_lst.size(); i++){
  855. Precomp(level, (Mat_Type)type, indx_lst[i]);
  856. }
  857. }
  858. //#pragma omp parallel for //lets use fine grained parallelism
  859. for(size_t mat_indx=0;mat_indx<mat_cnt;mat_indx++){
  860. Matrix<Real_t>& M0=interac_list.ClassMat(level,(Mat_Type)type,mat_indx);
  861. Permutation<Real_t>& pr=interac_list.Perm_R(level, (Mat_Type)type, mat_indx);
  862. Permutation<Real_t>& pc=interac_list.Perm_C(level, (Mat_Type)type, mat_indx);
  863. if(pr.Dim()!=M0.Dim(0) || pc.Dim()!=M0.Dim(1)) Precomp(level, (Mat_Type)type, mat_indx);
  864. }
  865. }
  866. }
  867. template <class FMMNode>
  868. void FMM_Pts<FMMNode>::CollectNodeData(std::vector<FMMNode*>& node, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, std::vector<size_t> extra_size){
  869. if( buff.size()<7) buff.resize(7);
  870. if( n_list.size()<7) n_list.resize(7);
  871. if(node.size()==0) return;
  872. {// 0. upward_equiv
  873. int indx=0;
  874. Matrix<Real_t>& M_uc2ue = this->interac_list.ClassMat(0, UC2UE_Type, 0);
  875. size_t vec_sz=M_uc2ue.Dim(1);
  876. std::vector< FMMNode* > node_lst;
  877. {
  878. std::vector<std::vector< FMMNode* > > node_lst_(MAX_DEPTH+1);
  879. FMMNode_t* r_node=NULL;
  880. for(size_t i=0;i<node.size();i++){
  881. if(!node[i]->IsLeaf())
  882. node_lst_[node[i]->Depth()].push_back(node[i]);
  883. if(node[i]->Depth()==0) r_node=node[i];
  884. }
  885. size_t chld_cnt=1UL<<COORD_DIM;
  886. for(int i=0;i<=MAX_DEPTH;i++)
  887. for(size_t j=0;j<node_lst_[i].size();j++)
  888. for(size_t k=0;k<chld_cnt;k++)
  889. node_lst.push_back((FMMNode_t*)node_lst_[i][j]->Child(k));
  890. if(r_node!=NULL) node_lst.push_back(r_node);
  891. }
  892. n_list[indx]=node_lst;
  893. size_t buff_size=node_lst.size()*vec_sz;
  894. buff_size+=(extra_size.size()>indx?extra_size[indx]:0);
  895. buff[indx].Resize(1,buff_size);
  896. #pragma omp parallel for
  897. for(size_t i=0;i<node_lst.size();i++){
  898. Vector<Real_t>& upward_equiv=node_lst[i]->FMMData()->upward_equiv;
  899. upward_equiv.ReInit(vec_sz, buff[indx][0]+i*vec_sz, false);
  900. upward_equiv.SetZero();
  901. }
  902. buff[indx].AllocDevice(true);
  903. }
  904. {// 1. dnward_equiv
  905. int indx=1;
  906. Matrix<Real_t>& M_dc2de = this->interac_list.ClassMat(0, DC2DE_Type, 0);
  907. size_t vec_sz=M_dc2de.Dim(1);
  908. std::vector< FMMNode* > node_lst;
  909. {
  910. std::vector<std::vector< FMMNode* > > node_lst_(MAX_DEPTH+1);
  911. FMMNode_t* r_node=NULL;
  912. for(size_t i=0;i<node.size();i++){
  913. if(!node[i]->IsLeaf() && !node[i]->IsGhost())
  914. node_lst_[node[i]->Depth()].push_back(node[i]);
  915. if(node[i]->Depth()==0) r_node=node[i];
  916. }
  917. size_t chld_cnt=1UL<<COORD_DIM;
  918. for(int i=0;i<=MAX_DEPTH;i++)
  919. for(size_t j=0;j<node_lst_[i].size();j++)
  920. for(size_t k=0;k<chld_cnt;k++)
  921. node_lst.push_back((FMMNode_t*)node_lst_[i][j]->Child(k));
  922. if(r_node!=NULL) node_lst.push_back(r_node);
  923. }
  924. n_list[indx]=node_lst;
  925. size_t buff_size=node_lst.size()*vec_sz;
  926. buff_size+=(extra_size.size()>indx?extra_size[indx]:0);
  927. buff[indx].Resize(1,buff_size);
  928. #pragma omp parallel for
  929. for(size_t i=0;i<node_lst.size();i++){
  930. Vector<Real_t>& dnward_equiv=node_lst[i]->FMMData()->dnward_equiv;
  931. dnward_equiv.ReInit(vec_sz, buff[indx][0]+i*vec_sz, false);
  932. dnward_equiv.SetZero();
  933. }
  934. buff[indx].AllocDevice(true);
  935. }
  936. {// 2. upward_equiv_fft
  937. int indx=2;
  938. std::vector< FMMNode* > node_lst;
  939. {
  940. std::vector<std::vector< FMMNode* > > node_lst_(MAX_DEPTH+1);
  941. for(size_t i=0;i<node.size();i++)
  942. if(!node[i]->IsLeaf())
  943. node_lst_[node[i]->Depth()].push_back(node[i]);
  944. for(int i=0;i<=MAX_DEPTH;i++)
  945. for(size_t j=0;j<node_lst_[i].size();j++)
  946. node_lst.push_back(node_lst_[i][j]);
  947. }
  948. n_list[indx]=node_lst;
  949. buff[indx].AllocDevice(true);
  950. }
  951. {// 3. dnward_check_fft
  952. int indx=3;
  953. std::vector< FMMNode* > node_lst;
  954. {
  955. std::vector<std::vector< FMMNode* > > node_lst_(MAX_DEPTH+1);
  956. for(size_t i=0;i<node.size();i++)
  957. if(!node[i]->IsLeaf() && !node[i]->IsGhost())
  958. node_lst_[node[i]->Depth()].push_back(node[i]);
  959. for(int i=0;i<=MAX_DEPTH;i++)
  960. for(size_t j=0;j<node_lst_[i].size();j++)
  961. node_lst.push_back(node_lst_[i][j]);
  962. }
  963. n_list[indx]=node_lst;
  964. buff[indx].AllocDevice(true);
  965. }
  966. {// 4. src_val
  967. int indx=4;
  968. int src_dof=kernel.ker_dim[0];
  969. int surf_dof=COORD_DIM+src_dof;
  970. std::vector< FMMNode* > node_lst;
  971. size_t buff_size=0;
  972. for(size_t i=0;i<node.size();i++)
  973. if(node[i]->IsLeaf()){
  974. node_lst.push_back(node[i]);
  975. buff_size+=(node[i]->src_coord.Dim()/COORD_DIM)*src_dof;
  976. buff_size+=(node[i]->surf_coord.Dim()/COORD_DIM)*surf_dof;
  977. }
  978. n_list[indx]=node_lst;
  979. #pragma omp parallel for
  980. for(size_t i=0;i<node_lst.size();i++){ // Move data before resizing buff[indx]
  981. { // src_value
  982. Vector<Real_t>& src_value=node_lst[i]->src_value;
  983. Vector<Real_t> new_buff=src_value;
  984. src_value.ReInit(new_buff.Dim(), &new_buff[0]);
  985. }
  986. { // surf_value
  987. Vector<Real_t>& surf_value=node_lst[i]->surf_value;
  988. Vector<Real_t> new_buff=surf_value;
  989. surf_value.ReInit(new_buff.Dim(), &new_buff[0]);
  990. }
  991. }
  992. buff[indx].Resize(1,buff_size+(extra_size.size()>indx?extra_size[indx]:0));
  993. Real_t* buff_ptr=&buff[indx][0][0];
  994. for(size_t i=0;i<node_lst.size();i++){
  995. FMMData* fmm_data=((FMMData*)node_lst[i]->FMMData());
  996. { // src_value
  997. Vector<Real_t>& src_value=node_lst[i]->src_value;
  998. mem::memcopy(buff_ptr,&src_value[0],src_value.Dim()*sizeof(Real_t));
  999. src_value.ReInit((node_lst[i]->src_coord.Dim()/COORD_DIM)*src_dof, buff_ptr, false);
  1000. buff_ptr+=(node_lst[i]->src_coord.Dim()/COORD_DIM)*src_dof;
  1001. }
  1002. { // surf_value
  1003. Vector<Real_t>& surf_value=node_lst[i]->surf_value;
  1004. mem::memcopy(buff_ptr,&surf_value[0],surf_value.Dim()*sizeof(Real_t));
  1005. surf_value.ReInit((node_lst[i]->surf_coord.Dim()/COORD_DIM)*surf_dof, buff_ptr, false);
  1006. buff_ptr+=(node_lst[i]->surf_coord.Dim()/COORD_DIM)*surf_dof;
  1007. }
  1008. }
  1009. buff[indx].AllocDevice(true);
  1010. }
  1011. {// 5. trg_val
  1012. int indx=5;
  1013. int trg_dof=kernel.ker_dim[1];
  1014. std::vector< FMMNode* > node_lst;
  1015. size_t buff_size=0;
  1016. for(size_t i=0;i<node.size();i++)
  1017. if(node[i]->IsLeaf() && !node[i]->IsGhost()){
  1018. node_lst.push_back(node[i]);
  1019. FMMData* fmm_data=((FMMData*)node[i]->FMMData());
  1020. buff_size+=(node[i]->trg_coord.Dim()/COORD_DIM)*trg_dof;
  1021. }
  1022. n_list[indx]=node_lst;
  1023. buff[indx].Resize(1,buff_size+(extra_size.size()>indx?extra_size[indx]:0));
  1024. Real_t* buff_ptr=&buff[indx][0][0];
  1025. for(size_t i=0;i<node_lst.size();i++){
  1026. FMMData* fmm_data=((FMMData*)node_lst[i]->FMMData());
  1027. { // trg_value
  1028. Vector<Real_t>& trg_value=node_lst[i]->trg_value;
  1029. trg_value.ReInit((node_lst[i]->trg_coord.Dim()/COORD_DIM)*trg_dof, buff_ptr, false);
  1030. buff_ptr+=(node_lst[i]->trg_coord.Dim()/COORD_DIM)*trg_dof;
  1031. }
  1032. }
  1033. #pragma omp parallel for
  1034. for(size_t i=0;i<node_lst.size();i++){
  1035. FMMData* fmm_data=((FMMData*)node_lst[i]->FMMData());
  1036. Vector<Real_t>& trg_value=node_lst[i]->trg_value;
  1037. trg_value.SetZero();
  1038. }
  1039. buff[indx].AllocDevice(true);
  1040. }
  1041. {// 6. pts_coord
  1042. int indx=6;
  1043. size_t m=MultipoleOrder();
  1044. std::vector< FMMNode* > node_lst;
  1045. size_t buff_size=0;
  1046. for(size_t i=0;i<node.size();i++)
  1047. if(node[i]->IsLeaf()){
  1048. node_lst.push_back(node[i]);
  1049. FMMData* fmm_data=((FMMData*)node[i]->FMMData());
  1050. buff_size+=node[i]->src_coord.Dim();
  1051. buff_size+=node[i]->surf_coord.Dim();
  1052. buff_size+=node[i]->trg_coord.Dim();
  1053. }
  1054. n_list[indx]=node_lst;
  1055. #pragma omp parallel for
  1056. for(size_t i=0;i<node_lst.size();i++){ // Move data before resizing buff[indx]
  1057. FMMData* fmm_data=((FMMData*)node_lst[i]->FMMData());
  1058. { // src_coord
  1059. Vector<Real_t>& src_coord=node_lst[i]->src_coord;
  1060. Vector<Real_t> new_buff=src_coord;
  1061. src_coord.ReInit(new_buff.Dim(), &new_buff[0]);
  1062. }
  1063. { // surf_coord
  1064. Vector<Real_t>& surf_coord=node_lst[i]->surf_coord;
  1065. Vector<Real_t> new_buff=surf_coord;
  1066. surf_coord.ReInit(new_buff.Dim(), &new_buff[0]);
  1067. }
  1068. { // trg_coord
  1069. Vector<Real_t>& trg_coord=node_lst[i]->trg_coord;
  1070. Vector<Real_t> new_buff=trg_coord;
  1071. trg_coord.ReInit(new_buff.Dim(), &new_buff[0]);
  1072. }
  1073. }
  1074. buff_size+=(extra_size.size()>indx?extra_size[indx]:0);
  1075. buff_size+=4*MAX_DEPTH*(6*(m-1)*(m-1)+2)*COORD_DIM;
  1076. buff[indx].Resize(1,buff_size);
  1077. Real_t* buff_ptr=&buff[indx][0][0];
  1078. for(size_t i=0;i<node_lst.size();i++){
  1079. FMMData* fmm_data=((FMMData*)node_lst[i]->FMMData());
  1080. { // src_coord
  1081. Vector<Real_t>& src_coord=node_lst[i]->src_coord;
  1082. mem::memcopy(buff_ptr,&src_coord[0],src_coord.Dim()*sizeof(Real_t));
  1083. src_coord.ReInit(node_lst[i]->src_coord.Dim(), buff_ptr, false);
  1084. buff_ptr+=node_lst[i]->src_coord.Dim();
  1085. }
  1086. { // surf_coord
  1087. Vector<Real_t>& surf_coord=node_lst[i]->surf_coord;
  1088. mem::memcopy(buff_ptr,&surf_coord[0],surf_coord.Dim()*sizeof(Real_t));
  1089. surf_coord.ReInit(node_lst[i]->surf_coord.Dim(), buff_ptr, false);
  1090. buff_ptr+=node_lst[i]->surf_coord.Dim();
  1091. }
  1092. { // trg_coord
  1093. Vector<Real_t>& trg_coord=node_lst[i]->trg_coord;
  1094. mem::memcopy(buff_ptr,&trg_coord[0],trg_coord.Dim()*sizeof(Real_t));
  1095. trg_coord.ReInit(node_lst[i]->trg_coord.Dim(), buff_ptr, false);
  1096. buff_ptr+=node_lst[i]->trg_coord.Dim();
  1097. }
  1098. }
  1099. { // check and equiv surfaces.
  1100. upwd_check_surf.resize(MAX_DEPTH);
  1101. upwd_equiv_surf.resize(MAX_DEPTH);
  1102. dnwd_check_surf.resize(MAX_DEPTH);
  1103. dnwd_equiv_surf.resize(MAX_DEPTH);
  1104. for(size_t depth=0;depth<MAX_DEPTH;depth++){
  1105. Real_t c[3]={0.0,0.0,0.0};
  1106. upwd_check_surf[depth].ReInit((6*(m-1)*(m-1)+2)*COORD_DIM, buff_ptr, false); buff_ptr+=(6*(m-1)*(m-1)+2)*COORD_DIM;
  1107. upwd_equiv_surf[depth].ReInit((6*(m-1)*(m-1)+2)*COORD_DIM, buff_ptr, false); buff_ptr+=(6*(m-1)*(m-1)+2)*COORD_DIM;
  1108. dnwd_check_surf[depth].ReInit((6*(m-1)*(m-1)+2)*COORD_DIM, buff_ptr, false); buff_ptr+=(6*(m-1)*(m-1)+2)*COORD_DIM;
  1109. dnwd_equiv_surf[depth].ReInit((6*(m-1)*(m-1)+2)*COORD_DIM, buff_ptr, false); buff_ptr+=(6*(m-1)*(m-1)+2)*COORD_DIM;
  1110. upwd_check_surf[depth]=u_check_surf(m,c,depth);
  1111. upwd_equiv_surf[depth]=u_equiv_surf(m,c,depth);
  1112. dnwd_check_surf[depth]=d_check_surf(m,c,depth);
  1113. dnwd_equiv_surf[depth]=d_equiv_surf(m,c,depth);
  1114. }
  1115. }
  1116. buff[indx].AllocDevice(true);
  1117. }
  1118. }
  1119. template <class FMMNode>
  1120. void FMM_Pts<FMMNode>::SetupPrecomp(SetupData<Real_t>& setup_data, bool device){
  1121. if(setup_data.precomp_data==NULL) return;
  1122. Profile::Tic("SetupPrecomp",&this->comm,true,25);
  1123. { // Build precomp_data
  1124. size_t precomp_offset=0;
  1125. int level=setup_data.level;
  1126. Matrix<char>& precomp_data=*setup_data.precomp_data;
  1127. std::vector<Mat_Type>& interac_type_lst=setup_data.interac_type;
  1128. for(size_t type_indx=0; type_indx<interac_type_lst.size(); type_indx++){
  1129. Mat_Type& interac_type=interac_type_lst[type_indx];
  1130. this->PrecompAll(interac_type, level); // Compute matrices.
  1131. precomp_offset=this->mat->CompactData(level, interac_type, precomp_data, precomp_offset);
  1132. }
  1133. }
  1134. Profile::Toc();
  1135. if(device){ // Host2Device
  1136. Profile::Tic("Host2Device",&this->comm,false,25);
  1137. setup_data.precomp_data->AllocDevice(true);
  1138. Profile::Toc();
  1139. }
  1140. }
  1141. template <class FMMNode>
  1142. void FMM_Pts<FMMNode>::SetupInterac(SetupData<Real_t>& setup_data, bool device){
  1143. int level=setup_data.level;
  1144. std::vector<Mat_Type>& interac_type_lst=setup_data.interac_type;
  1145. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  1146. std::vector<void*>& nodes_out=setup_data.nodes_out;
  1147. Matrix<Real_t>& input_data=*setup_data. input_data;
  1148. Matrix<Real_t>& output_data=*setup_data.output_data;
  1149. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector;
  1150. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector;
  1151. size_t n_in =nodes_in .size();
  1152. size_t n_out=nodes_out.size();
  1153. // Setup precomputed data.
  1154. SetupPrecomp(setup_data,device);
  1155. // Build interac_data
  1156. Profile::Tic("Interac-Data",&this->comm,true,25);
  1157. Matrix<char>& interac_data=setup_data.interac_data;
  1158. if(n_out>0 && n_in >0){ // Build precomp_data, interac_data
  1159. std::vector<size_t> interac_mat;
  1160. std::vector<size_t> interac_cnt;
  1161. std::vector<size_t> interac_blk;
  1162. std::vector<size_t> input_perm;
  1163. std::vector<size_t> output_perm;
  1164. std::vector<Real_t> scaling;
  1165. size_t dof=0, M_dim0=0, M_dim1=0;
  1166. size_t precomp_offset=0;
  1167. size_t buff_size=DEVICE_BUFFER_SIZE*1024l*1024l;
  1168. for(size_t type_indx=0; type_indx<interac_type_lst.size(); type_indx++){
  1169. Mat_Type& interac_type=interac_type_lst[type_indx];
  1170. size_t mat_cnt=this->interac_list.ListCount(interac_type);
  1171. Vector<size_t> precomp_data_offset;
  1172. { // Load precomp_data for interac_type.
  1173. Matrix<char>& precomp_data=*setup_data.precomp_data;
  1174. char* indx_ptr=precomp_data[0]+precomp_offset;
  1175. size_t total_size=((size_t*)indx_ptr)[0]; indx_ptr+=sizeof(size_t);
  1176. /*size_t mat_cnt_ =((size_t*)indx_ptr)[0];*/ indx_ptr+=sizeof(size_t);
  1177. precomp_data_offset.ReInit((1+2+2)*mat_cnt, (size_t*)indx_ptr, false);
  1178. precomp_offset+=total_size;
  1179. }
  1180. Matrix<FMMNode*> src_interac_list(n_in ,mat_cnt); src_interac_list.SetZero();
  1181. Matrix<FMMNode*> trg_interac_list(n_out,mat_cnt); trg_interac_list.SetZero();
  1182. { // Build trg_interac_list
  1183. for(size_t i=0;i<n_out;i++){
  1184. if(!((FMMNode*)nodes_out[i])->IsGhost() && (level==-1 || ((FMMNode*)nodes_out[i])->Depth()==level)){
  1185. std::vector<FMMNode*>& lst=((FMMNode*)nodes_out[i])->interac_list[interac_type];
  1186. mem::memcopy(&trg_interac_list[i][0], &lst[0], lst.size()*sizeof(FMMNode*));
  1187. assert(lst.size()==mat_cnt);
  1188. }
  1189. }
  1190. }
  1191. { // Build src_interac_list
  1192. for(size_t i=0;i<n_in ;i++) ((FMMNode*)nodes_in [i])->node_id=i;
  1193. for(size_t i=0;i<n_out;i++){
  1194. for(size_t j=0;j<mat_cnt;j++)
  1195. if(trg_interac_list[i][j]!=NULL){
  1196. src_interac_list[trg_interac_list[i][j]->node_id][j]=(FMMNode*)nodes_out[i];
  1197. }
  1198. }
  1199. }
  1200. Matrix<size_t> interac_dsp(n_out,mat_cnt);
  1201. std::vector<size_t> interac_blk_dsp(1,0);
  1202. { // Determine dof, M_dim0, M_dim1
  1203. dof=1;
  1204. Matrix<Real_t>& M0 = this->interac_list.ClassMat(level, interac_type_lst[0], 0);
  1205. M_dim0=M0.Dim(0); M_dim1=M0.Dim(1);
  1206. }
  1207. { // Determine interaction blocks which fit in memory.
  1208. size_t vec_size=(M_dim0+M_dim1)*sizeof(Real_t)*dof;
  1209. for(size_t j=0;j<mat_cnt;j++){// Determine minimum buff_size
  1210. size_t vec_cnt=0;
  1211. for(size_t i=0;i<n_out;i++){
  1212. if(trg_interac_list[i][j]!=NULL) vec_cnt++;
  1213. }
  1214. if(buff_size<vec_cnt*vec_size)
  1215. buff_size=vec_cnt*vec_size;
  1216. }
  1217. size_t interac_dsp_=0;
  1218. for(size_t j=0;j<mat_cnt;j++){
  1219. for(size_t i=0;i<n_out;i++){
  1220. interac_dsp[i][j]=interac_dsp_;
  1221. if(trg_interac_list[i][j]!=NULL) interac_dsp_++;
  1222. }
  1223. if(interac_dsp_*vec_size>buff_size){
  1224. interac_blk.push_back(j-interac_blk_dsp.back());
  1225. interac_blk_dsp.push_back(j);
  1226. size_t offset=interac_dsp[0][j];
  1227. for(size_t i=0;i<n_out;i++) interac_dsp[i][j]-=offset;
  1228. interac_dsp_-=offset;
  1229. assert(interac_dsp_*vec_size<=buff_size); // Problem too big for buff_size.
  1230. }
  1231. interac_mat.push_back(precomp_data_offset[5*this->interac_list.InteracClass(interac_type,j)+0]);
  1232. interac_cnt.push_back(interac_dsp_-interac_dsp[0][j]);
  1233. }
  1234. interac_blk.push_back(mat_cnt-interac_blk_dsp.back());
  1235. interac_blk_dsp.push_back(mat_cnt);
  1236. }
  1237. { // Determine input_perm.
  1238. size_t vec_size=M_dim0*dof;
  1239. for(size_t i=0;i<n_out;i++) ((FMMNode*)nodes_out[i])->node_id=i;
  1240. for(size_t k=1;k<interac_blk_dsp.size();k++){
  1241. for(size_t i=0;i<n_in ;i++){
  1242. for(size_t j=interac_blk_dsp[k-1];j<interac_blk_dsp[k];j++){
  1243. FMMNode_t* trg_node=src_interac_list[i][j];
  1244. if(trg_node!=NULL){
  1245. input_perm .push_back(precomp_data_offset[5*j+1]); // prem
  1246. input_perm .push_back(precomp_data_offset[5*j+2]); // scal
  1247. input_perm .push_back(interac_dsp[trg_node->node_id][j]*vec_size*sizeof(Real_t)); // trg_ptr
  1248. input_perm .push_back((size_t)(& input_vector[i][0][0]- input_data[0])); // src_ptr
  1249. assert(input_vector[i]->Dim()==vec_size);
  1250. }
  1251. }
  1252. }
  1253. }
  1254. }
  1255. { // Determine scaling and output_perm
  1256. size_t vec_size=M_dim1*dof;
  1257. for(size_t k=1;k<interac_blk_dsp.size();k++){
  1258. for(size_t i=0;i<n_out;i++){
  1259. Real_t scaling_=0.0;
  1260. if(!this->Homogen()) scaling_=1.0;
  1261. else if(interac_type==S2U_Type) scaling_=pow(0.5, COORD_DIM *((FMMNode*)nodes_out[i])->Depth());
  1262. else if(interac_type==U2U_Type) scaling_=1.0;
  1263. else if(interac_type==D2D_Type) scaling_=1.0;
  1264. else if(interac_type==D2T_Type) scaling_=pow(0.5, -setup_data.kernel->poten_scale *((FMMNode*)nodes_out[i])->Depth());
  1265. else if(interac_type== U0_Type) scaling_=pow(0.5,(COORD_DIM-setup_data.kernel->poten_scale)*((FMMNode*)nodes_out[i])->Depth());
  1266. else if(interac_type== U1_Type) scaling_=pow(0.5,(COORD_DIM-setup_data.kernel->poten_scale)*((FMMNode*)nodes_out[i])->Depth());
  1267. else if(interac_type== U2_Type) scaling_=pow(0.5,(COORD_DIM-setup_data.kernel->poten_scale)*((FMMNode*)nodes_out[i])->Depth());
  1268. else if(interac_type== W_Type) scaling_=pow(0.5, -setup_data.kernel->poten_scale *((FMMNode*)nodes_out[i])->Depth());
  1269. else if(interac_type== X_Type) scaling_=pow(0.5, COORD_DIM *((FMMNode*)nodes_out[i])->Depth());
  1270. for(size_t j=interac_blk_dsp[k-1];j<interac_blk_dsp[k];j++){
  1271. if(trg_interac_list[i][j]!=NULL){
  1272. scaling.push_back(scaling_); // scaling
  1273. output_perm.push_back(precomp_data_offset[5*j+3]); // prem
  1274. output_perm.push_back(precomp_data_offset[5*j+4]); // scal
  1275. output_perm.push_back(interac_dsp[ i ][j]*vec_size*sizeof(Real_t)); // src_ptr
  1276. output_perm.push_back((size_t)(&output_vector[i][0][0]-output_data[0])); // trg_ptr
  1277. assert(output_vector[i]->Dim()==vec_size);
  1278. }
  1279. }
  1280. }
  1281. }
  1282. }
  1283. }
  1284. if(this->dev_buffer.Dim()<buff_size) this->dev_buffer.Resize(buff_size);
  1285. if(this->cpu_buffer.Dim()<buff_size) this->cpu_buffer.Resize(buff_size);
  1286. { // Set interac_data.
  1287. size_t data_size=sizeof(size_t)*4;
  1288. data_size+=sizeof(size_t)+interac_blk.size()*sizeof(size_t);
  1289. data_size+=sizeof(size_t)+interac_cnt.size()*sizeof(size_t);
  1290. data_size+=sizeof(size_t)+interac_mat.size()*sizeof(size_t);
  1291. data_size+=sizeof(size_t)+ input_perm.size()*sizeof(size_t);
  1292. data_size+=sizeof(size_t)+output_perm.size()*sizeof(size_t);
  1293. data_size+=sizeof(size_t)+scaling.size()*sizeof(Real_t);
  1294. if(interac_data.Dim(0)*interac_data.Dim(1)<sizeof(size_t)){
  1295. data_size+=sizeof(size_t);
  1296. interac_data.Resize(1,data_size);
  1297. ((size_t*)&interac_data[0][0])[0]=sizeof(size_t);
  1298. }else{
  1299. size_t pts_data_size=*((size_t*)&interac_data[0][0]);
  1300. assert(interac_data.Dim(0)*interac_data.Dim(1)>=pts_data_size);
  1301. data_size+=pts_data_size;
  1302. if(data_size>interac_data.Dim(0)*interac_data.Dim(1)){ //Resize and copy interac_data.
  1303. Matrix< char> pts_interac_data=interac_data;
  1304. interac_data.Resize(1,data_size);
  1305. mem::memcopy(&interac_data[0][0],&pts_interac_data[0][0],pts_data_size);
  1306. }
  1307. }
  1308. char* data_ptr=&interac_data[0][0];
  1309. data_ptr+=((size_t*)data_ptr)[0];
  1310. ((size_t*)data_ptr)[0]=data_size; data_ptr+=sizeof(size_t);
  1311. ((size_t*)data_ptr)[0]= M_dim0; data_ptr+=sizeof(size_t);
  1312. ((size_t*)data_ptr)[0]= M_dim1; data_ptr+=sizeof(size_t);
  1313. ((size_t*)data_ptr)[0]= dof; data_ptr+=sizeof(size_t);
  1314. ((size_t*)data_ptr)[0]=interac_blk.size(); data_ptr+=sizeof(size_t);
  1315. mem::memcopy(data_ptr, &interac_blk[0], interac_blk.size()*sizeof(size_t));
  1316. data_ptr+=interac_blk.size()*sizeof(size_t);
  1317. ((size_t*)data_ptr)[0]=interac_cnt.size(); data_ptr+=sizeof(size_t);
  1318. mem::memcopy(data_ptr, &interac_cnt[0], interac_cnt.size()*sizeof(size_t));
  1319. data_ptr+=interac_cnt.size()*sizeof(size_t);
  1320. ((size_t*)data_ptr)[0]=interac_mat.size(); data_ptr+=sizeof(size_t);
  1321. mem::memcopy(data_ptr, &interac_mat[0], interac_mat.size()*sizeof(size_t));
  1322. data_ptr+=interac_mat.size()*sizeof(size_t);
  1323. ((size_t*)data_ptr)[0]= input_perm.size(); data_ptr+=sizeof(size_t);
  1324. mem::memcopy(data_ptr, & input_perm[0], input_perm.size()*sizeof(size_t));
  1325. data_ptr+= input_perm.size()*sizeof(size_t);
  1326. ((size_t*)data_ptr)[0]=output_perm.size(); data_ptr+=sizeof(size_t);
  1327. mem::memcopy(data_ptr, &output_perm[0], output_perm.size()*sizeof(size_t));
  1328. data_ptr+=output_perm.size()*sizeof(size_t);
  1329. ((size_t*)data_ptr)[0]=scaling.size(); data_ptr+=sizeof(size_t);
  1330. mem::memcopy(data_ptr, &scaling[0], scaling.size()*sizeof(Real_t));
  1331. data_ptr+=scaling.size()*sizeof(Real_t);
  1332. }
  1333. }
  1334. Profile::Toc();
  1335. if(device){ // Host2Device
  1336. Profile::Tic("Host2Device",&this->comm,false,25);
  1337. setup_data.interac_data .AllocDevice(true);
  1338. Profile::Toc();
  1339. }
  1340. }
  1341. template <class FMMNode>
  1342. void FMM_Pts<FMMNode>::EvalList(SetupData<Real_t>& setup_data, bool device){
  1343. if(setup_data.interac_data.Dim(0)==0 || setup_data.interac_data.Dim(1)==0){
  1344. Profile::Tic("Host2Device",&this->comm,false,25);
  1345. Profile::Toc();
  1346. Profile::Tic("DeviceComp",&this->comm,false,20);
  1347. Profile::Toc();
  1348. return;
  1349. }
  1350. Profile::Tic("Host2Device",&this->comm,false,25);
  1351. typename Vector<char>::Device buff;
  1352. typename Matrix<char>::Device precomp_data;
  1353. typename Matrix<char>::Device interac_data;
  1354. typename Matrix<Real_t>::Device input_data;
  1355. typename Matrix<Real_t>::Device output_data;
  1356. if(device){
  1357. buff = this-> dev_buffer. AllocDevice(false);
  1358. precomp_data= setup_data.precomp_data->AllocDevice(false);
  1359. interac_data= setup_data.interac_data. AllocDevice(false);
  1360. input_data = setup_data. input_data->AllocDevice(false);
  1361. output_data = setup_data. output_data->AllocDevice(false);
  1362. }else{
  1363. buff = this-> cpu_buffer;
  1364. precomp_data=*setup_data.precomp_data;
  1365. interac_data= setup_data.interac_data;
  1366. input_data =*setup_data. input_data;
  1367. output_data =*setup_data. output_data;
  1368. }
  1369. Profile::Toc();
  1370. Profile::Tic("DeviceComp",&this->comm,false,20);
  1371. #ifdef __INTEL_OFFLOAD
  1372. int lock_idx=-1;
  1373. int wait_lock_idx=-1;
  1374. if(device) wait_lock_idx=MIC_Lock::curr_lock();
  1375. if(device) lock_idx=MIC_Lock::get_lock();
  1376. #pragma offload if(device) target(mic:0) signal(&MIC_Lock::lock_vec[device?lock_idx:0])
  1377. #endif
  1378. { // Offloaded computation.
  1379. // Set interac_data.
  1380. size_t data_size, M_dim0, M_dim1, dof;
  1381. Vector<size_t> interac_blk;
  1382. Vector<size_t> interac_cnt;
  1383. Vector<size_t> interac_mat;
  1384. Vector<size_t> input_perm;
  1385. Vector<size_t> output_perm;
  1386. Vector<Real_t> scaling;
  1387. { // Set interac_data.
  1388. char* data_ptr=&interac_data[0][0];
  1389. data_size=((size_t*)data_ptr)[0]; data_ptr+=data_size;
  1390. data_size=((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  1391. M_dim0 =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  1392. M_dim1 =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  1393. dof =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  1394. interac_blk.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  1395. data_ptr+=sizeof(size_t)+interac_blk.Dim()*sizeof(size_t);
  1396. interac_cnt.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  1397. data_ptr+=sizeof(size_t)+interac_cnt.Dim()*sizeof(size_t);
  1398. interac_mat.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  1399. data_ptr+=sizeof(size_t)+interac_mat.Dim()*sizeof(size_t);
  1400. input_perm .ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  1401. data_ptr+=sizeof(size_t)+ input_perm.Dim()*sizeof(size_t);
  1402. output_perm.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  1403. data_ptr+=sizeof(size_t)+output_perm.Dim()*sizeof(size_t);
  1404. scaling.ReInit(((size_t*)data_ptr)[0],(Real_t*)(data_ptr+sizeof(size_t)),false);
  1405. data_ptr+=sizeof(size_t)+scaling.Dim()*sizeof(Real_t);
  1406. }
  1407. #ifdef __INTEL_OFFLOAD
  1408. if(device) MIC_Lock::wait_lock(wait_lock_idx);
  1409. #endif
  1410. //Compute interaction from Chebyshev source density.
  1411. { // interactions
  1412. int omp_p=omp_get_max_threads();
  1413. size_t interac_indx=0;
  1414. size_t interac_blk_dsp=0;
  1415. for(size_t k=0;k<interac_blk.Dim();k++){
  1416. size_t vec_cnt=0;
  1417. for(size_t j=interac_blk_dsp;j<interac_blk_dsp+interac_blk[k];j++) vec_cnt+=interac_cnt[j];
  1418. char* buff_in =&buff[0];
  1419. char* buff_out=&buff[vec_cnt*dof*M_dim0*sizeof(Real_t)];
  1420. // Input permutation.
  1421. #pragma omp parallel for
  1422. for(int tid=0;tid<omp_p;tid++){
  1423. size_t a=( tid *vec_cnt)/omp_p;
  1424. size_t b=((tid+1)*vec_cnt)/omp_p;
  1425. for(size_t i=a;i<b;i++){
  1426. const PERM_INT_T* perm=(PERM_INT_T*)(precomp_data[0]+input_perm[(interac_indx+i)*4+0]);
  1427. const Real_t* scal=( Real_t*)(precomp_data[0]+input_perm[(interac_indx+i)*4+1]);
  1428. const Real_t* v_in =( Real_t*)( input_data[0]+input_perm[(interac_indx+i)*4+3]);
  1429. Real_t* v_out=( Real_t*)( buff_in +input_perm[(interac_indx+i)*4+2]);
  1430. // TODO: Fix for dof>1
  1431. #ifdef __MIC__
  1432. {
  1433. __m512d v8;
  1434. size_t j_start=(((uintptr_t)(v_out ) + (uintptr_t)(MEM_ALIGN-1)) & ~ (uintptr_t)(MEM_ALIGN-1))-((uintptr_t)v_out);
  1435. size_t j_end =(((uintptr_t)(v_out+M_dim0) ) & ~ (uintptr_t)(MEM_ALIGN-1))-((uintptr_t)v_out);
  1436. j_start/=sizeof(Real_t);
  1437. j_end /=sizeof(Real_t);
  1438. assert(((uintptr_t)(v_out))%sizeof(Real_t)==0);
  1439. assert(((uintptr_t)(v_out+j_start))%64==0);
  1440. assert(((uintptr_t)(v_out+j_end ))%64==0);
  1441. size_t j=0;
  1442. for(;j<j_start;j++ ){
  1443. v_out[j]=v_in[perm[j]]*scal[j];
  1444. }
  1445. for(;j<j_end ;j+=8){
  1446. v8=_mm512_setr_pd(
  1447. v_in[perm[j+0]]*scal[j+0],
  1448. v_in[perm[j+1]]*scal[j+1],
  1449. v_in[perm[j+2]]*scal[j+2],
  1450. v_in[perm[j+3]]*scal[j+3],
  1451. v_in[perm[j+4]]*scal[j+4],
  1452. v_in[perm[j+5]]*scal[j+5],
  1453. v_in[perm[j+6]]*scal[j+6],
  1454. v_in[perm[j+7]]*scal[j+7]);
  1455. _mm512_storenrngo_pd(v_out+j,v8);
  1456. }
  1457. for(;j<M_dim0 ;j++ ){
  1458. v_out[j]=v_in[perm[j]]*scal[j];
  1459. }
  1460. }
  1461. #else
  1462. for(size_t j=0;j<M_dim0;j++ ){
  1463. v_out[j]=v_in[perm[j]]*scal[j];
  1464. }
  1465. #endif
  1466. }
  1467. }
  1468. size_t vec_cnt0=0;
  1469. for(size_t j=interac_blk_dsp;j<interac_blk_dsp+interac_blk[k];){
  1470. size_t vec_cnt1=0;
  1471. size_t interac_mat0=interac_mat[j];
  1472. for(;j<interac_blk_dsp+interac_blk[k] && interac_mat[j]==interac_mat0;j++) vec_cnt1+=interac_cnt[j];
  1473. Matrix<Real_t> M(M_dim0, M_dim1, (Real_t*)(precomp_data[0]+interac_mat0), false);
  1474. #ifdef __MIC__
  1475. {
  1476. Matrix<Real_t> Ms(dof*vec_cnt1, M_dim0, (Real_t*)(buff_in +M_dim0*vec_cnt0*dof*sizeof(Real_t)), false);
  1477. Matrix<Real_t> Mt(dof*vec_cnt1, M_dim1, (Real_t*)(buff_out+M_dim1*vec_cnt0*dof*sizeof(Real_t)), false);
  1478. Matrix<Real_t>::DGEMM(Mt,Ms,M);
  1479. }
  1480. #else
  1481. #pragma omp parallel for
  1482. for(int tid=0;tid<omp_p;tid++){
  1483. size_t a=(dof*vec_cnt1*(tid ))/omp_p;
  1484. size_t b=(dof*vec_cnt1*(tid+1))/omp_p;
  1485. Matrix<Real_t> Ms(b-a, M_dim0, (Real_t*)(buff_in +M_dim0*vec_cnt0*dof*sizeof(Real_t))+M_dim0*a, false);
  1486. Matrix<Real_t> Mt(b-a, M_dim1, (Real_t*)(buff_out+M_dim1*vec_cnt0*dof*sizeof(Real_t))+M_dim1*a, false);
  1487. Matrix<Real_t>::DGEMM(Mt,Ms,M);
  1488. }
  1489. #endif
  1490. vec_cnt0+=vec_cnt1;
  1491. }
  1492. // Output permutation.
  1493. #pragma omp parallel for
  1494. for(int tid=0;tid<omp_p;tid++){
  1495. size_t a=( tid *vec_cnt)/omp_p;
  1496. size_t b=((tid+1)*vec_cnt)/omp_p;
  1497. if(tid> 0 && a<vec_cnt){ // Find 'a' independent of other threads.
  1498. size_t out_ptr=output_perm[(interac_indx+a)*4+3];
  1499. if(tid> 0) while(a<vec_cnt && out_ptr==output_perm[(interac_indx+a)*4+3]) a++;
  1500. }
  1501. if(tid<omp_p-1 && b<vec_cnt){ // Find 'b' independent of other threads.
  1502. size_t out_ptr=output_perm[(interac_indx+b)*4+3];
  1503. if(tid<omp_p-1) while(b<vec_cnt && out_ptr==output_perm[(interac_indx+b)*4+3]) b++;
  1504. }
  1505. for(size_t i=a;i<b;i++){ // Compute permutations.
  1506. Real_t scaling_factor=scaling[interac_indx+i];
  1507. const PERM_INT_T* perm=(PERM_INT_T*)(precomp_data[0]+output_perm[(interac_indx+i)*4+0]);
  1508. const Real_t* scal=( Real_t*)(precomp_data[0]+output_perm[(interac_indx+i)*4+1]);
  1509. const Real_t* v_in =( Real_t*)( buff_out +output_perm[(interac_indx+i)*4+2]);
  1510. Real_t* v_out=( Real_t*)( output_data[0]+output_perm[(interac_indx+i)*4+3]);
  1511. // TODO: Fix for dof>1
  1512. #ifdef __MIC__
  1513. {
  1514. __m512d v8;
  1515. __m512d v_old;
  1516. size_t j_start=(((uintptr_t)(v_out ) + (uintptr_t)(MEM_ALIGN-1)) & ~ (uintptr_t)(MEM_ALIGN-1))-((uintptr_t)v_out);
  1517. size_t j_end =(((uintptr_t)(v_out+M_dim1) ) & ~ (uintptr_t)(MEM_ALIGN-1))-((uintptr_t)v_out);
  1518. j_start/=sizeof(Real_t);
  1519. j_end /=sizeof(Real_t);
  1520. assert(((uintptr_t)(v_out))%sizeof(Real_t)==0);
  1521. assert(((uintptr_t)(v_out+j_start))%64==0);
  1522. assert(((uintptr_t)(v_out+j_end ))%64==0);
  1523. size_t j=0;
  1524. for(;j<j_start;j++ ){
  1525. v_out[j]+=v_in[perm[j]]*scal[j]*scaling_factor;
  1526. }
  1527. for(;j<j_end ;j+=8){
  1528. v_old=_mm512_load_pd(v_out+j);
  1529. v8=_mm512_setr_pd(
  1530. v_in[perm[j+0]]*scal[j+0]*scaling_factor,
  1531. v_in[perm[j+1]]*scal[j+1]*scaling_factor,
  1532. v_in[perm[j+2]]*scal[j+2]*scaling_factor,
  1533. v_in[perm[j+3]]*scal[j+3]*scaling_factor,
  1534. v_in[perm[j+4]]*scal[j+4]*scaling_factor,
  1535. v_in[perm[j+5]]*scal[j+5]*scaling_factor,
  1536. v_in[perm[j+6]]*scal[j+6]*scaling_factor,
  1537. v_in[perm[j+7]]*scal[j+7]*scaling_factor);
  1538. v_old=_mm512_add_pd(v_old, v8);
  1539. _mm512_storenrngo_pd(v_out+j,v_old);
  1540. }
  1541. for(;j<M_dim1 ;j++ ){
  1542. v_out[j]+=v_in[perm[j]]*scal[j]*scaling_factor;
  1543. }
  1544. }
  1545. #else
  1546. for(size_t j=0;j<M_dim1;j++ ){
  1547. v_out[j]+=v_in[perm[j]]*scal[j]*scaling_factor;
  1548. }
  1549. #endif
  1550. }
  1551. }
  1552. interac_indx+=vec_cnt;
  1553. interac_blk_dsp+=interac_blk[k];
  1554. }
  1555. }
  1556. #ifdef __INTEL_OFFLOAD
  1557. if(device) MIC_Lock::release_lock(lock_idx);
  1558. #endif
  1559. }
  1560. #ifndef __MIC_ASYNCH__
  1561. #ifdef __INTEL_OFFLOAD
  1562. #pragma offload if(device) target(mic:0)
  1563. {if(device) MIC_Lock::wait_lock(lock_idx);}
  1564. #endif
  1565. #endif
  1566. Profile::Toc();
  1567. }
  1568. template <class FMMNode>
  1569. void FMM_Pts<FMMNode>::Source2UpSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  1570. if(this->MultipoleOrder()==0) return;
  1571. { // Set setup_data
  1572. setup_data.level=level;
  1573. setup_data.kernel=&aux_kernel;
  1574. setup_data.interac_type.resize(1);
  1575. setup_data.interac_type[0]=S2U_Type;
  1576. setup_data. input_data=&buff[4];
  1577. setup_data.output_data=&buff[0];
  1578. setup_data. coord_data=&buff[6];
  1579. Vector<FMMNode_t*>& nodes_in =n_list[4];
  1580. Vector<FMMNode_t*>& nodes_out=n_list[0];
  1581. setup_data.nodes_in .clear();
  1582. setup_data.nodes_out.clear();
  1583. for(size_t i=0;i<nodes_in .Dim();i++) if(nodes_in [i]->Depth()==level || level==-1) setup_data.nodes_in .push_back(nodes_in [i]);
  1584. for(size_t i=0;i<nodes_out.Dim();i++) if(nodes_out[i]->Depth()==level || level==-1) setup_data.nodes_out.push_back(nodes_out[i]);
  1585. }
  1586. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  1587. std::vector<void*>& nodes_out=setup_data.nodes_out;
  1588. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  1589. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  1590. for(size_t i=0;i<nodes_in .size();i++){
  1591. input_vector .push_back(&((FMMNode*)nodes_in [i])->src_coord);
  1592. input_vector .push_back(&((FMMNode*)nodes_in [i])->src_value);
  1593. input_vector .push_back(&((FMMNode*)nodes_in [i])->surf_coord);
  1594. input_vector .push_back(&((FMMNode*)nodes_in [i])->surf_value);
  1595. }
  1596. for(size_t i=0;i<nodes_out.size();i++){
  1597. output_vector.push_back(&upwd_check_surf[((FMMNode*)nodes_out[i])->Depth()]);
  1598. output_vector.push_back(&((FMMData*)((FMMNode*)nodes_out[i])->FMMData())->upward_equiv);
  1599. }
  1600. //Upward check to upward equivalent matrix.
  1601. Matrix<Real_t>& M_uc2ue = this->mat->Mat(level, UC2UE_Type, 0);
  1602. this->SetupInteracPts(setup_data, false, true, &M_uc2ue,device);
  1603. { // Resize device buffer
  1604. size_t n=setup_data.output_data->Dim(0)*setup_data.output_data->Dim(1)*sizeof(Real_t);
  1605. if(this->dev_buffer.Dim()<n) this->dev_buffer.Resize(n);
  1606. }
  1607. }
  1608. template <class FMMNode>
  1609. void FMM_Pts<FMMNode>::Source2Up(SetupData<Real_t>& setup_data, bool device){
  1610. //Add Source2Up contribution.
  1611. this->EvalListPts(setup_data, device);
  1612. }
  1613. template <class FMMNode>
  1614. void FMM_Pts<FMMNode>::Up2UpSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  1615. if(this->MultipoleOrder()==0) return;
  1616. { // Set setup_data
  1617. setup_data.level=level;
  1618. setup_data.kernel=&aux_kernel;
  1619. setup_data.interac_type.resize(1);
  1620. setup_data.interac_type[0]=U2U_Type;
  1621. setup_data. input_data=&buff[0];
  1622. setup_data.output_data=&buff[0];
  1623. Vector<FMMNode_t*>& nodes_in =n_list[0];
  1624. Vector<FMMNode_t*>& nodes_out=n_list[0];
  1625. setup_data.nodes_in .clear();
  1626. setup_data.nodes_out.clear();
  1627. for(size_t i=0;i<nodes_in .Dim();i++) if(nodes_in [i]->Depth()==level+1) setup_data.nodes_in .push_back(nodes_in [i]);
  1628. for(size_t i=0;i<nodes_out.Dim();i++) if(nodes_out[i]->Depth()==level ) setup_data.nodes_out.push_back(nodes_out[i]);
  1629. }
  1630. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  1631. std::vector<void*>& nodes_out=setup_data.nodes_out;
  1632. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  1633. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  1634. for(size_t i=0;i<nodes_in .size();i++) input_vector.push_back(&((FMMData*)((FMMNode*)nodes_in [i])->FMMData())->upward_equiv);
  1635. for(size_t i=0;i<nodes_out.size();i++) output_vector.push_back(&((FMMData*)((FMMNode*)nodes_out[i])->FMMData())->upward_equiv);
  1636. SetupInterac(setup_data,device);
  1637. }
  1638. template <class FMMNode>
  1639. void FMM_Pts<FMMNode>::Up2Up (SetupData<Real_t>& setup_data, bool device){
  1640. //Add Up2Up contribution.
  1641. EvalList(setup_data, device);
  1642. }
  1643. template <class FMMNode>
  1644. void FMM_Pts<FMMNode>::PeriodicBC(FMMNode* node){
  1645. if(this->MultipoleOrder()==0) return;
  1646. Matrix<Real_t>& M = Precomp(0, BC_Type, 0);
  1647. assert(node->FMMData()->upward_equiv.Dim()>0);
  1648. int dof=1;
  1649. Vector<Real_t>& upward_equiv=node->FMMData()->upward_equiv;
  1650. Vector<Real_t>& dnward_equiv=node->FMMData()->dnward_equiv;
  1651. assert(upward_equiv.Dim()==M.Dim(0)*dof);
  1652. assert(dnward_equiv.Dim()==M.Dim(1)*dof);
  1653. Matrix<Real_t> d_equiv(dof,M.Dim(0),&dnward_equiv[0],false);
  1654. Matrix<Real_t> u_equiv(dof,M.Dim(1),&upward_equiv[0],false);
  1655. Matrix<Real_t>::DGEMM(d_equiv,u_equiv,M);
  1656. }
  1657. template <class FMMNode>
  1658. void FMM_Pts<FMMNode>::FFT_UpEquiv(size_t dof, size_t m, size_t ker_dim0, Vector<size_t>& fft_vec,
  1659. Vector<Real_t>& input_data, Vector<Real_t>& output_data, Vector<Real_t>& buffer_){
  1660. size_t n1=m*2;
  1661. size_t n2=n1*n1;
  1662. size_t n3=n1*n2;
  1663. size_t n3_=n2*(n1/2+1);
  1664. size_t chld_cnt=1UL<<COORD_DIM;
  1665. size_t fftsize_in =2*n3_*chld_cnt*ker_dim0*dof;
  1666. int omp_p=omp_get_max_threads();
  1667. //Load permutation map.
  1668. size_t n=6*(m-1)*(m-1)+2;
  1669. static Vector<size_t> map;
  1670. { // Build map to reorder upward_equiv
  1671. size_t n_old=map.Dim();
  1672. if(n_old!=n){
  1673. Real_t c[3]={0,0,0};
  1674. Vector<Real_t> surf=surface(m, c, (Real_t)(m-1), 0);
  1675. map.Resize(surf.Dim()/COORD_DIM);
  1676. for(size_t i=0;i<map.Dim();i++)
  1677. map[i]=((size_t)(m-1-surf[i*3]+0.5))+((size_t)(m-1-surf[i*3+1]+0.5))*n1+((size_t)(m-1-surf[i*3+2]+0.5))*n2;
  1678. }
  1679. }
  1680. { // Build FFTW plan.
  1681. if(!vlist_fft_flag){
  1682. int nnn[3]={(int)n1,(int)n1,(int)n1};
  1683. void *fftw_in, *fftw_out;
  1684. fftw_in = mem::aligned_malloc<Real_t>( n3 *ker_dim0*chld_cnt);
  1685. fftw_out = mem::aligned_malloc<Real_t>(2*n3_*ker_dim0*chld_cnt);
  1686. vlist_fftplan = FFTW_t<Real_t>::fft_plan_many_dft_r2c(COORD_DIM,nnn,ker_dim0*chld_cnt,
  1687. (Real_t*)fftw_in, NULL, 1, n3, (typename FFTW_t<Real_t>::cplx*)(fftw_out),NULL, 1, n3_, FFTW_ESTIMATE);
  1688. mem::aligned_free<Real_t>((Real_t*)fftw_in );
  1689. mem::aligned_free<Real_t>((Real_t*)fftw_out);
  1690. vlist_fft_flag=true;
  1691. }
  1692. }
  1693. { // Offload section
  1694. size_t n_in = fft_vec.Dim();
  1695. #pragma omp parallel for
  1696. for(int pid=0; pid<omp_p; pid++){
  1697. size_t node_start=(n_in*(pid ))/omp_p;
  1698. size_t node_end =(n_in*(pid+1))/omp_p;
  1699. Vector<Real_t> buffer(fftsize_in, &buffer_[fftsize_in*pid], false);
  1700. for(size_t node_idx=node_start; node_idx<node_end; node_idx++){
  1701. Vector<Real_t*> upward_equiv(chld_cnt);
  1702. for(size_t i=0;i<chld_cnt;i++) upward_equiv[i]=&input_data[0] + fft_vec[node_idx] + n*ker_dim0*dof*i;
  1703. Vector<Real_t> upward_equiv_fft(fftsize_in, &output_data[fftsize_in *node_idx], false);
  1704. upward_equiv_fft.SetZero();
  1705. // Rearrange upward equivalent data.
  1706. for(size_t k=0;k<n;k++){
  1707. size_t idx=map[k];
  1708. for(int j1=0;j1<dof;j1++)
  1709. for(int j0=0;j0<(int)chld_cnt;j0++)
  1710. for(int i=0;i<ker_dim0;i++)
  1711. upward_equiv_fft[idx+(j0+(i+j1*ker_dim0)*chld_cnt)*n3]=upward_equiv[j0][ker_dim0*(n*j1+k)+i];
  1712. }
  1713. // Compute FFT.
  1714. for(int i=0;i<dof;i++)
  1715. FFTW_t<Real_t>::fft_execute_dft_r2c(vlist_fftplan, (Real_t*)&upward_equiv_fft[i* n3 *ker_dim0*chld_cnt],
  1716. (typename FFTW_t<Real_t>::cplx*)&buffer [i*2*n3_*ker_dim0*chld_cnt]);
  1717. //Compute flops.
  1718. #ifndef FFTW3_MKL
  1719. double add, mul, fma;
  1720. fftw_flops(vlist_fftplan, &add, &mul, &fma);
  1721. #ifndef __INTEL_OFFLOAD0
  1722. Profile::Add_FLOP((long long)(add+mul+2*fma));
  1723. #endif
  1724. #endif
  1725. for(int i=0;i<ker_dim0*dof;i++)
  1726. for(size_t j=0;j<n3_;j++)
  1727. for(size_t k=0;k<chld_cnt;k++){
  1728. upward_equiv_fft[2*(chld_cnt*(n3_*i+j)+k)+0]=buffer[2*(n3_*(chld_cnt*i+k)+j)+0];
  1729. upward_equiv_fft[2*(chld_cnt*(n3_*i+j)+k)+1]=buffer[2*(n3_*(chld_cnt*i+k)+j)+1];
  1730. }
  1731. }
  1732. }
  1733. }
  1734. }
  1735. template <class FMMNode>
  1736. void FMM_Pts<FMMNode>::FFT_Check2Equiv(size_t dof, size_t m, size_t ker_dim1, Vector<size_t>& ifft_vec,
  1737. Vector<Real_t>& input_data, Vector<Real_t>& output_data, Vector<Real_t>& buffer_, Matrix<Real_t>& M){
  1738. size_t n1=m*2;
  1739. size_t n2=n1*n1;
  1740. size_t n3=n1*n2;
  1741. size_t n3_=n2*(n1/2+1);
  1742. size_t chld_cnt=1UL<<COORD_DIM;
  1743. size_t fftsize_out=2*n3_*dof*ker_dim1*chld_cnt;
  1744. int omp_p=omp_get_max_threads();
  1745. //Load permutation map.
  1746. size_t n=6*(m-1)*(m-1)+2;
  1747. static Vector<size_t> map;
  1748. { // Build map to reorder dnward_check
  1749. size_t n_old=map.Dim();
  1750. if(n_old!=n){
  1751. Real_t c[3]={0,0,0};
  1752. Vector<Real_t> surf=surface(m, c, (Real_t)(m-1), 0);
  1753. map.Resize(surf.Dim()/COORD_DIM);
  1754. for(size_t i=0;i<map.Dim();i++)
  1755. map[i]=((size_t)(m*2-0.5-surf[i*3]))+((size_t)(m*2-0.5-surf[i*3+1]))*n1+((size_t)(m*2-0.5-surf[i*3+2]))*n2;
  1756. //map;//.AllocDevice(true);
  1757. }
  1758. }
  1759. { // Build FFTW plan.
  1760. if(!vlist_ifft_flag){
  1761. //Build FFTW plan.
  1762. int nnn[3]={(int)n1,(int)n1,(int)n1};
  1763. void *fftw_in, *fftw_out;
  1764. fftw_in = fftw_malloc(2*n3_*ker_dim1*sizeof(Real_t)*chld_cnt);
  1765. fftw_out = fftw_malloc( n3 *ker_dim1*sizeof(Real_t)*chld_cnt);
  1766. vlist_ifftplan = FFTW_t<Real_t>::fft_plan_many_dft_c2r(COORD_DIM,nnn,ker_dim1*chld_cnt,
  1767. (typename FFTW_t<Real_t>::cplx*)fftw_in, NULL, 1, n3_, (Real_t*)(fftw_out),NULL, 1, n3, FFTW_ESTIMATE);
  1768. fftw_free(fftw_in);
  1769. fftw_free(fftw_out);
  1770. vlist_ifft_flag=true;
  1771. }
  1772. }
  1773. { // Offload section
  1774. size_t n_out=ifft_vec.Dim();
  1775. #pragma omp parallel for
  1776. for(int pid=0; pid<omp_p; pid++){
  1777. size_t node_start=(n_out*(pid ))/omp_p;
  1778. size_t node_end =(n_out*(pid+1))/omp_p;
  1779. Vector<Real_t> buffer(fftsize_out, &buffer_[fftsize_out*pid], false);
  1780. for(size_t node_idx=node_start; node_idx<node_end; node_idx++){
  1781. Vector<Real_t> dnward_check_fft(fftsize_out, &input_data[fftsize_out*node_idx], false);
  1782. //De-interleave data.
  1783. for(int i=0;i<ker_dim1*dof;i++)
  1784. for(size_t j=0;j<n3_;j++)
  1785. for(size_t k=0;k<chld_cnt;k++){
  1786. buffer[2*(n3_*(ker_dim1*dof*k+i)+j)+0]=dnward_check_fft[2*(chld_cnt*(n3_*i+j)+k)+0];
  1787. buffer[2*(n3_*(ker_dim1*dof*k+i)+j)+1]=dnward_check_fft[2*(chld_cnt*(n3_*i+j)+k)+1];
  1788. }
  1789. // Compute FFT.
  1790. for(int i=0;i<dof;i++)
  1791. FFTW_t<Real_t>::fft_execute_dft_c2r(vlist_ifftplan, (typename FFTW_t<Real_t>::cplx*)&buffer [i*2*n3_*ker_dim1*chld_cnt],
  1792. (Real_t*)&dnward_check_fft[i* n3 *ker_dim1*chld_cnt]);
  1793. //Compute flops.
  1794. #ifndef FFTW3_MKL
  1795. double add, mul, fma;
  1796. fftw_flops(vlist_ifftplan, &add, &mul, &fma);
  1797. #ifndef __INTEL_OFFLOAD0
  1798. Profile::Add_FLOP((long long)(add+mul+2*fma));
  1799. #endif
  1800. #endif
  1801. // Rearrange downward check data.
  1802. for(size_t k=0;k<n;k++){
  1803. size_t idx=map[k];
  1804. for(int j1=0;j1<dof;j1++)
  1805. for(int j0=0;j0<(int)chld_cnt;j0++)
  1806. for(int i=0;i<ker_dim1;i++)
  1807. buffer[ker_dim1*(n*(dof*j0+j1)+k)+i]=dnward_check_fft[idx+(j1+(i+j0*ker_dim1)*dof)*n3];
  1808. }
  1809. // Compute check to equiv.
  1810. for(size_t j=0;j<chld_cnt;j++){
  1811. Matrix<Real_t> d_check(dof,M.Dim(0),&buffer[n*ker_dim1*dof*j],false);
  1812. Matrix<Real_t> d_equiv(dof,M.Dim(1),&output_data[0] + ifft_vec[node_idx] + M.Dim(1)*dof*j,false);
  1813. Matrix<Real_t>::DGEMM(d_equiv,d_check,M,1.0);
  1814. }
  1815. }
  1816. }
  1817. }
  1818. }
  1819. template <class Real_t>
  1820. void VListHadamard(size_t dof, size_t M_dim, size_t ker_dim0, size_t ker_dim1, Vector<size_t>& interac_dsp,
  1821. Vector<size_t>& interac_vec, Vector<Real_t*>& precomp_mat, Vector<Real_t>& fft_in, Vector<Real_t>& fft_out){
  1822. size_t chld_cnt=1UL<<COORD_DIM;
  1823. size_t fftsize_in =M_dim*ker_dim0*chld_cnt*2;
  1824. size_t fftsize_out=M_dim*ker_dim1*chld_cnt*2;
  1825. Real_t* zero_vec0=mem::aligned_malloc<Real_t>(fftsize_in );
  1826. Real_t* zero_vec1=mem::aligned_malloc<Real_t>(fftsize_out);
  1827. size_t n_out=fft_out.Dim()/fftsize_out;
  1828. // Set buff_out to zero.
  1829. #pragma omp parallel for
  1830. for(size_t k=0;k<n_out;k++){
  1831. Vector<Real_t> dnward_check_fft(fftsize_out, &fft_out[k*fftsize_out], false);
  1832. dnward_check_fft.SetZero();
  1833. }
  1834. // Build list of interaction pairs (in, out vectors).
  1835. size_t mat_cnt=precomp_mat.Dim();
  1836. size_t blk1_cnt=interac_dsp.Dim()/mat_cnt;
  1837. Real_t** IN_ =new Real_t*[2*V_BLK_SIZE*blk1_cnt*mat_cnt];
  1838. Real_t** OUT_=new Real_t*[2*V_BLK_SIZE*blk1_cnt*mat_cnt];
  1839. #pragma omp parallel for
  1840. for(size_t interac_blk1=0; interac_blk1<blk1_cnt*mat_cnt; interac_blk1++){
  1841. size_t interac_dsp0 = (interac_blk1==0?0:interac_dsp[interac_blk1-1]);
  1842. size_t interac_dsp1 = interac_dsp[interac_blk1 ] ;
  1843. size_t interac_cnt = interac_dsp1-interac_dsp0;
  1844. for(size_t j=0;j<interac_cnt;j++){
  1845. IN_ [2*V_BLK_SIZE*interac_blk1 +j]=&fft_in [interac_vec[(interac_dsp0+j)*2+0]];
  1846. OUT_[2*V_BLK_SIZE*interac_blk1 +j]=&fft_out[interac_vec[(interac_dsp0+j)*2+1]];
  1847. }
  1848. IN_ [2*V_BLK_SIZE*interac_blk1 +interac_cnt]=zero_vec0;
  1849. OUT_[2*V_BLK_SIZE*interac_blk1 +interac_cnt]=zero_vec1;
  1850. }
  1851. int omp_p=omp_get_max_threads();
  1852. #pragma omp parallel for
  1853. for(int pid=0; pid<omp_p; pid++){
  1854. size_t a=( pid *M_dim)/omp_p;
  1855. size_t b=((pid+1)*M_dim)/omp_p;
  1856. for(size_t blk1=0; blk1<blk1_cnt; blk1++)
  1857. for(size_t k=a; k< b; k++)
  1858. for(size_t mat_indx=0; mat_indx< mat_cnt;mat_indx++){
  1859. size_t interac_blk1 = blk1*mat_cnt+mat_indx;
  1860. size_t interac_dsp0 = (interac_blk1==0?0:interac_dsp[interac_blk1-1]);
  1861. size_t interac_dsp1 = interac_dsp[interac_blk1 ] ;
  1862. size_t interac_cnt = interac_dsp1-interac_dsp0;
  1863. Real_t** IN = IN_ + 2*V_BLK_SIZE*interac_blk1;
  1864. Real_t** OUT= OUT_+ 2*V_BLK_SIZE*interac_blk1;
  1865. Real_t* M = precomp_mat[mat_indx] + k*chld_cnt*chld_cnt*2;
  1866. #ifdef __SSE__
  1867. if (mat_indx +1 < mat_cnt){ // Prefetch
  1868. _mm_prefetch(((char *)(precomp_mat[mat_indx+1] + k*chld_cnt*chld_cnt*2)), _MM_HINT_T0);
  1869. _mm_prefetch(((char *)(precomp_mat[mat_indx+1] + k*chld_cnt*chld_cnt*2) + 64), _MM_HINT_T0);
  1870. }
  1871. #endif
  1872. for(int in_dim=0;in_dim<ker_dim0;in_dim++)
  1873. for(int ot_dim=0;ot_dim<ker_dim1;ot_dim++){
  1874. for(size_t j=0;j<interac_cnt;j+=2){
  1875. Real_t* M_ = M;
  1876. Real_t* IN0 = IN [j+0] + (in_dim*M_dim+k)*chld_cnt*2;
  1877. Real_t* IN1 = IN [j+1] + (in_dim*M_dim+k)*chld_cnt*2;
  1878. Real_t* OUT0 = OUT[j+0] + (ot_dim*M_dim+k)*chld_cnt*2;
  1879. Real_t* OUT1 = OUT[j+1] + (ot_dim*M_dim+k)*chld_cnt*2;
  1880. #ifdef __SSE__
  1881. if (j+2 < interac_cnt) { // Prefetch
  1882. _mm_prefetch(((char *)(IN[j+2] + (in_dim*M_dim+k)*chld_cnt*2)), _MM_HINT_T0);
  1883. _mm_prefetch(((char *)(IN[j+2] + (in_dim*M_dim+k)*chld_cnt*2) + 64), _MM_HINT_T0);
  1884. _mm_prefetch(((char *)(IN[j+3] + (in_dim*M_dim+k)*chld_cnt*2)), _MM_HINT_T0);
  1885. _mm_prefetch(((char *)(IN[j+3] + (in_dim*M_dim+k)*chld_cnt*2) + 64), _MM_HINT_T0);
  1886. _mm_prefetch(((char *)(OUT[j+2] + (ot_dim*M_dim+k)*chld_cnt*2)), _MM_HINT_T0);
  1887. _mm_prefetch(((char *)(OUT[j+2] + (ot_dim*M_dim+k)*chld_cnt*2) + 64), _MM_HINT_T0);
  1888. _mm_prefetch(((char *)(OUT[j+3] + (ot_dim*M_dim+k)*chld_cnt*2)), _MM_HINT_T0);
  1889. _mm_prefetch(((char *)(OUT[j+3] + (ot_dim*M_dim+k)*chld_cnt*2) + 64), _MM_HINT_T0);
  1890. }
  1891. #endif
  1892. #ifdef __AVX__ //AVX code.
  1893. __m256d out00,out01,out10,out11;
  1894. __m256d out20,out21,out30,out31;
  1895. double* in0__ = IN0;
  1896. double* in1__ = IN1;
  1897. out00 = _mm256_load_pd(OUT0);
  1898. out01 = _mm256_load_pd(OUT1);
  1899. out10 = _mm256_load_pd(OUT0+4);
  1900. out11 = _mm256_load_pd(OUT1+4);
  1901. out20 = _mm256_load_pd(OUT0+8);
  1902. out21 = _mm256_load_pd(OUT1+8);
  1903. out30 = _mm256_load_pd(OUT0+12);
  1904. out31 = _mm256_load_pd(OUT1+12);
  1905. for(int i2=0;i2<8;i2+=2){
  1906. __m256d m00;
  1907. __m256d ot00;
  1908. __m256d mt0,mtt0;
  1909. __m256d in00,in00_r,in01,in01_r;
  1910. in00 = _mm256_broadcast_pd((const __m128d*)in0__);
  1911. in00_r = _mm256_permute_pd(in00,5);
  1912. in01 = _mm256_broadcast_pd((const __m128d*)in1__);
  1913. in01_r = _mm256_permute_pd(in01,5);
  1914. m00 = _mm256_load_pd(M_);
  1915. mt0 = _mm256_unpacklo_pd(m00,m00);
  1916. ot00 = _mm256_mul_pd(mt0,in00);
  1917. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1918. out00 = _mm256_add_pd(out00,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1919. ot00 = _mm256_mul_pd(mt0,in01);
  1920. out01 = _mm256_add_pd(out01,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1921. m00 = _mm256_load_pd(M_+4);
  1922. mt0 = _mm256_unpacklo_pd(m00,m00);
  1923. ot00 = _mm256_mul_pd(mt0,in00);
  1924. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1925. out10 = _mm256_add_pd(out10,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1926. ot00 = _mm256_mul_pd(mt0,in01);
  1927. out11 = _mm256_add_pd(out11,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1928. m00 = _mm256_load_pd(M_+8);
  1929. mt0 = _mm256_unpacklo_pd(m00,m00);
  1930. ot00 = _mm256_mul_pd(mt0,in00);
  1931. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1932. out20 = _mm256_add_pd(out20,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1933. ot00 = _mm256_mul_pd(mt0,in01);
  1934. out21 = _mm256_add_pd(out21,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1935. m00 = _mm256_load_pd(M_+12);
  1936. mt0 = _mm256_unpacklo_pd(m00,m00);
  1937. ot00 = _mm256_mul_pd(mt0,in00);
  1938. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1939. out30 = _mm256_add_pd(out30,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1940. ot00 = _mm256_mul_pd(mt0,in01);
  1941. out31 = _mm256_add_pd(out31,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1942. in00 = _mm256_broadcast_pd((const __m128d*) (in0__+2));
  1943. in00_r = _mm256_permute_pd(in00,5);
  1944. in01 = _mm256_broadcast_pd((const __m128d*) (in1__+2));
  1945. in01_r = _mm256_permute_pd(in01,5);
  1946. m00 = _mm256_load_pd(M_+16);
  1947. mt0 = _mm256_unpacklo_pd(m00,m00);
  1948. ot00 = _mm256_mul_pd(mt0,in00);
  1949. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1950. out00 = _mm256_add_pd(out00,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1951. ot00 = _mm256_mul_pd(mt0,in01);
  1952. out01 = _mm256_add_pd(out01,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1953. m00 = _mm256_load_pd(M_+20);
  1954. mt0 = _mm256_unpacklo_pd(m00,m00);
  1955. ot00 = _mm256_mul_pd(mt0,in00);
  1956. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1957. out10 = _mm256_add_pd(out10,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1958. ot00 = _mm256_mul_pd(mt0,in01);
  1959. out11 = _mm256_add_pd(out11,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1960. m00 = _mm256_load_pd(M_+24);
  1961. mt0 = _mm256_unpacklo_pd(m00,m00);
  1962. ot00 = _mm256_mul_pd(mt0,in00);
  1963. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1964. out20 = _mm256_add_pd(out20,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1965. ot00 = _mm256_mul_pd(mt0,in01);
  1966. out21 = _mm256_add_pd(out21,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1967. m00 = _mm256_load_pd(M_+28);
  1968. mt0 = _mm256_unpacklo_pd(m00,m00);
  1969. ot00 = _mm256_mul_pd(mt0,in00);
  1970. mtt0 = _mm256_unpackhi_pd(m00,m00);
  1971. out30 = _mm256_add_pd(out30,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in00_r)));
  1972. ot00 = _mm256_mul_pd(mt0,in01);
  1973. out31 = _mm256_add_pd(out31,_mm256_addsub_pd(ot00,_mm256_mul_pd(mtt0,in01_r)));
  1974. M_ += 32;
  1975. in0__ += 4;
  1976. in1__ += 4;
  1977. }
  1978. _mm256_store_pd(OUT0,out00);
  1979. _mm256_store_pd(OUT1,out01);
  1980. _mm256_store_pd(OUT0+4,out10);
  1981. _mm256_store_pd(OUT1+4,out11);
  1982. _mm256_store_pd(OUT0+8,out20);
  1983. _mm256_store_pd(OUT1+8,out21);
  1984. _mm256_store_pd(OUT0+12,out30);
  1985. _mm256_store_pd(OUT1+12,out31);
  1986. #elif defined __SSE3__ // SSE code.
  1987. __m128d out00, out01, out10, out11;
  1988. __m128d in00, in01, in10, in11;
  1989. __m128d m00, m01, m10, m11;
  1990. //#pragma unroll
  1991. for(int i1=0;i1<8;i1+=2){
  1992. double* IN0_=IN0;
  1993. double* IN1_=IN1;
  1994. out00 =_mm_load_pd (OUT0 );
  1995. out10 =_mm_load_pd (OUT0+2);
  1996. out01 =_mm_load_pd (OUT1 );
  1997. out11 =_mm_load_pd (OUT1+2);
  1998. //#pragma unroll
  1999. for(int i2=0;i2<8;i2+=2){
  2000. m00 =_mm_load1_pd (M_ );
  2001. m10 =_mm_load1_pd (M_+ 2);
  2002. m01 =_mm_load1_pd (M_+16);
  2003. m11 =_mm_load1_pd (M_+18);
  2004. in00 =_mm_load_pd (IN0_ );
  2005. in10 =_mm_load_pd (IN0_+2);
  2006. in01 =_mm_load_pd (IN1_ );
  2007. in11 =_mm_load_pd (IN1_+2);
  2008. out00 = _mm_add_pd (out00, _mm_mul_pd(m00 , in00 ));
  2009. out00 = _mm_add_pd (out00, _mm_mul_pd(m01 , in10 ));
  2010. out01 = _mm_add_pd (out01, _mm_mul_pd(m00 , in01 ));
  2011. out01 = _mm_add_pd (out01, _mm_mul_pd(m01 , in11 ));
  2012. out10 = _mm_add_pd (out10, _mm_mul_pd(m10 , in00 ));
  2013. out10 = _mm_add_pd (out10, _mm_mul_pd(m11 , in10 ));
  2014. out11 = _mm_add_pd (out11, _mm_mul_pd(m10 , in01 ));
  2015. out11 = _mm_add_pd (out11, _mm_mul_pd(m11 , in11 ));
  2016. m00 =_mm_load1_pd (M_+ 1);
  2017. m10 =_mm_load1_pd (M_+ 2+1);
  2018. m01 =_mm_load1_pd (M_+16+1);
  2019. m11 =_mm_load1_pd (M_+18+1);
  2020. in00 =_mm_shuffle_pd (in00,in00,_MM_SHUFFLE2(0,1));
  2021. in01 =_mm_shuffle_pd (in01,in01,_MM_SHUFFLE2(0,1));
  2022. in10 =_mm_shuffle_pd (in10,in10,_MM_SHUFFLE2(0,1));
  2023. in11 =_mm_shuffle_pd (in11,in11,_MM_SHUFFLE2(0,1));
  2024. out00 = _mm_addsub_pd(out00, _mm_mul_pd(m00, in00));
  2025. out00 = _mm_addsub_pd(out00, _mm_mul_pd(m01, in10));
  2026. out01 = _mm_addsub_pd(out01, _mm_mul_pd(m00, in01));
  2027. out01 = _mm_addsub_pd(out01, _mm_mul_pd(m01, in11));
  2028. out10 = _mm_addsub_pd(out10, _mm_mul_pd(m10, in00));
  2029. out10 = _mm_addsub_pd(out10, _mm_mul_pd(m11, in10));
  2030. out11 = _mm_addsub_pd(out11, _mm_mul_pd(m10, in01));
  2031. out11 = _mm_addsub_pd(out11, _mm_mul_pd(m11, in11));
  2032. M_+=32; // Jump to (column+2).
  2033. IN0_+=4;
  2034. IN1_+=4;
  2035. }
  2036. _mm_store_pd (OUT0 ,out00);
  2037. _mm_store_pd (OUT0+2,out10);
  2038. _mm_store_pd (OUT1 ,out01);
  2039. _mm_store_pd (OUT1+2,out11);
  2040. M_+=4-64*2; // Jump back to first column (row+2).
  2041. OUT0+=4;
  2042. OUT1+=4;
  2043. }
  2044. #else // Generic code.
  2045. Real_t out_reg000, out_reg001, out_reg010, out_reg011;
  2046. Real_t out_reg100, out_reg101, out_reg110, out_reg111;
  2047. Real_t in_reg000, in_reg001, in_reg010, in_reg011;
  2048. Real_t in_reg100, in_reg101, in_reg110, in_reg111;
  2049. Real_t m_reg000, m_reg001, m_reg010, m_reg011;
  2050. Real_t m_reg100, m_reg101, m_reg110, m_reg111;
  2051. //#pragma unroll
  2052. for(int i1=0;i1<8;i1+=2){
  2053. Real_t* IN0_=IN0;
  2054. Real_t* IN1_=IN1;
  2055. out_reg000=OUT0[ 0]; out_reg001=OUT0[ 1];
  2056. out_reg010=OUT0[ 2]; out_reg011=OUT0[ 3];
  2057. out_reg100=OUT1[ 0]; out_reg101=OUT1[ 1];
  2058. out_reg110=OUT1[ 2]; out_reg111=OUT1[ 3];
  2059. //#pragma unroll
  2060. for(int i2=0;i2<8;i2+=2){
  2061. m_reg000=M_[ 0]; m_reg001=M_[ 1];
  2062. m_reg010=M_[ 2]; m_reg011=M_[ 3];
  2063. m_reg100=M_[16]; m_reg101=M_[17];
  2064. m_reg110=M_[18]; m_reg111=M_[19];
  2065. in_reg000=IN0_[0]; in_reg001=IN0_[1];
  2066. in_reg010=IN0_[2]; in_reg011=IN0_[3];
  2067. in_reg100=IN1_[0]; in_reg101=IN1_[1];
  2068. in_reg110=IN1_[2]; in_reg111=IN1_[3];
  2069. out_reg000 += m_reg000*in_reg000 - m_reg001*in_reg001;
  2070. out_reg001 += m_reg000*in_reg001 + m_reg001*in_reg000;
  2071. out_reg010 += m_reg010*in_reg000 - m_reg011*in_reg001;
  2072. out_reg011 += m_reg010*in_reg001 + m_reg011*in_reg000;
  2073. out_reg000 += m_reg100*in_reg010 - m_reg101*in_reg011;
  2074. out_reg001 += m_reg100*in_reg011 + m_reg101*in_reg010;
  2075. out_reg010 += m_reg110*in_reg010 - m_reg111*in_reg011;
  2076. out_reg011 += m_reg110*in_reg011 + m_reg111*in_reg010;
  2077. out_reg100 += m_reg000*in_reg100 - m_reg001*in_reg101;
  2078. out_reg101 += m_reg000*in_reg101 + m_reg001*in_reg100;
  2079. out_reg110 += m_reg010*in_reg100 - m_reg011*in_reg101;
  2080. out_reg111 += m_reg010*in_reg101 + m_reg011*in_reg100;
  2081. out_reg100 += m_reg100*in_reg110 - m_reg101*in_reg111;
  2082. out_reg101 += m_reg100*in_reg111 + m_reg101*in_reg110;
  2083. out_reg110 += m_reg110*in_reg110 - m_reg111*in_reg111;
  2084. out_reg111 += m_reg110*in_reg111 + m_reg111*in_reg110;
  2085. M_+=32; // Jump to (column+2).
  2086. IN0_+=4;
  2087. IN1_+=4;
  2088. }
  2089. OUT0[ 0]=out_reg000; OUT0[ 1]=out_reg001;
  2090. OUT0[ 2]=out_reg010; OUT0[ 3]=out_reg011;
  2091. OUT1[ 0]=out_reg100; OUT1[ 1]=out_reg101;
  2092. OUT1[ 2]=out_reg110; OUT1[ 3]=out_reg111;
  2093. M_+=4-64*2; // Jump back to first column (row+2).
  2094. OUT0+=4;
  2095. OUT1+=4;
  2096. }
  2097. #endif
  2098. }
  2099. M += M_dim*128;
  2100. }
  2101. }
  2102. }
  2103. // Compute flops.
  2104. {
  2105. Profile::Add_FLOP(8*8*8*(interac_vec.Dim()/2)*M_dim*ker_dim0*ker_dim1*dof);
  2106. }
  2107. // Free memory
  2108. delete[] IN_ ;
  2109. delete[] OUT_;
  2110. mem::aligned_free<Real_t>(zero_vec0);
  2111. mem::aligned_free<Real_t>(zero_vec1);
  2112. }
  2113. template <class FMMNode>
  2114. void FMM_Pts<FMMNode>::V_ListSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  2115. if(this->MultipoleOrder()==0) return;
  2116. if(level==0) return;
  2117. { // Set setup_data
  2118. setup_data.level=level;
  2119. setup_data.kernel=&aux_kernel;
  2120. setup_data.interac_type.resize(1);
  2121. setup_data.interac_type[0]=V1_Type;
  2122. setup_data. input_data=&buff[0];
  2123. setup_data.output_data=&buff[1];
  2124. Vector<FMMNode_t*>& nodes_in =n_list[2];
  2125. Vector<FMMNode_t*>& nodes_out=n_list[3];
  2126. setup_data.nodes_in .clear();
  2127. setup_data.nodes_out.clear();
  2128. for(size_t i=0;i<nodes_in .Dim();i++) if(nodes_in [i]->Depth()==level-1 || level==-1) setup_data.nodes_in .push_back(nodes_in [i]);
  2129. for(size_t i=0;i<nodes_out.Dim();i++) if(nodes_out[i]->Depth()==level-1 || level==-1) setup_data.nodes_out.push_back(nodes_out[i]);
  2130. }
  2131. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  2132. std::vector<void*>& nodes_out=setup_data.nodes_out;
  2133. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  2134. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  2135. for(size_t i=0;i<nodes_in .size();i++) input_vector.push_back(&((FMMData*)((FMMNode*)((FMMNode*)nodes_in [i])->Child(0))->FMMData())->upward_equiv);
  2136. for(size_t i=0;i<nodes_out.size();i++) output_vector.push_back(&((FMMData*)((FMMNode*)((FMMNode*)nodes_out[i])->Child(0))->FMMData())->dnward_equiv);
  2137. /////////////////////////////////////////////////////////////////////////////
  2138. size_t n_in =nodes_in .size();
  2139. size_t n_out=nodes_out.size();
  2140. // Setup precomputed data.
  2141. SetupPrecomp(setup_data,device);
  2142. // Build interac_data
  2143. Profile::Tic("Interac-Data",&this->comm,true,25);
  2144. Matrix<char>& interac_data=setup_data.interac_data;
  2145. if(n_out>0 && n_in >0){ // Build precomp_data, interac_data
  2146. size_t precomp_offset=0;
  2147. Mat_Type& interac_type=setup_data.interac_type[0];
  2148. size_t mat_cnt=this->interac_list.ListCount(interac_type);
  2149. Vector<size_t> precomp_data_offset;
  2150. std::vector<size_t> interac_mat;
  2151. { // Load precomp_data for interac_type.
  2152. Matrix<char>& precomp_data=*setup_data.precomp_data;
  2153. char* indx_ptr=precomp_data[0]+precomp_offset;
  2154. size_t total_size=((size_t*)indx_ptr)[0]; indx_ptr+=sizeof(size_t);
  2155. /*size_t mat_cnt_ =((size_t*)indx_ptr)[0];*/ indx_ptr+=sizeof(size_t);
  2156. precomp_data_offset.ReInit((1+2+2)*mat_cnt, (size_t*)indx_ptr, false);
  2157. precomp_offset+=total_size;
  2158. for(size_t mat_id=0;mat_id<mat_cnt;mat_id++){
  2159. Matrix<Real_t>& M0 = this->mat->Mat(level, interac_type, mat_id);
  2160. assert(M0.Dim(0)>0 && M0.Dim(1)>0); UNUSED(M0);
  2161. interac_mat.push_back(precomp_data_offset[5*mat_id]);
  2162. }
  2163. }
  2164. size_t dof;
  2165. size_t m=MultipoleOrder();
  2166. size_t ker_dim0=setup_data.kernel->ker_dim[0];
  2167. size_t ker_dim1=setup_data.kernel->ker_dim[1];
  2168. size_t fftsize;
  2169. {
  2170. size_t n1=m*2;
  2171. size_t n2=n1*n1;
  2172. size_t n3_=n2*(n1/2+1);
  2173. size_t chld_cnt=1UL<<COORD_DIM;
  2174. fftsize=2*n3_*chld_cnt;
  2175. dof=1;
  2176. }
  2177. int omp_p=omp_get_max_threads();
  2178. size_t buff_size=DEVICE_BUFFER_SIZE*1024l*1024l;
  2179. size_t n_blk0=2*fftsize*dof*(ker_dim0*n_in +ker_dim1*n_out)*sizeof(Real_t)/buff_size;
  2180. if(n_blk0==0) n_blk0=1;
  2181. std::vector<std::vector<size_t> > fft_vec(n_blk0);
  2182. std::vector<std::vector<size_t> > ifft_vec(n_blk0);
  2183. std::vector<std::vector<size_t> > interac_vec(n_blk0);
  2184. std::vector<std::vector<size_t> > interac_dsp(n_blk0);
  2185. {
  2186. Matrix<Real_t>& input_data=*setup_data. input_data;
  2187. Matrix<Real_t>& output_data=*setup_data.output_data;
  2188. std::vector<std::vector<FMMNode*> > nodes_blk_in (n_blk0);
  2189. std::vector<std::vector<FMMNode*> > nodes_blk_out(n_blk0);
  2190. for(size_t i=0;i<n_in;i++) ((FMMNode*)nodes_in[i])->node_id=i;
  2191. for(size_t blk0=0;blk0<n_blk0;blk0++){
  2192. size_t blk0_start=(n_out* blk0 )/n_blk0;
  2193. size_t blk0_end =(n_out*(blk0+1))/n_blk0;
  2194. std::vector<FMMNode*>& nodes_in_ =nodes_blk_in [blk0];
  2195. std::vector<FMMNode*>& nodes_out_=nodes_blk_out[blk0];
  2196. { // Build node list for blk0.
  2197. std::set<void*> nodes_in;
  2198. for(size_t i=blk0_start;i<blk0_end;i++){
  2199. nodes_out_.push_back((FMMNode*)nodes_out[i]);
  2200. std::vector<FMMNode*>& lst=((FMMNode*)nodes_out[i])->interac_list[interac_type];
  2201. for(size_t k=0;k<mat_cnt;k++) if(lst[k]!=NULL) nodes_in.insert(lst[k]);
  2202. }
  2203. for(std::set<void*>::iterator node=nodes_in.begin(); node != nodes_in.end(); node++){
  2204. nodes_in_.push_back((FMMNode*)*node);
  2205. }
  2206. size_t input_dim=nodes_in_ .size()*ker_dim0*dof*fftsize;
  2207. size_t output_dim=nodes_out_.size()*ker_dim1*dof*fftsize;
  2208. size_t buffer_dim=(ker_dim0+ker_dim1)*dof*fftsize*omp_p;
  2209. if(buff_size<(input_dim + output_dim + buffer_dim)*sizeof(Real_t))
  2210. buff_size=(input_dim + output_dim + buffer_dim)*sizeof(Real_t);
  2211. }
  2212. { // Set fft vectors.
  2213. for(size_t i=0;i<nodes_in_ .size();i++) fft_vec[blk0].push_back((size_t)(& input_vector[nodes_in_[i]->node_id][0][0]- input_data[0]));
  2214. for(size_t i=0;i<nodes_out_.size();i++)ifft_vec[blk0].push_back((size_t)(&output_vector[blk0_start + i ][0][0]-output_data[0]));
  2215. }
  2216. }
  2217. for(size_t blk0=0;blk0<n_blk0;blk0++){ // Hadamard interactions.
  2218. std::vector<FMMNode*>& nodes_in_ =nodes_blk_in [blk0];
  2219. std::vector<FMMNode*>& nodes_out_=nodes_blk_out[blk0];
  2220. for(size_t i=0;i<nodes_in_.size();i++) nodes_in_[i]->node_id=i;
  2221. { // Next blocking level.
  2222. size_t n_blk1=nodes_out_.size()*(ker_dim0+ker_dim1)/V_BLK_SIZE; //64 vectors should fit in L1.
  2223. if(n_blk1==0) n_blk1=1;
  2224. size_t interac_dsp_=0;
  2225. for(size_t blk1=0;blk1<n_blk1;blk1++){
  2226. size_t blk1_start=(nodes_out_.size()* blk1 )/n_blk1;
  2227. size_t blk1_end =(nodes_out_.size()*(blk1+1))/n_blk1;
  2228. for(size_t k=0;k<mat_cnt;k++){
  2229. for(size_t i=blk1_start;i<blk1_end;i++){
  2230. std::vector<FMMNode*>& lst=((FMMNode*)nodes_out_[i])->interac_list[interac_type];
  2231. if(lst[k]!=NULL){
  2232. interac_vec[blk0].push_back(lst[k]->node_id*fftsize*ker_dim0*dof);
  2233. interac_vec[blk0].push_back( i *fftsize*ker_dim1*dof);
  2234. interac_dsp_++;
  2235. }
  2236. }
  2237. interac_dsp[blk0].push_back(interac_dsp_);
  2238. }
  2239. }
  2240. }
  2241. }
  2242. }
  2243. { // Set interac_data.
  2244. size_t data_size=sizeof(size_t)*5; // m, dof, ker_dim0, ker_dim1, n_blk0
  2245. for(size_t blk0=0;blk0<n_blk0;blk0++){
  2246. data_size+=sizeof(size_t)+ fft_vec[blk0].size()*sizeof(size_t);
  2247. data_size+=sizeof(size_t)+ ifft_vec[blk0].size()*sizeof(size_t);
  2248. data_size+=sizeof(size_t)+interac_vec[blk0].size()*sizeof(size_t);
  2249. data_size+=sizeof(size_t)+interac_dsp[blk0].size()*sizeof(size_t);
  2250. }
  2251. data_size+=sizeof(size_t)+interac_mat.size()*sizeof(size_t);
  2252. if(data_size>interac_data.Dim(0)*interac_data.Dim(1))
  2253. interac_data.Resize(1,data_size);
  2254. char* data_ptr=&interac_data[0][0];
  2255. ((size_t*)data_ptr)[0]=buff_size; data_ptr+=sizeof(size_t);
  2256. ((size_t*)data_ptr)[0]= m; data_ptr+=sizeof(size_t);
  2257. ((size_t*)data_ptr)[0]= dof; data_ptr+=sizeof(size_t);
  2258. ((size_t*)data_ptr)[0]= ker_dim0; data_ptr+=sizeof(size_t);
  2259. ((size_t*)data_ptr)[0]= ker_dim1; data_ptr+=sizeof(size_t);
  2260. ((size_t*)data_ptr)[0]= n_blk0; data_ptr+=sizeof(size_t);
  2261. ((size_t*)data_ptr)[0]= interac_mat.size(); data_ptr+=sizeof(size_t);
  2262. mem::memcopy(data_ptr, &interac_mat[0], interac_mat.size()*sizeof(size_t));
  2263. data_ptr+=interac_mat.size()*sizeof(size_t);
  2264. for(size_t blk0=0;blk0<n_blk0;blk0++){
  2265. ((size_t*)data_ptr)[0]= fft_vec[blk0].size(); data_ptr+=sizeof(size_t);
  2266. mem::memcopy(data_ptr, & fft_vec[blk0][0], fft_vec[blk0].size()*sizeof(size_t));
  2267. data_ptr+= fft_vec[blk0].size()*sizeof(size_t);
  2268. ((size_t*)data_ptr)[0]=ifft_vec[blk0].size(); data_ptr+=sizeof(size_t);
  2269. mem::memcopy(data_ptr, &ifft_vec[blk0][0], ifft_vec[blk0].size()*sizeof(size_t));
  2270. data_ptr+=ifft_vec[blk0].size()*sizeof(size_t);
  2271. ((size_t*)data_ptr)[0]=interac_vec[blk0].size(); data_ptr+=sizeof(size_t);
  2272. mem::memcopy(data_ptr, &interac_vec[blk0][0], interac_vec[blk0].size()*sizeof(size_t));
  2273. data_ptr+=interac_vec[blk0].size()*sizeof(size_t);
  2274. ((size_t*)data_ptr)[0]=interac_dsp[blk0].size(); data_ptr+=sizeof(size_t);
  2275. mem::memcopy(data_ptr, &interac_dsp[blk0][0], interac_dsp[blk0].size()*sizeof(size_t));
  2276. data_ptr+=interac_dsp[blk0].size()*sizeof(size_t);
  2277. }
  2278. }
  2279. }
  2280. Profile::Toc();
  2281. Profile::Tic("Host2Device",&this->comm,false,25);
  2282. if(device){ // Host2Device
  2283. setup_data.interac_data. AllocDevice(true);
  2284. }
  2285. Profile::Toc();
  2286. }
  2287. template <class FMMNode>
  2288. void FMM_Pts<FMMNode>::V_List (SetupData<Real_t>& setup_data, bool device){
  2289. assert(!device); //Can not run on accelerator yet.
  2290. if(setup_data.interac_data.Dim(0)==0 || setup_data.interac_data.Dim(1)==0){
  2291. Profile::Tic("Host2Device",&this->comm,false,25);
  2292. Profile::Toc();
  2293. Profile::Tic("FFT",&comm,false,100);
  2294. Profile::Toc();
  2295. Profile::Tic("HadamardProduct",&comm,false,100);
  2296. Profile::Toc();
  2297. Profile::Tic("IFFT",&comm,false,100);
  2298. Profile::Toc();
  2299. return;
  2300. }
  2301. Profile::Tic("Host2Device",&this->comm,false,25);
  2302. int level=setup_data.level;
  2303. size_t buff_size=*((size_t*)&setup_data.interac_data[0][0]);
  2304. typename Matrix<Real_t>::Device M_d;
  2305. typename Vector<char>::Device buff;
  2306. typename Matrix<char>::Device precomp_data;
  2307. typename Matrix<char>::Device interac_data;
  2308. typename Matrix<Real_t>::Device input_data;
  2309. typename Matrix<Real_t>::Device output_data;
  2310. Matrix<Real_t>& M = this->mat->Mat(level, DC2DE_Type, 0);
  2311. if(device){
  2312. if(this->dev_buffer.Dim()<buff_size) this->dev_buffer.Resize(buff_size);
  2313. M_d = M. AllocDevice(false);
  2314. buff = this-> dev_buffer. AllocDevice(false);
  2315. precomp_data= setup_data.precomp_data->AllocDevice(false);
  2316. interac_data= setup_data.interac_data. AllocDevice(false);
  2317. input_data = setup_data. input_data->AllocDevice(false);
  2318. output_data = setup_data. output_data->AllocDevice(false);
  2319. }else{
  2320. if(this->cpu_buffer.Dim()<buff_size) this->cpu_buffer.Resize(buff_size);
  2321. M_d = M;
  2322. buff = this-> cpu_buffer;
  2323. precomp_data=*setup_data.precomp_data;
  2324. interac_data= setup_data.interac_data;
  2325. input_data =*setup_data. input_data;
  2326. output_data =*setup_data. output_data;
  2327. }
  2328. Profile::Toc();
  2329. { // Offloaded computation.
  2330. // Set interac_data.
  2331. size_t m, dof, ker_dim0, ker_dim1, n_blk0;
  2332. std::vector<Vector<size_t> > fft_vec;
  2333. std::vector<Vector<size_t> > ifft_vec;
  2334. std::vector<Vector<size_t> > interac_vec;
  2335. std::vector<Vector<size_t> > interac_dsp;
  2336. Vector<Real_t*> precomp_mat;
  2337. { // Set interac_data.
  2338. char* data_ptr=&interac_data[0][0];
  2339. buff_size=((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2340. m =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2341. dof =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2342. ker_dim0 =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2343. ker_dim1 =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2344. n_blk0 =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2345. fft_vec .resize(n_blk0);
  2346. ifft_vec.resize(n_blk0);
  2347. interac_vec.resize(n_blk0);
  2348. interac_dsp.resize(n_blk0);
  2349. Vector<size_t> interac_mat;
  2350. interac_mat.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2351. data_ptr+=sizeof(size_t)+interac_mat.Dim()*sizeof(size_t);
  2352. precomp_mat.Resize(interac_mat.Dim());
  2353. for(size_t i=0;i<interac_mat.Dim();i++){
  2354. precomp_mat[i]=(Real_t*)(precomp_data[0]+interac_mat[i]);
  2355. }
  2356. for(size_t blk0=0;blk0<n_blk0;blk0++){
  2357. fft_vec[blk0].ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2358. data_ptr+=sizeof(size_t)+fft_vec[blk0].Dim()*sizeof(size_t);
  2359. ifft_vec[blk0].ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2360. data_ptr+=sizeof(size_t)+ifft_vec[blk0].Dim()*sizeof(size_t);
  2361. interac_vec[blk0].ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2362. data_ptr+=sizeof(size_t)+interac_vec[blk0].Dim()*sizeof(size_t);
  2363. interac_dsp[blk0].ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2364. data_ptr+=sizeof(size_t)+interac_dsp[blk0].Dim()*sizeof(size_t);
  2365. }
  2366. }
  2367. int omp_p=omp_get_max_threads();
  2368. size_t M_dim, fftsize;
  2369. {
  2370. size_t n1=m*2;
  2371. size_t n2=n1*n1;
  2372. size_t n3_=n2*(n1/2+1);
  2373. size_t chld_cnt=1UL<<COORD_DIM;
  2374. fftsize=2*n3_*chld_cnt;
  2375. M_dim=n3_;
  2376. }
  2377. for(size_t blk0=0;blk0<n_blk0;blk0++){ // interactions
  2378. size_t n_in = fft_vec[blk0].Dim();
  2379. size_t n_out=ifft_vec[blk0].Dim();
  2380. size_t input_dim=n_in *ker_dim0*dof*fftsize;
  2381. size_t output_dim=n_out*ker_dim1*dof*fftsize;
  2382. size_t buffer_dim=(ker_dim0+ker_dim1)*dof*fftsize*omp_p;
  2383. Vector<Real_t> fft_in ( input_dim, (Real_t*)&buff[ 0 ],false);
  2384. Vector<Real_t> fft_out(output_dim, (Real_t*)&buff[ input_dim *sizeof(Real_t)],false);
  2385. Vector<Real_t> buffer(buffer_dim, (Real_t*)&buff[(input_dim+output_dim)*sizeof(Real_t)],false);
  2386. { // FFT
  2387. Profile::Tic("FFT",&comm,false,100);
  2388. Vector<Real_t> input_data_( input_data.dim[0]* input_data.dim[1], input_data[0], false);
  2389. FFT_UpEquiv(dof, m, ker_dim0, fft_vec[blk0], input_data_, fft_in, buffer);
  2390. Profile::Toc();
  2391. }
  2392. { // Hadamard
  2393. #ifdef PVFMM_HAVE_PAPI
  2394. #ifdef __VERBOSE__
  2395. std::cout << "Starting counters new\n";
  2396. if (PAPI_start(EventSet) != PAPI_OK) std::cout << "handle_error3" << std::endl;
  2397. #endif
  2398. #endif
  2399. Profile::Tic("HadamardProduct",&comm,false,100);
  2400. VListHadamard<Real_t>(dof, M_dim, ker_dim0, ker_dim1, interac_dsp[blk0], interac_vec[blk0], precomp_mat, fft_in, fft_out);
  2401. Profile::Toc();
  2402. #ifdef PVFMM_HAVE_PAPI
  2403. #ifdef __VERBOSE__
  2404. if (PAPI_stop(EventSet, values) != PAPI_OK) std::cout << "handle_error4" << std::endl;
  2405. std::cout << "Stopping counters\n";
  2406. #endif
  2407. #endif
  2408. }
  2409. { // IFFT
  2410. Profile::Tic("IFFT",&comm,false,100);
  2411. Matrix<Real_t> M(M_d.dim[0],M_d.dim[1],M_d[0],false);
  2412. Vector<Real_t> output_data_(output_data.dim[0]*output_data.dim[1], output_data[0], false);
  2413. FFT_Check2Equiv(dof, m, ker_dim1, ifft_vec[blk0], fft_out, output_data_, buffer, M);
  2414. Profile::Toc();
  2415. }
  2416. }
  2417. }
  2418. }
  2419. template <class FMMNode>
  2420. void FMM_Pts<FMMNode>::Down2DownSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  2421. if(this->MultipoleOrder()==0) return;
  2422. { // Set setup_data
  2423. setup_data.level=level;
  2424. setup_data.kernel=&aux_kernel;
  2425. setup_data.interac_type.resize(1);
  2426. setup_data.interac_type[0]=D2D_Type;
  2427. setup_data. input_data=&buff[1];
  2428. setup_data.output_data=&buff[1];
  2429. Vector<FMMNode_t*>& nodes_in =n_list[1];
  2430. Vector<FMMNode_t*>& nodes_out=n_list[1];
  2431. setup_data.nodes_in .clear();
  2432. setup_data.nodes_out.clear();
  2433. for(size_t i=0;i<nodes_in .Dim();i++) if(nodes_in [i]->Depth()==level-1) setup_data.nodes_in .push_back(nodes_in [i]);
  2434. for(size_t i=0;i<nodes_out.Dim();i++) if(nodes_out[i]->Depth()==level ) setup_data.nodes_out.push_back(nodes_out[i]);
  2435. }
  2436. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  2437. std::vector<void*>& nodes_out=setup_data.nodes_out;
  2438. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  2439. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  2440. for(size_t i=0;i<nodes_in .size();i++) input_vector.push_back(&((FMMData*)((FMMNode*)nodes_in [i])->FMMData())->dnward_equiv);
  2441. for(size_t i=0;i<nodes_out.size();i++) output_vector.push_back(&((FMMData*)((FMMNode*)nodes_out[i])->FMMData())->dnward_equiv);
  2442. SetupInterac(setup_data,device);
  2443. }
  2444. template <class FMMNode>
  2445. void FMM_Pts<FMMNode>::Down2Down (SetupData<Real_t>& setup_data, bool device){
  2446. //Add Down2Down contribution.
  2447. EvalList(setup_data, device);
  2448. }
  2449. template <class FMMNode>
  2450. void FMM_Pts<FMMNode>::SetupInteracPts(SetupData<Real_t>& setup_data, bool shift_src, bool shift_trg, Matrix<Real_t>* M, bool device){
  2451. int level=setup_data.level;
  2452. std::vector<Mat_Type>& interac_type_lst=setup_data.interac_type;
  2453. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  2454. std::vector<void*>& nodes_out=setup_data.nodes_out;
  2455. Matrix<Real_t>& output_data=*setup_data.output_data;
  2456. Matrix<Real_t>& input_data=*setup_data. input_data;
  2457. Matrix<Real_t>& coord_data=*setup_data. coord_data;
  2458. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector;
  2459. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector;
  2460. size_t n_in =nodes_in .size();
  2461. size_t n_out=nodes_out.size();
  2462. //setup_data.precomp_data=NULL;
  2463. // Build interac_data
  2464. Profile::Tic("Interac-Data",&this->comm,true,25);
  2465. Matrix<char>& interac_data=setup_data.interac_data;
  2466. if(n_out>0 && n_in >0){
  2467. size_t ker_dim0=setup_data.kernel->ker_dim[0];
  2468. size_t ker_dim1=setup_data.kernel->ker_dim[1];
  2469. size_t dof=1;
  2470. for(size_t i=0;i<n_in ;i++) ((FMMNode*)nodes_in [i])->node_id=i;
  2471. std::vector<size_t> trg_interac_cnt(n_out,0);
  2472. std::vector<size_t> trg_coord(n_out);
  2473. std::vector<size_t> trg_value(n_out);
  2474. std::vector<size_t> trg_cnt(n_out);
  2475. std::vector<Real_t> scaling(n_out,0);
  2476. { // Set trg data
  2477. Mat_Type& interac_type=interac_type_lst[0];
  2478. for(size_t i=0;i<n_out;i++){
  2479. if(!((FMMNode*)nodes_out[i])->IsGhost() && (level==-1 || ((FMMNode*)nodes_out[i])->Depth()==level)){
  2480. trg_cnt [i]=output_vector[i*2+0]->Dim()/COORD_DIM;
  2481. trg_coord[i]=(size_t)(&output_vector[i*2+0][0][0]- coord_data[0]);
  2482. trg_value[i]=(size_t)(&output_vector[i*2+1][0][0]-output_data[0]);
  2483. if(!this->Homogen()) scaling[i]=1.0;
  2484. else if(interac_type==S2U_Type) scaling[i]=pow(0.5, setup_data.kernel->poten_scale *((FMMNode*)nodes_out[i])->Depth());
  2485. else if(interac_type== X_Type) scaling[i]=pow(0.5, setup_data.kernel->poten_scale *((FMMNode*)nodes_out[i])->Depth());
  2486. }
  2487. }
  2488. }
  2489. std::vector<std::vector<size_t> > src_cnt(n_out);
  2490. std::vector<std::vector<size_t> > src_coord(n_out);
  2491. std::vector<std::vector<size_t> > src_value(n_out);
  2492. std::vector<std::vector<Real_t> > shift_coord(n_out);
  2493. for(size_t type_indx=0; type_indx<interac_type_lst.size(); type_indx++){
  2494. Mat_Type& interac_type=interac_type_lst[type_indx];
  2495. size_t mat_cnt=this->interac_list.ListCount(interac_type);
  2496. for(size_t i=0;i<n_out;i++){ // For each target node.
  2497. if(!((FMMNode*)nodes_out[i])->IsGhost() && (level==-1 || ((FMMNode*)nodes_out[i])->Depth()==level)){
  2498. std::vector<FMMNode*>& lst=((FMMNode*)nodes_out[i])->interac_list[interac_type];
  2499. for(size_t mat_indx=0;mat_indx<mat_cnt;mat_indx++) if(lst[mat_indx]!=NULL){ // For each direction.
  2500. size_t j=lst[mat_indx]->node_id;
  2501. if(input_vector[j*4+0]->Dim()>0 || input_vector[j*4+2]->Dim()>0){
  2502. trg_interac_cnt[i]++;
  2503. { // Determine shift for periodic boundary condition
  2504. Real_t* sc=lst[mat_indx]->Coord();
  2505. Real_t* tc=((FMMNode*)nodes_out[i])->Coord();
  2506. int* rel_coord=this->interac_list.RelativeCoord(interac_type, mat_indx);
  2507. shift_coord[i].push_back((tc[0]>sc[0] && rel_coord[0]>0? 1.0:
  2508. (tc[0]<sc[0] && rel_coord[0]<0?-1.0:0.0)) +
  2509. (shift_src?sc[0]:0) - (shift_trg?tc[0]:0) );
  2510. shift_coord[i].push_back((tc[1]>sc[1] && rel_coord[1]>0? 1.0:
  2511. (tc[1]<sc[1] && rel_coord[1]<0?-1.0:0.0)) +
  2512. (shift_src?sc[1]:0) - (shift_trg?tc[1]:0) );
  2513. shift_coord[i].push_back((tc[2]>sc[2] && rel_coord[2]>0? 1.0:
  2514. (tc[2]<sc[2] && rel_coord[2]<0?-1.0:0.0)) +
  2515. (shift_src?sc[2]:0) - (shift_trg?tc[2]:0) );
  2516. }
  2517. { // Set src data
  2518. if(input_vector[j*4+0]!=NULL){
  2519. src_cnt [i].push_back(input_vector[j*4+0]->Dim()/COORD_DIM);
  2520. src_coord[i].push_back((size_t)(& input_vector[j*4+0][0][0]- coord_data[0]));
  2521. src_value[i].push_back((size_t)(& input_vector[j*4+1][0][0]- input_data[0]));
  2522. }else{
  2523. src_cnt [i].push_back(0);
  2524. src_coord[i].push_back(0);
  2525. src_value[i].push_back(0);
  2526. }
  2527. if(input_vector[j*4+2]!=NULL){
  2528. src_cnt [i].push_back(input_vector[j*4+2]->Dim()/COORD_DIM);
  2529. src_coord[i].push_back((size_t)(& input_vector[j*4+2][0][0]- coord_data[0]));
  2530. src_value[i].push_back((size_t)(& input_vector[j*4+3][0][0]- input_data[0]));
  2531. }else{
  2532. src_cnt [i].push_back(0);
  2533. src_coord[i].push_back(0);
  2534. src_value[i].push_back(0);
  2535. }
  2536. }
  2537. }
  2538. }
  2539. }
  2540. }
  2541. }
  2542. { // Set interac_data.
  2543. size_t data_size=sizeof(size_t)*4;
  2544. data_size+=sizeof(size_t)+trg_interac_cnt.size()*sizeof(size_t);
  2545. data_size+=sizeof(size_t)+trg_coord.size()*sizeof(size_t);
  2546. data_size+=sizeof(size_t)+trg_value.size()*sizeof(size_t);
  2547. data_size+=sizeof(size_t)+trg_cnt .size()*sizeof(size_t);
  2548. data_size+=sizeof(size_t)+scaling .size()*sizeof(Real_t);
  2549. data_size+=sizeof(size_t)*2+(M!=NULL?M->Dim(0)*M->Dim(1)*sizeof(Real_t):0);
  2550. for(size_t i=0;i<n_out;i++){
  2551. data_size+=sizeof(size_t)+src_cnt [i].size()*sizeof(size_t);
  2552. data_size+=sizeof(size_t)+src_coord[i].size()*sizeof(size_t);
  2553. data_size+=sizeof(size_t)+src_value[i].size()*sizeof(size_t);
  2554. data_size+=sizeof(size_t)+shift_coord[i].size()*sizeof(Real_t);
  2555. }
  2556. if(data_size>interac_data.Dim(0)*interac_data.Dim(1))
  2557. interac_data.Resize(1,data_size);
  2558. char* data_ptr=&interac_data[0][0];
  2559. ((size_t*)data_ptr)[0]=data_size; data_ptr+=sizeof(size_t);
  2560. ((size_t*)data_ptr)[0]= ker_dim0; data_ptr+=sizeof(size_t);
  2561. ((size_t*)data_ptr)[0]= ker_dim1; data_ptr+=sizeof(size_t);
  2562. ((size_t*)data_ptr)[0]= dof; data_ptr+=sizeof(size_t);
  2563. ((size_t*)data_ptr)[0]=trg_interac_cnt.size(); data_ptr+=sizeof(size_t);
  2564. mem::memcopy(data_ptr, &trg_interac_cnt[0], trg_interac_cnt.size()*sizeof(size_t));
  2565. data_ptr+=trg_interac_cnt.size()*sizeof(size_t);
  2566. ((size_t*)data_ptr)[0]=trg_coord.size(); data_ptr+=sizeof(size_t);
  2567. mem::memcopy(data_ptr, &trg_coord[0], trg_coord.size()*sizeof(size_t));
  2568. data_ptr+=trg_coord.size()*sizeof(size_t);
  2569. ((size_t*)data_ptr)[0]=trg_value.size(); data_ptr+=sizeof(size_t);
  2570. mem::memcopy(data_ptr, &trg_value[0], trg_value.size()*sizeof(size_t));
  2571. data_ptr+=trg_value.size()*sizeof(size_t);
  2572. ((size_t*)data_ptr)[0]=trg_cnt.size(); data_ptr+=sizeof(size_t);
  2573. mem::memcopy(data_ptr, &trg_cnt[0], trg_cnt.size()*sizeof(size_t));
  2574. data_ptr+=trg_cnt.size()*sizeof(size_t);
  2575. ((size_t*)data_ptr)[0]=scaling.size(); data_ptr+=sizeof(size_t);
  2576. mem::memcopy(data_ptr, &scaling[0], scaling.size()*sizeof(Real_t));
  2577. data_ptr+=scaling.size()*sizeof(Real_t);
  2578. if(M!=NULL){
  2579. ((size_t*)data_ptr)[0]=M->Dim(0); data_ptr+=sizeof(size_t);
  2580. ((size_t*)data_ptr)[0]=M->Dim(1); data_ptr+=sizeof(size_t);
  2581. mem::memcopy(data_ptr, M[0][0], M->Dim(0)*M->Dim(1)*sizeof(Real_t));
  2582. data_ptr+=M->Dim(0)*M->Dim(1)*sizeof(Real_t);
  2583. }else{
  2584. ((size_t*)data_ptr)[0]=0; data_ptr+=sizeof(size_t);
  2585. ((size_t*)data_ptr)[0]=0; data_ptr+=sizeof(size_t);
  2586. }
  2587. for(size_t i=0;i<n_out;i++){
  2588. ((size_t*)data_ptr)[0]=src_cnt[i].size(); data_ptr+=sizeof(size_t);
  2589. mem::memcopy(data_ptr, &src_cnt[i][0], src_cnt[i].size()*sizeof(size_t));
  2590. data_ptr+=src_cnt[i].size()*sizeof(size_t);
  2591. ((size_t*)data_ptr)[0]=src_coord[i].size(); data_ptr+=sizeof(size_t);
  2592. mem::memcopy(data_ptr, &src_coord[i][0], src_coord[i].size()*sizeof(size_t));
  2593. data_ptr+=src_coord[i].size()*sizeof(size_t);
  2594. ((size_t*)data_ptr)[0]=src_value[i].size(); data_ptr+=sizeof(size_t);
  2595. mem::memcopy(data_ptr, &src_value[i][0], src_value[i].size()*sizeof(size_t));
  2596. data_ptr+=src_value[i].size()*sizeof(size_t);
  2597. ((size_t*)data_ptr)[0]=shift_coord[i].size(); data_ptr+=sizeof(size_t);
  2598. mem::memcopy(data_ptr, &shift_coord[i][0], shift_coord[i].size()*sizeof(Real_t));
  2599. data_ptr+=shift_coord[i].size()*sizeof(Real_t);
  2600. }
  2601. }
  2602. size_t buff_size=DEVICE_BUFFER_SIZE*1024l*1024l;
  2603. if(this->dev_buffer.Dim()<buff_size) this->dev_buffer.Resize(buff_size);
  2604. if(this->cpu_buffer.Dim()<buff_size) this->cpu_buffer.Resize(buff_size);
  2605. }
  2606. Profile::Toc();
  2607. Profile::Tic("Host2Device",&this->comm,false,25);
  2608. if(device){ // Host2Device
  2609. setup_data.interac_data .AllocDevice(true);
  2610. }
  2611. Profile::Toc();
  2612. }
  2613. template <class FMMNode>
  2614. void FMM_Pts<FMMNode>::EvalListPts(SetupData<Real_t>& setup_data, bool device){
  2615. if(setup_data.interac_data.Dim(0)==0 || setup_data.interac_data.Dim(1)==0){
  2616. Profile::Tic("Host2Device",&this->comm,false,25);
  2617. Profile::Toc();
  2618. Profile::Tic("DeviceComp",&this->comm,false,20);
  2619. Profile::Toc();
  2620. return;
  2621. }
  2622. Profile::Tic("Host2Device",&this->comm,false,25);
  2623. typename Vector<char>::Device buff;
  2624. //typename Matrix<char>::Device precomp_data;
  2625. typename Matrix<char>::Device interac_data;
  2626. typename Matrix<Real_t>::Device coord_data;
  2627. typename Matrix<Real_t>::Device input_data;
  2628. typename Matrix<Real_t>::Device output_data;
  2629. if(device){
  2630. buff = this-> dev_buffer. AllocDevice(false);
  2631. interac_data= setup_data.interac_data. AllocDevice(false);
  2632. //if(setup_data.precomp_data!=NULL) precomp_data= setup_data.precomp_data->AllocDevice(false);
  2633. if(setup_data. coord_data!=NULL) coord_data = setup_data. coord_data->AllocDevice(false);
  2634. if(setup_data. input_data!=NULL) input_data = setup_data. input_data->AllocDevice(false);
  2635. if(setup_data. output_data!=NULL) output_data = setup_data. output_data->AllocDevice(false);
  2636. }else{
  2637. buff = this-> cpu_buffer;
  2638. interac_data= setup_data.interac_data;
  2639. //if(setup_data.precomp_data!=NULL) precomp_data=*setup_data.precomp_data;
  2640. if(setup_data. coord_data!=NULL) coord_data =*setup_data. coord_data;
  2641. if(setup_data. input_data!=NULL) input_data =*setup_data. input_data;
  2642. if(setup_data. output_data!=NULL) output_data =*setup_data. output_data;
  2643. }
  2644. Profile::Toc();
  2645. size_t ptr_single_layer_kernel=(size_t)setup_data.kernel->ker_poten;
  2646. size_t ptr_double_layer_kernel=(size_t)setup_data.kernel->dbl_layer_poten;
  2647. Profile::Tic("DeviceComp",&this->comm,false,20);
  2648. #ifdef __INTEL_OFFLOAD
  2649. int lock_idx=-1;
  2650. int wait_lock_idx=-1;
  2651. if(device) wait_lock_idx=MIC_Lock::curr_lock();
  2652. if(device) lock_idx=MIC_Lock::get_lock();
  2653. if(device) ptr_single_layer_kernel=setup_data.kernel->dev_ker_poten;
  2654. if(device) ptr_double_layer_kernel=setup_data.kernel->dev_dbl_layer_poten;
  2655. #pragma offload if(device) target(mic:0) signal(&MIC_Lock::lock_vec[device?lock_idx:0])
  2656. #endif
  2657. { // Offloaded computation.
  2658. // Set interac_data.
  2659. //size_t data_size;
  2660. //size_t ker_dim0;
  2661. size_t ker_dim1;
  2662. size_t dof, n_out;
  2663. Vector<size_t> trg_interac_cnt;
  2664. Vector<size_t> trg_coord;
  2665. Vector<size_t> trg_value;
  2666. Vector<size_t> trg_cnt;
  2667. Vector<Real_t> scaling;
  2668. Matrix<Real_t> M;
  2669. Vector< Vector<size_t> > src_cnt;
  2670. Vector< Vector<size_t> > src_coord;
  2671. Vector< Vector<size_t> > src_value;
  2672. Vector< Vector<Real_t> > shift_coord;
  2673. { // Set interac_data.
  2674. char* data_ptr=&interac_data[0][0];
  2675. /*data_size=((size_t*)data_ptr)[0];*/ data_ptr+=sizeof(size_t);
  2676. /*ker_dim0=((size_t*)data_ptr)[0];*/ data_ptr+=sizeof(size_t);
  2677. ker_dim1=((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2678. dof =((size_t*)data_ptr)[0]; data_ptr+=sizeof(size_t);
  2679. trg_interac_cnt.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2680. data_ptr+=sizeof(size_t)+trg_interac_cnt.Dim()*sizeof(size_t);
  2681. n_out=trg_interac_cnt.Dim();
  2682. trg_coord.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2683. data_ptr+=sizeof(size_t)+trg_coord.Dim()*sizeof(size_t);
  2684. trg_value.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2685. data_ptr+=sizeof(size_t)+trg_value.Dim()*sizeof(size_t);
  2686. trg_cnt.ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2687. data_ptr+=sizeof(size_t)+trg_cnt.Dim()*sizeof(size_t);
  2688. scaling.ReInit(((size_t*)data_ptr)[0],(Real_t*)(data_ptr+sizeof(size_t)),false);
  2689. data_ptr+=sizeof(size_t)+scaling.Dim()*sizeof(Real_t);
  2690. M.ReInit(((size_t*)data_ptr)[0],((size_t*)data_ptr)[1],(Real_t*)(data_ptr+2*sizeof(size_t)),false);
  2691. data_ptr+=sizeof(size_t)*2+M.Dim(0)*M.Dim(1)*sizeof(Real_t);
  2692. src_cnt.Resize(n_out);
  2693. src_coord.Resize(n_out);
  2694. src_value.Resize(n_out);
  2695. shift_coord.Resize(n_out);
  2696. for(size_t i=0;i<n_out;i++){
  2697. src_cnt[i].ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2698. data_ptr+=sizeof(size_t)+src_cnt[i].Dim()*sizeof(size_t);
  2699. src_coord[i].ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2700. data_ptr+=sizeof(size_t)+src_coord[i].Dim()*sizeof(size_t);
  2701. src_value[i].ReInit(((size_t*)data_ptr)[0],(size_t*)(data_ptr+sizeof(size_t)),false);
  2702. data_ptr+=sizeof(size_t)+src_value[i].Dim()*sizeof(size_t);
  2703. shift_coord[i].ReInit(((size_t*)data_ptr)[0],(Real_t*)(data_ptr+sizeof(size_t)),false);
  2704. data_ptr+=sizeof(size_t)+shift_coord[i].Dim()*sizeof(Real_t);
  2705. }
  2706. }
  2707. #ifdef __INTEL_OFFLOAD
  2708. if(device) MIC_Lock::wait_lock(wait_lock_idx);
  2709. #endif
  2710. //Compute interaction from point sources.
  2711. { // interactions
  2712. typename Kernel<Real_t>::Ker_t single_layer_kernel=(typename Kernel<Real_t>::Ker_t)ptr_single_layer_kernel;
  2713. typename Kernel<Real_t>::Ker_t double_layer_kernel=(typename Kernel<Real_t>::Ker_t)ptr_double_layer_kernel;
  2714. int omp_p=omp_get_max_threads();
  2715. Vector<Real_t*> thread_buff(omp_p);
  2716. size_t thread_buff_size=buff.dim/sizeof(Real_t)/omp_p;
  2717. for(int i=0;i<omp_p;i++) thread_buff[i]=(Real_t*)&buff[i*thread_buff_size*sizeof(Real_t)];
  2718. #pragma omp parallel for
  2719. for(size_t i=0;i<n_out;i++)
  2720. if(trg_interac_cnt[i]>0 && trg_cnt[i]>0){
  2721. int thread_id=omp_get_thread_num();
  2722. Real_t* s_coord=thread_buff[thread_id];
  2723. Real_t* t_value=output_data[0]+trg_value[i];
  2724. if(M.Dim(0)>0 && M.Dim(1)>0){
  2725. s_coord+=dof*M.Dim(0);
  2726. t_value=thread_buff[thread_id];
  2727. for(size_t j=0;j<dof*M.Dim(0);j++) t_value[j]=0;
  2728. }
  2729. size_t interac_cnt=0;
  2730. for(size_t j=0;j<trg_interac_cnt[i];j++){
  2731. if(ptr_single_layer_kernel!=(size_t)NULL){// Single layer kernel
  2732. Real_t* src_coord_=coord_data[0]+src_coord[i][2*j+0];
  2733. assert(thread_buff_size>=dof*M.Dim(0)+src_cnt[i][2*j+0]*COORD_DIM);
  2734. for(size_t k1=0;k1<src_cnt[i][2*j+0];k1++){ // Compute shifted source coordinates.
  2735. for(size_t k0=0;k0<COORD_DIM;k0++){
  2736. s_coord[k1*COORD_DIM+k0]=src_coord_[k1*COORD_DIM+k0]+shift_coord[i][j*COORD_DIM+k0];
  2737. }
  2738. }
  2739. single_layer_kernel( s_coord , src_cnt[i][2*j+0], input_data[0]+src_value[i][2*j+0], dof,
  2740. coord_data[0]+trg_coord[i], trg_cnt[i] , t_value);
  2741. interac_cnt+=src_cnt[i][2*j+0]*trg_cnt[i];
  2742. }
  2743. if(ptr_double_layer_kernel!=(size_t)NULL){// Double layer kernel
  2744. Real_t* src_coord_=coord_data[0]+src_coord[i][2*j+1];
  2745. assert(thread_buff_size>=dof*M.Dim(0)+src_cnt[i][2*j+1]*COORD_DIM);
  2746. for(size_t k1=0;k1<src_cnt[i][2*j+1];k1++){ // Compute shifted source coordinates.
  2747. for(size_t k0=0;k0<COORD_DIM;k0++){
  2748. s_coord[k1*COORD_DIM+k0]=src_coord_[k1*COORD_DIM+k0]+shift_coord[i][j*COORD_DIM+k0];
  2749. }
  2750. }
  2751. double_layer_kernel( s_coord , src_cnt[i][2*j+1], input_data[0]+src_value[i][2*j+1], dof,
  2752. coord_data[0]+trg_coord[i], trg_cnt[i] , t_value);
  2753. interac_cnt+=src_cnt[i][2*j+1]*trg_cnt[i];
  2754. }
  2755. }
  2756. if(M.Dim(0)>0 && M.Dim(1)>0 && interac_cnt>0){
  2757. assert(trg_cnt[i]*ker_dim1==M.Dim(0)); UNUSED(ker_dim1);
  2758. for(size_t j=0;j<dof*M.Dim(0);j++) t_value[j]*=scaling[i];
  2759. Matrix<Real_t> in_vec(dof, M.Dim(0), t_value , false);
  2760. Matrix<Real_t> out_vec(dof, M.Dim(1), output_data[0]+trg_value[i], false);
  2761. Matrix<Real_t>::DGEMM(out_vec, in_vec, M, 1.0);
  2762. }
  2763. }
  2764. }
  2765. #ifdef __INTEL_OFFLOAD
  2766. if(device) MIC_Lock::release_lock(lock_idx);
  2767. #endif
  2768. }
  2769. #ifndef __MIC_ASYNCH__
  2770. #ifdef __INTEL_OFFLOAD
  2771. #pragma offload if(device) target(mic:0)
  2772. {if(device) MIC_Lock::wait_lock(lock_idx);}
  2773. #endif
  2774. #endif
  2775. Profile::Toc();
  2776. }
  2777. template <class FMMNode>
  2778. void FMM_Pts<FMMNode>::X_ListSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  2779. if(this->MultipoleOrder()==0) return;
  2780. { // Set setup_data
  2781. setup_data.level=level;
  2782. setup_data.kernel=&aux_kernel;
  2783. setup_data.interac_type.resize(1);
  2784. setup_data.interac_type[0]=X_Type;
  2785. setup_data. input_data=&buff[4];
  2786. setup_data.output_data=&buff[1];
  2787. setup_data. coord_data=&buff[6];
  2788. Vector<FMMNode_t*>& nodes_in =n_list[4];
  2789. Vector<FMMNode_t*>& nodes_out=n_list[1];
  2790. setup_data.nodes_in .clear();
  2791. setup_data.nodes_out.clear();
  2792. for(size_t i=0;i<nodes_in .Dim();i++) if(nodes_in [i]->Depth()==level-1 || level==-1) setup_data.nodes_in .push_back(nodes_in [i]);
  2793. for(size_t i=0;i<nodes_out.Dim();i++) if(nodes_out[i]->Depth()==level || level==-1) setup_data.nodes_out.push_back(nodes_out[i]);
  2794. }
  2795. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  2796. std::vector<void*>& nodes_out=setup_data.nodes_out;
  2797. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  2798. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  2799. for(size_t i=0;i<nodes_in .size();i++){
  2800. input_vector .push_back(&((FMMNode*)nodes_in [i])->src_coord);
  2801. input_vector .push_back(&((FMMNode*)nodes_in [i])->src_value);
  2802. input_vector .push_back(&((FMMNode*)nodes_in [i])->surf_coord);
  2803. input_vector .push_back(&((FMMNode*)nodes_in [i])->surf_value);
  2804. }
  2805. for(size_t i=0;i<nodes_out.size();i++){
  2806. output_vector.push_back(&dnwd_check_surf[((FMMNode*)nodes_out[i])->Depth()]);
  2807. output_vector.push_back(&((FMMData*)((FMMNode*)nodes_out[i])->FMMData())->dnward_equiv);
  2808. }
  2809. //Downward check to downward equivalent matrix.
  2810. Matrix<Real_t>& M_dc2de = this->mat->Mat(level, DC2DE_Type, 0);
  2811. this->SetupInteracPts(setup_data, false, true, &M_dc2de,device);
  2812. { // Resize device buffer
  2813. size_t n=setup_data.output_data->Dim(0)*setup_data.output_data->Dim(1)*sizeof(Real_t);
  2814. if(this->dev_buffer.Dim()<n) this->dev_buffer.Resize(n);
  2815. }
  2816. }
  2817. template <class FMMNode>
  2818. void FMM_Pts<FMMNode>::X_List (SetupData<Real_t>& setup_data, bool device){
  2819. //Add X_List contribution.
  2820. this->EvalListPts(setup_data, device);
  2821. }
  2822. template <class FMMNode>
  2823. void FMM_Pts<FMMNode>::W_ListSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  2824. if(this->MultipoleOrder()==0) return;
  2825. { // Set setup_data
  2826. setup_data.level=level;
  2827. setup_data.kernel=&kernel;
  2828. setup_data.interac_type.resize(1);
  2829. setup_data.interac_type[0]=W_Type;
  2830. setup_data. input_data=&buff[0];
  2831. setup_data.output_data=&buff[5];
  2832. setup_data. coord_data=&buff[6];
  2833. Vector<FMMNode_t*>& nodes_in =n_list[0];
  2834. Vector<FMMNode_t*>& nodes_out=n_list[5];
  2835. setup_data.nodes_in .clear();
  2836. setup_data.nodes_out.clear();
  2837. for(size_t i=0;i<nodes_in .Dim();i++) if(nodes_in [i]->Depth()==level+1 || level==-1) setup_data.nodes_in .push_back(nodes_in [i]);
  2838. for(size_t i=0;i<nodes_out.Dim();i++) if(nodes_out[i]->Depth()==level || level==-1) setup_data.nodes_out.push_back(nodes_out[i]);
  2839. }
  2840. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  2841. std::vector<void*>& nodes_out=setup_data.nodes_out;
  2842. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  2843. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  2844. for(size_t i=0;i<nodes_in .size();i++){
  2845. input_vector .push_back(&upwd_equiv_surf[((FMMNode*)nodes_in [i])->Depth()]);
  2846. input_vector .push_back(&((FMMData*)((FMMNode*)nodes_in [i])->FMMData())->upward_equiv);
  2847. input_vector .push_back(NULL);
  2848. input_vector .push_back(NULL);
  2849. }
  2850. for(size_t i=0;i<nodes_out.size();i++){
  2851. output_vector.push_back(&((FMMNode*)nodes_out[i])->trg_coord);
  2852. output_vector.push_back(&((FMMNode*)nodes_out[i])->trg_value);
  2853. }
  2854. this->SetupInteracPts(setup_data, true, false, NULL, device);
  2855. { // Resize device buffer
  2856. size_t n=setup_data.output_data->Dim(0)*setup_data.output_data->Dim(1)*sizeof(Real_t);
  2857. if(this->dev_buffer.Dim()<n) this->dev_buffer.Resize(n);
  2858. }
  2859. }
  2860. template <class FMMNode>
  2861. void FMM_Pts<FMMNode>::W_List (SetupData<Real_t>& setup_data, bool device){
  2862. //Add W_List contribution.
  2863. this->EvalListPts(setup_data, device);
  2864. }
  2865. template <class FMMNode>
  2866. void FMM_Pts<FMMNode>::U_ListSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  2867. { // Set setup_data
  2868. setup_data.level=level;
  2869. setup_data.kernel=&kernel;
  2870. setup_data.interac_type.resize(3);
  2871. setup_data.interac_type[0]=U0_Type;
  2872. setup_data.interac_type[1]=U1_Type;
  2873. setup_data.interac_type[2]=U2_Type;
  2874. setup_data. input_data=&buff[4];
  2875. setup_data.output_data=&buff[5];
  2876. setup_data. coord_data=&buff[6];
  2877. Vector<FMMNode_t*>& nodes_in =n_list[4];
  2878. Vector<FMMNode_t*>& nodes_out=n_list[5];
  2879. setup_data.nodes_in .clear();
  2880. setup_data.nodes_out.clear();
  2881. for(size_t i=0;i<nodes_in .Dim();i++) if((level-1<=nodes_in [i]->Depth() && nodes_in [i]->Depth()<=level+1) || level==-1) setup_data.nodes_in .push_back(nodes_in [i]);
  2882. for(size_t i=0;i<nodes_out.Dim();i++) if(( nodes_out[i]->Depth()==level ) || level==-1) setup_data.nodes_out.push_back(nodes_out[i]);
  2883. }
  2884. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  2885. std::vector<void*>& nodes_out=setup_data.nodes_out;
  2886. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  2887. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  2888. for(size_t i=0;i<nodes_in .size();i++){
  2889. input_vector .push_back(&((FMMNode*)nodes_in [i])->src_coord);
  2890. input_vector .push_back(&((FMMNode*)nodes_in [i])->src_value);
  2891. input_vector .push_back(&((FMMNode*)nodes_in [i])->surf_coord);
  2892. input_vector .push_back(&((FMMNode*)nodes_in [i])->surf_value);
  2893. }
  2894. for(size_t i=0;i<nodes_out.size();i++){
  2895. output_vector.push_back(&((FMMNode*)nodes_out[i])->trg_coord);
  2896. output_vector.push_back(&((FMMNode*)nodes_out[i])->trg_value);
  2897. }
  2898. this->SetupInteracPts(setup_data, false, false, NULL, device);
  2899. { // Resize device buffer
  2900. size_t n=setup_data.output_data->Dim(0)*setup_data.output_data->Dim(1)*sizeof(Real_t);
  2901. if(this->dev_buffer.Dim()<n) this->dev_buffer.Resize(n);
  2902. }
  2903. }
  2904. template <class FMMNode>
  2905. void FMM_Pts<FMMNode>::U_List (SetupData<Real_t>& setup_data, bool device){
  2906. //Add U_List contribution.
  2907. this->EvalListPts(setup_data, device);
  2908. }
  2909. template <class FMMNode>
  2910. void FMM_Pts<FMMNode>::Down2TargetSetup(SetupData<Real_t>& setup_data, std::vector<Matrix<Real_t> >& buff, std::vector<Vector<FMMNode_t*> >& n_list, int level, bool device){
  2911. if(this->MultipoleOrder()==0) return;
  2912. { // Set setup_data
  2913. setup_data.level=level;
  2914. setup_data.kernel=&kernel;
  2915. setup_data.interac_type.resize(1);
  2916. setup_data.interac_type[0]=D2T_Type;
  2917. setup_data. input_data=&buff[1];
  2918. setup_data.output_data=&buff[5];
  2919. setup_data. coord_data=&buff[6];
  2920. Vector<FMMNode_t*>& nodes_in =n_list[1];
  2921. Vector<FMMNode_t*>& nodes_out=n_list[5];
  2922. setup_data.nodes_in .clear();
  2923. setup_data.nodes_out.clear();
  2924. for(size_t i=0;i<nodes_in .Dim();i++) if(nodes_in [i]->Depth()==level || level==-1) setup_data.nodes_in .push_back(nodes_in [i]);
  2925. for(size_t i=0;i<nodes_out.Dim();i++) if(nodes_out[i]->Depth()==level || level==-1) setup_data.nodes_out.push_back(nodes_out[i]);
  2926. }
  2927. std::vector<void*>& nodes_in =setup_data.nodes_in ;
  2928. std::vector<void*>& nodes_out=setup_data.nodes_out;
  2929. std::vector<Vector<Real_t>*>& input_vector=setup_data. input_vector; input_vector.clear();
  2930. std::vector<Vector<Real_t>*>& output_vector=setup_data.output_vector; output_vector.clear();
  2931. for(size_t i=0;i<nodes_in .size();i++){
  2932. input_vector .push_back(&dnwd_equiv_surf[((FMMNode*)nodes_in [i])->Depth()]);
  2933. input_vector .push_back(&((FMMData*)((FMMNode*)nodes_in [i])->FMMData())->dnward_equiv);
  2934. input_vector .push_back(NULL);
  2935. input_vector .push_back(NULL);
  2936. }
  2937. for(size_t i=0;i<nodes_out.size();i++){
  2938. output_vector.push_back(&((FMMNode*)nodes_out[i])->trg_coord);
  2939. output_vector.push_back(&((FMMNode*)nodes_out[i])->trg_value);
  2940. }
  2941. this->SetupInteracPts(setup_data, true, false, NULL, device);
  2942. { // Resize device buffer
  2943. size_t n=setup_data.output_data->Dim(0)*setup_data.output_data->Dim(1)*sizeof(Real_t);
  2944. if(this->dev_buffer.Dim()<n) this->dev_buffer.Resize(n);
  2945. }
  2946. }
  2947. template <class FMMNode>
  2948. void FMM_Pts<FMMNode>::Down2Target(SetupData<Real_t>& setup_data, bool device){
  2949. //Add Down2Target contribution.
  2950. this->EvalListPts(setup_data, device);
  2951. }
  2952. template <class FMMNode>
  2953. void FMM_Pts<FMMNode>::PostProcessing(std::vector<FMMNode_t*>& nodes){
  2954. }
  2955. template <class FMMNode>
  2956. void FMM_Pts<FMMNode>::CopyOutput(FMMNode** nodes, size_t n){
  2957. // for(size_t i=0;i<n;i++){
  2958. // FMMData* fmm_data=((FMMData*)nodes[i]->FMMData());
  2959. // }
  2960. }
  2961. }//end namespace