kernel.txx 91 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483248424852486248724882489249024912492249324942495249624972498249925002501250225032504250525062507250825092510251125122513251425152516251725182519252025212522252325242525252625272528252925302531253225332534253525362537253825392540254125422543254425452546254725482549255025512552255325542555255625572558255925602561256225632564256525662567256825692570257125722573257425752576257725782579258025812582258325842585258625872588258925902591259225932594259525962597259825992600260126022603260426052606260726082609261026112612261326142615261626172618261926202621262226232624262526262627262826292630263126322633263426352636
  1. /**
  2. * \file kernel.txx
  3. * \author Dhairya Malhotra, dhairya.malhotra@gmail.com
  4. * \date 12-20-2011
  5. * \brief This file contains the implementation of the struct Kernel and also the
  6. * implementation of various kernels for FMM.
  7. */
  8. #include <cmath>
  9. #include <cstdlib>
  10. #include <vector>
  11. #include <mem_mgr.hpp>
  12. #include <profile.hpp>
  13. #include <vector.hpp>
  14. #include <matrix.hpp>
  15. #include <precomp_mat.hpp>
  16. #include <intrin_wrapper.hpp>
  17. #include <cheb_utils.hpp>
  18. namespace pvfmm{
  19. /**
  20. * \brief Constructor.
  21. */
  22. template <class T>
  23. Kernel<T>::Kernel(Ker_t poten, Ker_t dbl_poten, const char* name, int dim_, std::pair<int,int> k_dim,
  24. size_t dev_poten, size_t dev_dbl_poten){
  25. dim=dim_;
  26. ker_dim[0]=k_dim.first;
  27. ker_dim[1]=k_dim.second;
  28. ker_poten=poten;
  29. dbl_layer_poten=dbl_poten;
  30. ker_name=std::string(name);
  31. dev_ker_poten=dev_poten;
  32. dev_dbl_layer_poten=dev_dbl_poten;
  33. k_s2m=NULL;
  34. k_s2l=NULL;
  35. k_s2t=NULL;
  36. k_m2m=NULL;
  37. k_m2l=NULL;
  38. k_m2t=NULL;
  39. k_l2l=NULL;
  40. k_l2t=NULL;
  41. vol_poten=NULL;
  42. scale_invar=true;
  43. src_scal.Resize(ker_dim[0]); src_scal.SetZero();
  44. trg_scal.Resize(ker_dim[1]); trg_scal.SetZero();
  45. perm_vec.Resize(Perm_Count);
  46. for(size_t p_type=0;p_type<C_Perm;p_type++){
  47. perm_vec[p_type ]=Permutation<T>(ker_dim[0]);
  48. perm_vec[p_type+C_Perm]=Permutation<T>(ker_dim[1]);
  49. }
  50. init=false;
  51. }
  52. /**
  53. * \brief Initialize the kernel.
  54. */
  55. template <class T>
  56. void Kernel<T>::Initialize(bool verbose) const{
  57. if(init) return;
  58. init=true;
  59. T eps=1.0;
  60. while(eps+(T)1.0>1.0) eps*=0.5;
  61. T scal=1.0;
  62. if(ker_dim[0]*ker_dim[1]>0){ // Determine scaling
  63. Matrix<T> M_scal(ker_dim[0],ker_dim[1]);
  64. size_t N=1024;
  65. T eps_=N*eps;
  66. T src_coord[3]={0,0,0};
  67. std::vector<T> trg_coord1(N*COORD_DIM);
  68. Matrix<T> M1(N,ker_dim[0]*ker_dim[1]);
  69. while(true){
  70. T abs_sum=0;
  71. for(size_t i=0;i<N/2;i++){
  72. T x,y,z,r;
  73. do{
  74. x=(drand48()-0.5);
  75. y=(drand48()-0.5);
  76. z=(drand48()-0.5);
  77. r=pvfmm::sqrt<T>(x*x+y*y+z*z);
  78. }while(r<0.25);
  79. trg_coord1[i*COORD_DIM+0]=x*scal;
  80. trg_coord1[i*COORD_DIM+1]=y*scal;
  81. trg_coord1[i*COORD_DIM+2]=z*scal;
  82. }
  83. for(size_t i=N/2;i<N;i++){
  84. T x,y,z,r;
  85. do{
  86. x=(drand48()-0.5);
  87. y=(drand48()-0.5);
  88. z=(drand48()-0.5);
  89. r=pvfmm::sqrt<T>(x*x+y*y+z*z);
  90. }while(r<0.25);
  91. trg_coord1[i*COORD_DIM+0]=x*1.0/scal;
  92. trg_coord1[i*COORD_DIM+1]=y*1.0/scal;
  93. trg_coord1[i*COORD_DIM+2]=z*1.0/scal;
  94. }
  95. for(size_t i=0;i<N;i++){
  96. BuildMatrix(&src_coord [ 0], 1,
  97. &trg_coord1[i*COORD_DIM], 1, &(M1[i][0]));
  98. for(size_t j=0;j<ker_dim[0]*ker_dim[1];j++){
  99. abs_sum+=pvfmm::fabs<T>(M1[i][j]);
  100. }
  101. }
  102. if(abs_sum>pvfmm::sqrt<T>(eps) || scal<eps) break;
  103. scal=scal*0.5;
  104. }
  105. std::vector<T> trg_coord2(N*COORD_DIM);
  106. Matrix<T> M2(N,ker_dim[0]*ker_dim[1]);
  107. for(size_t i=0;i<N*COORD_DIM;i++){
  108. trg_coord2[i]=trg_coord1[i]*0.5;
  109. }
  110. for(size_t i=0;i<N;i++){
  111. BuildMatrix(&src_coord [ 0], 1,
  112. &trg_coord2[i*COORD_DIM], 1, &(M2[i][0]));
  113. }
  114. for(size_t i=0;i<ker_dim[0]*ker_dim[1];i++){
  115. T dot11=0, dot12=0, dot22=0;
  116. for(size_t j=0;j<N;j++){
  117. dot11+=M1[j][i]*M1[j][i];
  118. dot12+=M1[j][i]*M2[j][i];
  119. dot22+=M2[j][i]*M2[j][i];
  120. }
  121. T max_val=std::max<T>(dot11,dot22);
  122. if(dot11>max_val*eps &&
  123. dot22>max_val*eps ){
  124. T s=dot12/dot11;
  125. M_scal[0][i]=pvfmm::log<T>(s)/pvfmm::log<T>(2.0);
  126. T err=pvfmm::sqrt<T>(0.5*(dot22/dot11)/(s*s)-0.5);
  127. if(err>eps_){
  128. scale_invar=false;
  129. M_scal[0][i]=0.0;
  130. }
  131. //assert(M_scal[0][i]>=0.0); // Kernel function must decay
  132. }else if(dot11>max_val*eps ||
  133. dot22>max_val*eps ){
  134. scale_invar=false;
  135. M_scal[0][i]=0.0;
  136. }else{
  137. M_scal[0][i]=-1;
  138. }
  139. }
  140. src_scal.Resize(ker_dim[0]); src_scal.SetZero();
  141. trg_scal.Resize(ker_dim[1]); trg_scal.SetZero();
  142. if(scale_invar){
  143. Matrix<T> b(ker_dim[0]*ker_dim[1]+1,1); b.SetZero();
  144. mem::memcopy(&b[0][0],&M_scal[0][0],ker_dim[0]*ker_dim[1]*sizeof(T));
  145. Matrix<T> M(ker_dim[0]*ker_dim[1]+1,ker_dim[0]+ker_dim[1]); M.SetZero();
  146. M[ker_dim[0]*ker_dim[1]][0]=1;
  147. for(size_t i0=0;i0<ker_dim[0];i0++)
  148. for(size_t i1=0;i1<ker_dim[1];i1++){
  149. size_t j=i0*ker_dim[1]+i1;
  150. if(fabs(b[j][0])>=0){
  151. M[j][ 0+ i0]=1;
  152. M[j][i1+ker_dim[0]]=1;
  153. }
  154. }
  155. Matrix<T> x=M.pinv()*b;
  156. for(size_t i=0;i<ker_dim[0];i++){
  157. src_scal[i]=x[i][0];
  158. }
  159. for(size_t i=0;i<ker_dim[1];i++){
  160. trg_scal[i]=x[ker_dim[0]+i][0];
  161. }
  162. for(size_t i0=0;i0<ker_dim[0];i0++)
  163. for(size_t i1=0;i1<ker_dim[1];i1++){
  164. if(M_scal[i0][i1]>=0){
  165. if(pvfmm::fabs<T>(src_scal[i0]+trg_scal[i1]-M_scal[i0][i1])>eps_){
  166. scale_invar=false;
  167. }
  168. }
  169. }
  170. }
  171. if(!scale_invar){
  172. src_scal.SetZero();
  173. trg_scal.SetZero();
  174. //std::cout<<ker_name<<" not-scale-invariant\n";
  175. }
  176. }
  177. if(ker_dim[0]*ker_dim[1]>0){ // Determine symmetry
  178. size_t N=1024;
  179. T eps_=N*eps;
  180. T src_coord[3]={0,0,0};
  181. std::vector<T> trg_coord1(N*COORD_DIM);
  182. std::vector<T> trg_coord2(N*COORD_DIM);
  183. for(size_t i=0;i<N/2;i++){
  184. T x,y,z,r;
  185. do{
  186. x=(drand48()-0.5);
  187. y=(drand48()-0.5);
  188. z=(drand48()-0.5);
  189. r=pvfmm::sqrt<T>(x*x+y*y+z*z);
  190. }while(r<0.25);
  191. trg_coord1[i*COORD_DIM+0]=x*scal;
  192. trg_coord1[i*COORD_DIM+1]=y*scal;
  193. trg_coord1[i*COORD_DIM+2]=z*scal;
  194. }
  195. for(size_t i=N/2;i<N;i++){
  196. T x,y,z,r;
  197. do{
  198. x=(drand48()-0.5);
  199. y=(drand48()-0.5);
  200. z=(drand48()-0.5);
  201. r=pvfmm::sqrt<T>(x*x+y*y+z*z);
  202. }while(r<0.25);
  203. trg_coord1[i*COORD_DIM+0]=x*1.0/scal;
  204. trg_coord1[i*COORD_DIM+1]=y*1.0/scal;
  205. trg_coord1[i*COORD_DIM+2]=z*1.0/scal;
  206. }
  207. for(size_t p_type=0;p_type<C_Perm;p_type++){ // For each symmetry transform
  208. switch(p_type){ // Set trg_coord2
  209. case ReflecX:
  210. for(size_t i=0;i<N;i++){
  211. trg_coord2[i*COORD_DIM+0]=-trg_coord1[i*COORD_DIM+0];
  212. trg_coord2[i*COORD_DIM+1]= trg_coord1[i*COORD_DIM+1];
  213. trg_coord2[i*COORD_DIM+2]= trg_coord1[i*COORD_DIM+2];
  214. }
  215. break;
  216. case ReflecY:
  217. for(size_t i=0;i<N;i++){
  218. trg_coord2[i*COORD_DIM+0]= trg_coord1[i*COORD_DIM+0];
  219. trg_coord2[i*COORD_DIM+1]=-trg_coord1[i*COORD_DIM+1];
  220. trg_coord2[i*COORD_DIM+2]= trg_coord1[i*COORD_DIM+2];
  221. }
  222. break;
  223. case ReflecZ:
  224. for(size_t i=0;i<N;i++){
  225. trg_coord2[i*COORD_DIM+0]= trg_coord1[i*COORD_DIM+0];
  226. trg_coord2[i*COORD_DIM+1]= trg_coord1[i*COORD_DIM+1];
  227. trg_coord2[i*COORD_DIM+2]=-trg_coord1[i*COORD_DIM+2];
  228. }
  229. break;
  230. case SwapXY:
  231. for(size_t i=0;i<N;i++){
  232. trg_coord2[i*COORD_DIM+0]= trg_coord1[i*COORD_DIM+1];
  233. trg_coord2[i*COORD_DIM+1]= trg_coord1[i*COORD_DIM+0];
  234. trg_coord2[i*COORD_DIM+2]= trg_coord1[i*COORD_DIM+2];
  235. }
  236. break;
  237. case SwapXZ:
  238. for(size_t i=0;i<N;i++){
  239. trg_coord2[i*COORD_DIM+0]= trg_coord1[i*COORD_DIM+2];
  240. trg_coord2[i*COORD_DIM+1]= trg_coord1[i*COORD_DIM+1];
  241. trg_coord2[i*COORD_DIM+2]= trg_coord1[i*COORD_DIM+0];
  242. }
  243. break;
  244. default:
  245. for(size_t i=0;i<N;i++){
  246. trg_coord2[i*COORD_DIM+0]= trg_coord1[i*COORD_DIM+0];
  247. trg_coord2[i*COORD_DIM+1]= trg_coord1[i*COORD_DIM+1];
  248. trg_coord2[i*COORD_DIM+2]= trg_coord1[i*COORD_DIM+2];
  249. }
  250. }
  251. Matrix<long long> M11, M22;
  252. {
  253. Matrix<T> M1(N,ker_dim[0]*ker_dim[1]); M1.SetZero();
  254. Matrix<T> M2(N,ker_dim[0]*ker_dim[1]); M2.SetZero();
  255. for(size_t i=0;i<N;i++){
  256. BuildMatrix(&src_coord [ 0], 1,
  257. &trg_coord1[i*COORD_DIM], 1, &(M1[i][0]));
  258. BuildMatrix(&src_coord [ 0], 1,
  259. &trg_coord2[i*COORD_DIM], 1, &(M2[i][0]));
  260. }
  261. Matrix<T> dot11(ker_dim[0]*ker_dim[1],ker_dim[0]*ker_dim[1]);dot11.SetZero();
  262. Matrix<T> dot12(ker_dim[0]*ker_dim[1],ker_dim[0]*ker_dim[1]);dot12.SetZero();
  263. Matrix<T> dot22(ker_dim[0]*ker_dim[1],ker_dim[0]*ker_dim[1]);dot22.SetZero();
  264. std::vector<T> norm1(ker_dim[0]*ker_dim[1]);
  265. std::vector<T> norm2(ker_dim[0]*ker_dim[1]);
  266. {
  267. for(size_t k=0;k<N;k++)
  268. for(size_t i=0;i<ker_dim[0]*ker_dim[1];i++)
  269. for(size_t j=0;j<ker_dim[0]*ker_dim[1];j++){
  270. dot11[i][j]+=M1[k][i]*M1[k][j];
  271. dot12[i][j]+=M1[k][i]*M2[k][j];
  272. dot22[i][j]+=M2[k][i]*M2[k][j];
  273. }
  274. for(size_t i=0;i<ker_dim[0]*ker_dim[1];i++){
  275. norm1[i]=pvfmm::sqrt<T>(dot11[i][i]);
  276. norm2[i]=pvfmm::sqrt<T>(dot22[i][i]);
  277. }
  278. for(size_t i=0;i<ker_dim[0]*ker_dim[1];i++)
  279. for(size_t j=0;j<ker_dim[0]*ker_dim[1];j++){
  280. dot11[i][j]/=(norm1[i]*norm1[j]);
  281. dot12[i][j]/=(norm1[i]*norm2[j]);
  282. dot22[i][j]/=(norm2[i]*norm2[j]);
  283. }
  284. }
  285. long long flag=1;
  286. M11.Resize(ker_dim[0],ker_dim[1]); M11.SetZero();
  287. M22.Resize(ker_dim[0],ker_dim[1]); M22.SetZero();
  288. for(size_t i=0;i<ker_dim[0]*ker_dim[1];i++){
  289. if(norm1[i]>eps_ && M11[0][i]==0){
  290. for(size_t j=0;j<ker_dim[0]*ker_dim[1];j++){
  291. if(pvfmm::fabs<T>(norm1[i]-norm1[j])<eps_ && pvfmm::fabs<T>(pvfmm::fabs<T>(dot11[i][j])-1.0)<eps_){
  292. M11[0][j]=(dot11[i][j]>0?flag:-flag);
  293. }
  294. if(pvfmm::fabs<T>(norm1[i]-norm2[j])<eps_ && pvfmm::fabs<T>(pvfmm::fabs<T>(dot12[i][j])-1.0)<eps_){
  295. M22[0][j]=(dot12[i][j]>0?flag:-flag);
  296. }
  297. }
  298. flag++;
  299. }
  300. }
  301. }
  302. Matrix<long long> P1, P2;
  303. { // P1
  304. Matrix<long long>& P=P1;
  305. Matrix<long long> M1=M11;
  306. Matrix<long long> M2=M22;
  307. for(size_t i=0;i<M1.Dim(0);i++){
  308. for(size_t j=0;j<M1.Dim(1);j++){
  309. if(M1[i][j]<0) M1[i][j]=-M1[i][j];
  310. if(M2[i][j]<0) M2[i][j]=-M2[i][j];
  311. }
  312. std::sort(&M1[i][0],&M1[i][M1.Dim(1)]);
  313. std::sort(&M2[i][0],&M2[i][M2.Dim(1)]);
  314. }
  315. P.Resize(M1.Dim(0),M1.Dim(0));
  316. for(size_t i=0;i<M1.Dim(0);i++)
  317. for(size_t j=0;j<M1.Dim(0);j++){
  318. P[i][j]=1;
  319. for(size_t k=0;k<M1.Dim(1);k++)
  320. if(M1[i][k]!=M2[j][k]){
  321. P[i][j]=0;
  322. break;
  323. }
  324. }
  325. }
  326. { // P2
  327. Matrix<long long>& P=P2;
  328. Matrix<long long> M1=M11.Transpose();
  329. Matrix<long long> M2=M22.Transpose();
  330. for(size_t i=0;i<M1.Dim(0);i++){
  331. for(size_t j=0;j<M1.Dim(1);j++){
  332. if(M1[i][j]<0) M1[i][j]=-M1[i][j];
  333. if(M2[i][j]<0) M2[i][j]=-M2[i][j];
  334. }
  335. std::sort(&M1[i][0],&M1[i][M1.Dim(1)]);
  336. std::sort(&M2[i][0],&M2[i][M2.Dim(1)]);
  337. }
  338. P.Resize(M1.Dim(0),M1.Dim(0));
  339. for(size_t i=0;i<M1.Dim(0);i++)
  340. for(size_t j=0;j<M1.Dim(0);j++){
  341. P[i][j]=1;
  342. for(size_t k=0;k<M1.Dim(1);k++)
  343. if(M1[i][k]!=M2[j][k]){
  344. P[i][j]=0;
  345. break;
  346. }
  347. }
  348. }
  349. std::vector<Permutation<long long> > P1vec, P2vec;
  350. { // P1vec
  351. Matrix<long long>& Pmat=P1;
  352. std::vector<Permutation<long long> >& Pvec=P1vec;
  353. Permutation<long long> P(Pmat.Dim(0));
  354. Vector<PERM_INT_T>& perm=P.perm;
  355. perm.SetZero();
  356. // First permutation
  357. for(size_t i=0;i<P.Dim();i++)
  358. for(size_t j=0;j<P.Dim();j++){
  359. if(Pmat[i][j]){
  360. perm[i]=j;
  361. break;
  362. }
  363. }
  364. Vector<PERM_INT_T> perm_tmp;
  365. while(true){ // Next permutation
  366. perm_tmp=perm;
  367. std::sort(&perm_tmp[0],&perm_tmp[0]+perm_tmp.Dim());
  368. for(size_t i=0;i<perm_tmp.Dim();i++){
  369. if(perm_tmp[i]!=i) break;
  370. if(i==perm_tmp.Dim()-1){
  371. Pvec.push_back(P);
  372. }
  373. }
  374. bool last=false;
  375. for(size_t i=0;i<P.Dim();i++){
  376. PERM_INT_T tmp=perm[i];
  377. for(size_t j=perm[i]+1;j<P.Dim();j++){
  378. if(Pmat[i][j]){
  379. perm[i]=j;
  380. break;
  381. }
  382. }
  383. if(perm[i]>tmp) break;
  384. for(size_t j=0;j<P.Dim();j++){
  385. if(Pmat[i][j]){
  386. perm[i]=j;
  387. break;
  388. }
  389. }
  390. if(i==P.Dim()-1) last=true;
  391. }
  392. if(last) break;
  393. }
  394. }
  395. { // P2vec
  396. Matrix<long long>& Pmat=P2;
  397. std::vector<Permutation<long long> >& Pvec=P2vec;
  398. Permutation<long long> P(Pmat.Dim(0));
  399. Vector<PERM_INT_T>& perm=P.perm;
  400. perm.SetZero();
  401. // First permutation
  402. for(size_t i=0;i<P.Dim();i++)
  403. for(size_t j=0;j<P.Dim();j++){
  404. if(Pmat[i][j]){
  405. perm[i]=j;
  406. break;
  407. }
  408. }
  409. Vector<PERM_INT_T> perm_tmp;
  410. while(true){ // Next permutation
  411. perm_tmp=perm;
  412. std::sort(&perm_tmp[0],&perm_tmp[0]+perm_tmp.Dim());
  413. for(size_t i=0;i<perm_tmp.Dim();i++){
  414. if(perm_tmp[i]!=i) break;
  415. if(i==perm_tmp.Dim()-1){
  416. Pvec.push_back(P);
  417. }
  418. }
  419. bool last=false;
  420. for(size_t i=0;i<P.Dim();i++){
  421. PERM_INT_T tmp=perm[i];
  422. for(size_t j=perm[i]+1;j<P.Dim();j++){
  423. if(Pmat[i][j]){
  424. perm[i]=j;
  425. break;
  426. }
  427. }
  428. if(perm[i]>tmp) break;
  429. for(size_t j=0;j<P.Dim();j++){
  430. if(Pmat[i][j]){
  431. perm[i]=j;
  432. break;
  433. }
  434. }
  435. if(i==P.Dim()-1) last=true;
  436. }
  437. if(last) break;
  438. }
  439. }
  440. { // Find pairs which acutally work (neglect scaling)
  441. std::vector<Permutation<long long> > P1vec_, P2vec_;
  442. Matrix<long long> M1=M11;
  443. Matrix<long long> M2=M22;
  444. for(size_t i=0;i<M1.Dim(0);i++){
  445. for(size_t j=0;j<M1.Dim(1);j++){
  446. if(M1[i][j]<0) M1[i][j]=-M1[i][j];
  447. if(M2[i][j]<0) M2[i][j]=-M2[i][j];
  448. }
  449. }
  450. Matrix<long long> M;
  451. for(size_t i=0;i<P1vec.size();i++)
  452. for(size_t j=0;j<P2vec.size();j++){
  453. M=P1vec[i]*M2*P2vec[j];
  454. for(size_t k=0;k<M.Dim(0)*M.Dim(1);k++){
  455. if(M[0][k]!=M1[0][k]) break;
  456. if(k==M.Dim(0)*M.Dim(1)-1){
  457. P1vec_.push_back(P1vec[i]);
  458. P2vec_.push_back(P2vec[j]);
  459. }
  460. }
  461. }
  462. P1vec=P1vec_;
  463. P2vec=P2vec_;
  464. }
  465. Permutation<T> P1_, P2_;
  466. { // Find pairs which acutally work
  467. for(size_t k=0;k<P1vec.size();k++){
  468. Permutation<long long> P1=P1vec[k];
  469. Permutation<long long> P2=P2vec[k];
  470. Matrix<long long> M1= M11 ;
  471. Matrix<long long> M2=P1*M22*P2;
  472. Matrix<T> M(M1.Dim(0)*M1.Dim(1)+1,M1.Dim(0)+M1.Dim(1));
  473. M.SetZero(); M[M1.Dim(0)*M1.Dim(1)][0]=1.0;
  474. for(size_t i=0;i<M1.Dim(0);i++)
  475. for(size_t j=0;j<M1.Dim(1);j++){
  476. size_t k=i*M1.Dim(1)+j;
  477. M[k][ i]= M1[i][j];
  478. M[k][M1.Dim(0)+j]=-M2[i][j];
  479. }
  480. M=M.pinv();
  481. { // Construct new permutation
  482. Permutation<long long> P1_(M1.Dim(0));
  483. Permutation<long long> P2_(M1.Dim(1));
  484. for(size_t i=0;i<M1.Dim(0);i++){
  485. P1_.scal[i]=(M[i][M1.Dim(0)*M1.Dim(1)]>0?1:-1);
  486. }
  487. for(size_t i=0;i<M1.Dim(1);i++){
  488. P2_.scal[i]=(M[M1.Dim(0)+i][M1.Dim(0)*M1.Dim(1)]>0?1:-1);
  489. }
  490. P1=P1_*P1 ;
  491. P2=P2 *P2_;
  492. }
  493. bool done=true;
  494. Matrix<long long> Merr=P1*M22*P2-M11;
  495. for(size_t i=0;i<Merr.Dim(0)*Merr.Dim(1);i++){
  496. if(Merr[0][i]){
  497. done=false;
  498. break;
  499. }
  500. }
  501. if(done){
  502. P1_=Permutation<T>(P1.Dim());
  503. P2_=Permutation<T>(P2.Dim());
  504. for(size_t i=0;i<P1.Dim();i++){
  505. P1_.perm[i]=P1.perm[i];
  506. P1_.scal[i]=P1.scal[i];
  507. }
  508. for(size_t i=0;i<P2.Dim();i++){
  509. P2_.perm[i]=P2.perm[i];
  510. P2_.scal[i]=P2.scal[i];
  511. }
  512. break;
  513. }
  514. }
  515. }
  516. //std::cout<<P1_<<'\n';
  517. //std::cout<<P2_<<'\n';
  518. perm_vec[p_type ]=P1_.Transpose();
  519. perm_vec[p_type+C_Perm]=P2_;
  520. }
  521. for(size_t i=0;i<2*C_Perm;i++){
  522. if(perm_vec[i].Dim()==0){
  523. perm_vec.Resize(0);
  524. std::cout<<"no-symmetry for: "<<ker_name<<'\n';
  525. break;
  526. }
  527. }
  528. }
  529. if(verbose){ // Display kernel information
  530. std::cout<<"\n";
  531. std::cout<<"Kernel Name : "<<ker_name<<'\n';
  532. std::cout<<"Precision : "<<(double)eps<<'\n';
  533. std::cout<<"Symmetry : "<<(perm_vec.Dim()>0?"yes":"no")<<'\n';
  534. std::cout<<"Scale Invariant: "<<(scale_invar?"yes":"no")<<'\n';
  535. if(scale_invar && ker_dim[0]*ker_dim[1]>0){
  536. std::cout<<"Scaling Matrix :\n";
  537. Matrix<T> Src(ker_dim[0],1);
  538. Matrix<T> Trg(1,ker_dim[1]);
  539. for(size_t i=0;i<ker_dim[0];i++) Src[i][0]=pvfmm::pow<T>(2.0,src_scal[i]);
  540. for(size_t i=0;i<ker_dim[1];i++) Trg[0][i]=pvfmm::pow<T>(2.0,trg_scal[i]);
  541. std::cout<<Src*Trg;
  542. }
  543. if(ker_dim[0]*ker_dim[1]>0){ // Accuracy of multipole expansion
  544. std::cout<<"Multipole Error: ";
  545. for(T rad=1.0; rad>1.0e-2; rad*=0.5){
  546. int m=8; // multipole order
  547. std::vector<T> equiv_surf;
  548. std::vector<T> check_surf;
  549. for(int i0=0;i0<m;i0++){
  550. for(int i1=0;i1<m;i1++){
  551. for(int i2=0;i2<m;i2++){
  552. if(i0== 0 || i1== 0 || i2== 0 ||
  553. i0==m-1 || i1==m-1 || i2==m-1){
  554. // Range: [-1/3,1/3]^3
  555. T x=((T)2*i0-(m-1))/(m-1)/3;
  556. T y=((T)2*i1-(m-1))/(m-1)/3;
  557. T z=((T)2*i2-(m-1))/(m-1)/3;
  558. equiv_surf.push_back(x*RAD0*rad);
  559. equiv_surf.push_back(y*RAD0*rad);
  560. equiv_surf.push_back(z*RAD0*rad);
  561. check_surf.push_back(x*RAD1*rad);
  562. check_surf.push_back(y*RAD1*rad);
  563. check_surf.push_back(z*RAD1*rad);
  564. }
  565. }
  566. }
  567. }
  568. size_t n_equiv=equiv_surf.size()/COORD_DIM;
  569. size_t n_check=equiv_surf.size()/COORD_DIM;
  570. size_t n_src=m*m;
  571. size_t n_trg=m*m;
  572. std::vector<T> src_coord;
  573. std::vector<T> trg_coord;
  574. for(size_t i=0;i<n_src*COORD_DIM;i++){
  575. src_coord.push_back((2*drand48()-1)/3*rad);
  576. }
  577. for(size_t i=0;i<n_trg;i++){
  578. T x,y,z,r;
  579. do{
  580. x=(drand48()-0.5);
  581. y=(drand48()-0.5);
  582. z=(drand48()-0.5);
  583. r=pvfmm::sqrt<T>(x*x+y*y+z*z);
  584. }while(r==0.0);
  585. trg_coord.push_back(x/r*pvfmm::sqrt<T>((T)COORD_DIM)*rad*(1.0+drand48()));
  586. trg_coord.push_back(y/r*pvfmm::sqrt<T>((T)COORD_DIM)*rad*(1.0+drand48()));
  587. trg_coord.push_back(z/r*pvfmm::sqrt<T>((T)COORD_DIM)*rad*(1.0+drand48()));
  588. }
  589. Matrix<T> M_s2c(n_src*ker_dim[0],n_check*ker_dim[1]);
  590. BuildMatrix( &src_coord[0], n_src,
  591. &check_surf[0], n_check, &(M_s2c[0][0]));
  592. Matrix<T> M_e2c(n_equiv*ker_dim[0],n_check*ker_dim[1]);
  593. BuildMatrix(&equiv_surf[0], n_equiv,
  594. &check_surf[0], n_check, &(M_e2c[0][0]));
  595. Matrix<T> M_c2e0, M_c2e1;
  596. {
  597. Matrix<T> U,S,V;
  598. M_e2c.SVD(U,S,V);
  599. T eps=1, max_S=0;
  600. while(eps*(T)0.5+(T)1.0>1.0) eps*=0.5;
  601. for(size_t i=0;i<std::min(S.Dim(0),S.Dim(1));i++){
  602. if(pvfmm::fabs<T>(S[i][i])>max_S) max_S=pvfmm::fabs<T>(S[i][i]);
  603. }
  604. for(size_t i=0;i<S.Dim(0);i++) S[i][i]=(S[i][i]>eps*max_S*4?1.0/S[i][i]:0.0);
  605. M_c2e0=V.Transpose()*S;
  606. M_c2e1=U.Transpose();
  607. }
  608. Matrix<T> M_e2t(n_equiv*ker_dim[0],n_trg*ker_dim[1]);
  609. BuildMatrix(&equiv_surf[0], n_equiv,
  610. &trg_coord[0], n_trg , &(M_e2t[0][0]));
  611. Matrix<T> M_s2t(n_src*ker_dim[0],n_trg*ker_dim[1]);
  612. BuildMatrix( &src_coord[0], n_src,
  613. &trg_coord[0], n_trg , &(M_s2t[0][0]));
  614. Matrix<T> M=(M_s2c*M_c2e0)*(M_c2e1*M_e2t)-M_s2t;
  615. T max_error=0, max_value=0;
  616. for(size_t i=0;i<M.Dim(0);i++)
  617. for(size_t j=0;j<M.Dim(1);j++){
  618. max_error=std::max<T>(max_error,pvfmm::fabs<T>(M [i][j]));
  619. max_value=std::max<T>(max_value,pvfmm::fabs<T>(M_s2t[i][j]));
  620. }
  621. std::cout<<(double)(max_error/max_value)<<' ';
  622. if(scale_invar) break;
  623. }
  624. std::cout<<"\n";
  625. }
  626. if(ker_dim[0]*ker_dim[1]>0){ // Accuracy of local expansion
  627. std::cout<<"Local-exp Error: ";
  628. for(T rad=1.0; rad>1.0e-2; rad*=0.5){
  629. int m=8; // multipole order
  630. std::vector<T> equiv_surf;
  631. std::vector<T> check_surf;
  632. for(int i0=0;i0<m;i0++){
  633. for(int i1=0;i1<m;i1++){
  634. for(int i2=0;i2<m;i2++){
  635. if(i0== 0 || i1== 0 || i2== 0 ||
  636. i0==m-1 || i1==m-1 || i2==m-1){
  637. // Range: [-1/3,1/3]^3
  638. T x=((T)2*i0-(m-1))/(m-1)/3;
  639. T y=((T)2*i1-(m-1))/(m-1)/3;
  640. T z=((T)2*i2-(m-1))/(m-1)/3;
  641. equiv_surf.push_back(x*RAD1*rad);
  642. equiv_surf.push_back(y*RAD1*rad);
  643. equiv_surf.push_back(z*RAD1*rad);
  644. check_surf.push_back(x*RAD0*rad);
  645. check_surf.push_back(y*RAD0*rad);
  646. check_surf.push_back(z*RAD0*rad);
  647. }
  648. }
  649. }
  650. }
  651. size_t n_equiv=equiv_surf.size()/COORD_DIM;
  652. size_t n_check=equiv_surf.size()/COORD_DIM;
  653. size_t n_src=m*m;
  654. size_t n_trg=m*m;
  655. std::vector<T> src_coord;
  656. std::vector<T> trg_coord;
  657. for(size_t i=0;i<n_trg*COORD_DIM;i++){
  658. trg_coord.push_back((2*drand48()-1)/3*rad);
  659. }
  660. for(size_t i=0;i<n_src;i++){
  661. T x,y,z,r;
  662. do{
  663. x=(drand48()-0.5);
  664. y=(drand48()-0.5);
  665. z=(drand48()-0.5);
  666. r=pvfmm::sqrt<T>(x*x+y*y+z*z);
  667. }while(r==0.0);
  668. src_coord.push_back(x/r*pvfmm::sqrt<T>((T)COORD_DIM)*rad*(1.0+drand48()));
  669. src_coord.push_back(y/r*pvfmm::sqrt<T>((T)COORD_DIM)*rad*(1.0+drand48()));
  670. src_coord.push_back(z/r*pvfmm::sqrt<T>((T)COORD_DIM)*rad*(1.0+drand48()));
  671. }
  672. Matrix<T> M_s2c(n_src*ker_dim[0],n_check*ker_dim[1]);
  673. BuildMatrix( &src_coord[0], n_src,
  674. &check_surf[0], n_check, &(M_s2c[0][0]));
  675. Matrix<T> M_e2c(n_equiv*ker_dim[0],n_check*ker_dim[1]);
  676. BuildMatrix(&equiv_surf[0], n_equiv,
  677. &check_surf[0], n_check, &(M_e2c[0][0]));
  678. Matrix<T> M_c2e0, M_c2e1;
  679. {
  680. Matrix<T> U,S,V;
  681. M_e2c.SVD(U,S,V);
  682. T eps=1, max_S=0;
  683. while(eps*(T)0.5+(T)1.0>1.0) eps*=0.5;
  684. for(size_t i=0;i<std::min(S.Dim(0),S.Dim(1));i++){
  685. if(pvfmm::fabs<T>(S[i][i])>max_S) max_S=pvfmm::fabs<T>(S[i][i]);
  686. }
  687. for(size_t i=0;i<S.Dim(0);i++) S[i][i]=(S[i][i]>eps*max_S*4?1.0/S[i][i]:0.0);
  688. M_c2e0=V.Transpose()*S;
  689. M_c2e1=U.Transpose();
  690. }
  691. Matrix<T> M_e2t(n_equiv*ker_dim[0],n_trg*ker_dim[1]);
  692. BuildMatrix(&equiv_surf[0], n_equiv,
  693. &trg_coord[0], n_trg , &(M_e2t[0][0]));
  694. Matrix<T> M_s2t(n_src*ker_dim[0],n_trg*ker_dim[1]);
  695. BuildMatrix( &src_coord[0], n_src,
  696. &trg_coord[0], n_trg , &(M_s2t[0][0]));
  697. Matrix<T> M=(M_s2c*M_c2e0)*(M_c2e1*M_e2t)-M_s2t;
  698. T max_error=0, max_value=0;
  699. for(size_t i=0;i<M.Dim(0);i++)
  700. for(size_t j=0;j<M.Dim(1);j++){
  701. max_error=std::max<T>(max_error,pvfmm::fabs<T>(M [i][j]));
  702. max_value=std::max<T>(max_value,pvfmm::fabs<T>(M_s2t[i][j]));
  703. }
  704. std::cout<<(double)(max_error/max_value)<<' ';
  705. if(scale_invar) break;
  706. }
  707. std::cout<<"\n";
  708. }
  709. if(vol_poten && ker_dim[0]*ker_dim[1]>0){ // Check if the volume potential is consistent with integral of kernel.
  710. int m=8; // multipole order
  711. std::vector<T> equiv_surf;
  712. std::vector<T> check_surf;
  713. std::vector<T> trg_coord;
  714. for(size_t i=0;i<m*COORD_DIM;i++){
  715. trg_coord.push_back(drand48()+1.0);
  716. }
  717. for(int i0=0;i0<m;i0++){
  718. for(int i1=0;i1<m;i1++){
  719. for(int i2=0;i2<m;i2++){
  720. if(i0== 0 || i1== 0 || i2== 0 ||
  721. i0==m-1 || i1==m-1 || i2==m-1){
  722. // Range: [-1/2,1/2]^3
  723. T x=((T)2*i0-(m-1))/(m-1)/2;
  724. T y=((T)2*i1-(m-1))/(m-1)/2;
  725. T z=((T)2*i2-(m-1))/(m-1)/2;
  726. equiv_surf.push_back(x*RAD1+1.5);
  727. equiv_surf.push_back(y*RAD1+1.5);
  728. equiv_surf.push_back(z*RAD1+1.5);
  729. check_surf.push_back(x*RAD0+1.5);
  730. check_surf.push_back(y*RAD0+1.5);
  731. check_surf.push_back(z*RAD0+1.5);
  732. }
  733. }
  734. }
  735. }
  736. size_t n_equiv=equiv_surf.size()/COORD_DIM;
  737. size_t n_check=equiv_surf.size()/COORD_DIM;
  738. size_t n_trg =trg_coord .size()/COORD_DIM;
  739. Matrix<T> M_local, M_analytic;
  740. Matrix<T> T_local, T_analytic;
  741. { // Compute local expansions M_local, T_local
  742. Matrix<T> M_near(ker_dim[0],n_check*ker_dim[1]);
  743. Matrix<T> T_near(ker_dim[0],n_trg *ker_dim[1]);
  744. #pragma omp parallel for schedule(dynamic)
  745. for(size_t i=0;i<n_check;i++){ // Compute near-interaction for operator M_near
  746. std::vector<T> M_=cheb_integ<T>(0, &check_surf[i*3], 3.0, *this);
  747. for(size_t j=0; j<ker_dim[0]; j++)
  748. for(int k=0; k<ker_dim[1]; k++)
  749. M_near[j][i*ker_dim[1]+k] = M_[j+k*ker_dim[0]];
  750. }
  751. #pragma omp parallel for schedule(dynamic)
  752. for(size_t i=0;i<n_trg;i++){ // Compute near-interaction for targets T_near
  753. std::vector<T> M_=cheb_integ<T>(0, &trg_coord[i*3], 3.0, *this);
  754. for(size_t j=0; j<ker_dim[0]; j++)
  755. for(int k=0; k<ker_dim[1]; k++)
  756. T_near[j][i*ker_dim[1]+k] = M_[j+k*ker_dim[0]];
  757. }
  758. { // M_local = M_analytic - M_near
  759. M_analytic.ReInit(ker_dim[0],n_check*ker_dim[1]); M_analytic.SetZero();
  760. vol_poten(&check_surf[0],n_check,&M_analytic[0][0]);
  761. M_local=M_analytic-M_near;
  762. }
  763. { // T_local = T_analytic - T_near
  764. T_analytic.ReInit(ker_dim[0],n_trg *ker_dim[1]); T_analytic.SetZero();
  765. vol_poten(&trg_coord[0],n_trg,&T_analytic[0][0]);
  766. T_local=T_analytic-T_near;
  767. }
  768. }
  769. Matrix<T> T_err;
  770. { // Now we should be able to compute T_local from M_local
  771. Matrix<T> M_e2c(n_equiv*ker_dim[0],n_check*ker_dim[1]);
  772. BuildMatrix(&equiv_surf[0], n_equiv,
  773. &check_surf[0], n_check, &(M_e2c[0][0]));
  774. Matrix<T> M_e2t(n_equiv*ker_dim[0],n_trg *ker_dim[1]);
  775. BuildMatrix(&equiv_surf[0], n_equiv,
  776. &trg_coord [0], n_trg , &(M_e2t[0][0]));
  777. Matrix<T> M_c2e0, M_c2e1;
  778. {
  779. Matrix<T> U,S,V;
  780. M_e2c.SVD(U,S,V);
  781. T eps=1, max_S=0;
  782. while(eps*(T)0.5+(T)1.0>1.0) eps*=0.5;
  783. for(size_t i=0;i<std::min(S.Dim(0),S.Dim(1));i++){
  784. if(pvfmm::fabs<T>(S[i][i])>max_S) max_S=pvfmm::fabs<T>(S[i][i]);
  785. }
  786. for(size_t i=0;i<S.Dim(0);i++) S[i][i]=(S[i][i]>eps*max_S*4?1.0/S[i][i]:0.0);
  787. M_c2e0=V.Transpose()*S;
  788. M_c2e1=U.Transpose();
  789. }
  790. T_err=(M_local*M_c2e0)*(M_c2e1*M_e2t)-T_local;
  791. }
  792. { // Print relative error
  793. T err_sum=0, analytic_sum=0;
  794. for(size_t i=0;i<T_err .Dim(0)*T_err .Dim(1);i++) err_sum+=pvfmm::fabs<T>(T_err [0][i]);
  795. for(size_t i=0;i<T_analytic.Dim(0)*T_analytic.Dim(1);i++) analytic_sum+=pvfmm::fabs<T>(T_analytic[0][i]);
  796. std::cout<<"Volume Error : "<<err_sum/analytic_sum<<"\n";
  797. }
  798. }
  799. std::cout<<"\n";
  800. }
  801. { // Initialize auxiliary FMM kernels
  802. if(!k_s2m) k_s2m=this;
  803. if(!k_s2l) k_s2l=this;
  804. if(!k_s2t) k_s2t=this;
  805. if(!k_m2m) k_m2m=this;
  806. if(!k_m2l) k_m2l=this;
  807. if(!k_m2t) k_m2t=this;
  808. if(!k_l2l) k_l2l=this;
  809. if(!k_l2t) k_l2t=this;
  810. assert(k_s2t->ker_dim[0]==ker_dim[0]);
  811. assert(k_s2m->ker_dim[0]==k_s2l->ker_dim[0]);
  812. assert(k_s2m->ker_dim[0]==k_s2t->ker_dim[0]);
  813. assert(k_m2m->ker_dim[0]==k_m2l->ker_dim[0]);
  814. assert(k_m2m->ker_dim[0]==k_m2t->ker_dim[0]);
  815. assert(k_l2l->ker_dim[0]==k_l2t->ker_dim[0]);
  816. assert(k_s2t->ker_dim[1]==ker_dim[1]);
  817. assert(k_s2m->ker_dim[1]==k_m2m->ker_dim[1]);
  818. assert(k_s2l->ker_dim[1]==k_l2l->ker_dim[1]);
  819. assert(k_m2l->ker_dim[1]==k_l2l->ker_dim[1]);
  820. assert(k_s2t->ker_dim[1]==k_m2t->ker_dim[1]);
  821. assert(k_s2t->ker_dim[1]==k_l2t->ker_dim[1]);
  822. k_s2m->Initialize(verbose);
  823. k_s2l->Initialize(verbose);
  824. k_s2t->Initialize(verbose);
  825. k_m2m->Initialize(verbose);
  826. k_m2l->Initialize(verbose);
  827. k_m2t->Initialize(verbose);
  828. k_l2l->Initialize(verbose);
  829. k_l2t->Initialize(verbose);
  830. }
  831. }
  832. /**
  833. * \brief Compute the transformation matrix (on the source strength vector)
  834. * to get potential at target coordinates due to sources at the given
  835. * coordinates.
  836. * \param[in] r_src Coordinates of source points.
  837. * \param[in] src_cnt Number of source points.
  838. * \param[in] r_trg Coordinates of target points.
  839. * \param[in] trg_cnt Number of target points.
  840. * \param[out] k_out Output array with potential values.
  841. */
  842. template <class T>
  843. void Kernel<T>::BuildMatrix(T* r_src, int src_cnt,
  844. T* r_trg, int trg_cnt, T* k_out) const{
  845. int dim=3; //Only supporting 3D
  846. memset(k_out, 0, src_cnt*ker_dim[0]*trg_cnt*ker_dim[1]*sizeof(T));
  847. for(int i=0;i<src_cnt;i++) //TODO Optimize this.
  848. for(int j=0;j<ker_dim[0];j++){
  849. std::vector<T> v_src(ker_dim[0],0);
  850. v_src[j]=1.0;
  851. ker_poten(&r_src[i*dim], 1, &v_src[0], 1, r_trg, trg_cnt,
  852. &k_out[(i*ker_dim[0]+j)*trg_cnt*ker_dim[1]], NULL);
  853. }
  854. }
  855. /**
  856. * \brief Generic kernel which rearranges data for vectorization, calls the
  857. * actual uKernel and copies data to the output array in the original order.
  858. */
  859. template <class Real_t, int SRC_DIM, int TRG_DIM, void (*uKernel)(Matrix<Real_t>&, Matrix<Real_t>&, Matrix<Real_t>&, Matrix<Real_t>&)>
  860. void generic_kernel(Real_t* r_src, int src_cnt, Real_t* v_src, int dof, Real_t* r_trg, int trg_cnt, Real_t* v_trg, mem::MemoryManager* mem_mgr){
  861. assert(dof==1);
  862. int VecLen=8;
  863. if(sizeof(Real_t)==sizeof( float)) VecLen=8;
  864. if(sizeof(Real_t)==sizeof(double)) VecLen=4;
  865. #define STACK_BUFF_SIZE 4096
  866. Real_t stack_buff[STACK_BUFF_SIZE+MEM_ALIGN];
  867. Real_t* buff=NULL;
  868. Matrix<Real_t> src_coord;
  869. Matrix<Real_t> src_value;
  870. Matrix<Real_t> trg_coord;
  871. Matrix<Real_t> trg_value;
  872. { // Rearrange data in src_coord, src_coord, trg_coord, trg_value
  873. size_t src_cnt_, trg_cnt_; // counts after zero padding
  874. src_cnt_=((src_cnt+VecLen-1)/VecLen)*VecLen;
  875. trg_cnt_=((trg_cnt+VecLen-1)/VecLen)*VecLen;
  876. size_t buff_size=src_cnt_*(COORD_DIM+SRC_DIM)+
  877. trg_cnt_*(COORD_DIM+TRG_DIM);
  878. if(buff_size>STACK_BUFF_SIZE){ // Allocate buff
  879. buff=mem::aligned_new<Real_t>(buff_size, mem_mgr);
  880. }
  881. Real_t* buff_ptr=buff;
  882. if(!buff_ptr){ // use stack_buff
  883. uintptr_t ptr=(uintptr_t)stack_buff;
  884. static uintptr_t ALIGN_MINUS_ONE=MEM_ALIGN-1;
  885. static uintptr_t NOT_ALIGN_MINUS_ONE=~ALIGN_MINUS_ONE;
  886. ptr=((ptr+ALIGN_MINUS_ONE) & NOT_ALIGN_MINUS_ONE);
  887. buff_ptr=(Real_t*)ptr;
  888. }
  889. src_coord.ReInit(COORD_DIM, src_cnt_,buff_ptr,false); buff_ptr+=COORD_DIM*src_cnt_;
  890. src_value.ReInit( SRC_DIM, src_cnt_,buff_ptr,false); buff_ptr+= SRC_DIM*src_cnt_;
  891. trg_coord.ReInit(COORD_DIM, trg_cnt_,buff_ptr,false); buff_ptr+=COORD_DIM*trg_cnt_;
  892. trg_value.ReInit( TRG_DIM, trg_cnt_,buff_ptr,false);//buff_ptr+= TRG_DIM*trg_cnt_;
  893. { // Set src_coord
  894. size_t i=0;
  895. for( ;i<src_cnt ;i++){
  896. for(size_t j=0;j<COORD_DIM;j++){
  897. src_coord[j][i]=r_src[i*COORD_DIM+j];
  898. }
  899. }
  900. for( ;i<src_cnt_;i++){
  901. for(size_t j=0;j<COORD_DIM;j++){
  902. src_coord[j][i]=0;
  903. }
  904. }
  905. }
  906. { // Set src_value
  907. size_t i=0;
  908. for( ;i<src_cnt ;i++){
  909. for(size_t j=0;j<SRC_DIM;j++){
  910. src_value[j][i]=v_src[i*SRC_DIM+j];
  911. }
  912. }
  913. for( ;i<src_cnt_;i++){
  914. for(size_t j=0;j<SRC_DIM;j++){
  915. src_value[j][i]=0;
  916. }
  917. }
  918. }
  919. { // Set trg_coord
  920. size_t i=0;
  921. for( ;i<trg_cnt ;i++){
  922. for(size_t j=0;j<COORD_DIM;j++){
  923. trg_coord[j][i]=r_trg[i*COORD_DIM+j];
  924. }
  925. }
  926. for( ;i<trg_cnt_;i++){
  927. for(size_t j=0;j<COORD_DIM;j++){
  928. trg_coord[j][i]=0;
  929. }
  930. }
  931. }
  932. { // Set trg_value
  933. size_t i=0;
  934. for( ;i<trg_cnt_;i++){
  935. for(size_t j=0;j<TRG_DIM;j++){
  936. trg_value[j][i]=0;
  937. }
  938. }
  939. }
  940. }
  941. uKernel(src_coord,src_value,trg_coord,trg_value);
  942. { // Set v_trg
  943. for(size_t i=0;i<trg_cnt ;i++){
  944. for(size_t j=0;j<TRG_DIM;j++){
  945. v_trg[i*TRG_DIM+j]+=trg_value[j][i];
  946. }
  947. }
  948. }
  949. if(buff){ // Free memory: buff
  950. mem::aligned_delete<Real_t>(buff);
  951. }
  952. }
  953. ////////////////////////////////////////////////////////////////////////////////
  954. //////// LAPLACE KERNEL ////////
  955. ////////////////////////////////////////////////////////////////////////////////
  956. /**
  957. * \brief Green's function for the Poisson's equation. Kernel tensor
  958. * dimension = 1x1.
  959. */
  960. template <class Real_t, class Vec_t=Real_t, Vec_t (*RSQRT_INTRIN)(Vec_t)=rsqrt_intrin0<Vec_t> >
  961. void laplace_poten_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value, Matrix<Real_t>& trg_coord, Matrix<Real_t>& trg_value){
  962. #define SRC_BLK 1000
  963. size_t VecLen=sizeof(Vec_t)/sizeof(Real_t);
  964. //// Number of newton iterations
  965. size_t NWTN_ITER=0;
  966. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin0<Vec_t,Real_t>) NWTN_ITER=0;
  967. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin1<Vec_t,Real_t>) NWTN_ITER=1;
  968. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin2<Vec_t,Real_t>) NWTN_ITER=2;
  969. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin3<Vec_t,Real_t>) NWTN_ITER=3;
  970. Real_t nwtn_scal=1; // scaling factor for newton iterations
  971. for(int i=0;i<NWTN_ITER;i++){
  972. nwtn_scal=2*nwtn_scal*nwtn_scal*nwtn_scal;
  973. }
  974. const Real_t OOFP = 1.0/(4*nwtn_scal*const_pi<Real_t>());
  975. size_t src_cnt_=src_coord.Dim(1);
  976. size_t trg_cnt_=trg_coord.Dim(1);
  977. for(size_t sblk=0;sblk<src_cnt_;sblk+=SRC_BLK){
  978. size_t src_cnt=src_cnt_-sblk;
  979. if(src_cnt>SRC_BLK) src_cnt=SRC_BLK;
  980. for(size_t t=0;t<trg_cnt_;t+=VecLen){
  981. Vec_t tx=load_intrin<Vec_t>(&trg_coord[0][t]);
  982. Vec_t ty=load_intrin<Vec_t>(&trg_coord[1][t]);
  983. Vec_t tz=load_intrin<Vec_t>(&trg_coord[2][t]);
  984. Vec_t tv=zero_intrin<Vec_t>();
  985. for(size_t s=sblk;s<sblk+src_cnt;s++){
  986. Vec_t dx=sub_intrin(tx,bcast_intrin<Vec_t>(&src_coord[0][s]));
  987. Vec_t dy=sub_intrin(ty,bcast_intrin<Vec_t>(&src_coord[1][s]));
  988. Vec_t dz=sub_intrin(tz,bcast_intrin<Vec_t>(&src_coord[2][s]));
  989. Vec_t sv= bcast_intrin<Vec_t>(&src_value[0][s]) ;
  990. Vec_t r2= mul_intrin(dx,dx) ;
  991. r2=add_intrin(r2,mul_intrin(dy,dy));
  992. r2=add_intrin(r2,mul_intrin(dz,dz));
  993. Vec_t rinv=RSQRT_INTRIN(r2);
  994. tv=add_intrin(tv,mul_intrin(rinv,sv));
  995. }
  996. Vec_t oofp=set_intrin<Vec_t,Real_t>(OOFP);
  997. tv=add_intrin(mul_intrin(tv,oofp),load_intrin<Vec_t>(&trg_value[0][t]));
  998. store_intrin(&trg_value[0][t],tv);
  999. }
  1000. }
  1001. { // Add FLOPS
  1002. #ifndef __MIC__
  1003. Profile::Add_FLOP((long long)trg_cnt_*(long long)src_cnt_*(12+4*(NWTN_ITER)));
  1004. #endif
  1005. }
  1006. #undef SRC_BLK
  1007. }
  1008. template <class T, int newton_iter=0>
  1009. void laplace_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
  1010. #define LAP_KER_NWTN(nwtn) if(newton_iter==nwtn) \
  1011. generic_kernel<Real_t, 1, 1, laplace_poten_uKernel<Real_t,Vec_t, rsqrt_intrin##nwtn<Vec_t,Real_t> > > \
  1012. ((Real_t*)r_src, src_cnt, (Real_t*)v_src, dof, (Real_t*)r_trg, trg_cnt, (Real_t*)v_trg, mem_mgr)
  1013. #define LAPLACE_KERNEL LAP_KER_NWTN(0); LAP_KER_NWTN(1); LAP_KER_NWTN(2); LAP_KER_NWTN(3);
  1014. if(mem::TypeTraits<T>::ID()==mem::TypeTraits<float>::ID()){
  1015. typedef float Real_t;
  1016. #if defined __MIC__
  1017. #define Vec_t Real_t
  1018. #elif defined __AVX__
  1019. #define Vec_t __m256
  1020. #elif defined __SSE3__
  1021. #define Vec_t __m128
  1022. #else
  1023. #define Vec_t Real_t
  1024. #endif
  1025. LAPLACE_KERNEL;
  1026. #undef Vec_t
  1027. }else if(mem::TypeTraits<T>::ID()==mem::TypeTraits<double>::ID()){
  1028. typedef double Real_t;
  1029. #if defined __MIC__
  1030. #define Vec_t Real_t
  1031. #elif defined __AVX__
  1032. #define Vec_t __m256d
  1033. #elif defined __SSE3__
  1034. #define Vec_t __m128d
  1035. #else
  1036. #define Vec_t Real_t
  1037. #endif
  1038. LAPLACE_KERNEL;
  1039. #undef Vec_t
  1040. }else{
  1041. typedef T Real_t;
  1042. #define Vec_t Real_t
  1043. LAPLACE_KERNEL;
  1044. #undef Vec_t
  1045. }
  1046. #undef LAP_KER_NWTN
  1047. #undef LAPLACE_KERNEL
  1048. }
  1049. template <class Real_t>
  1050. void laplace_vol_poten(const Real_t* coord, int n, Real_t* out){
  1051. for(int i=0;i<n;i++){
  1052. const Real_t* c=&coord[i*COORD_DIM];
  1053. Real_t r_2=c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
  1054. out[i]=-r_2/6;
  1055. }
  1056. }
  1057. // Laplace double layer potential.
  1058. template <class Real_t, class Vec_t=Real_t, Vec_t (*RSQRT_INTRIN)(Vec_t)=rsqrt_intrin0<Vec_t> >
  1059. void laplace_dbl_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value, Matrix<Real_t>& trg_coord, Matrix<Real_t>& trg_value){
  1060. #define SRC_BLK 500
  1061. size_t VecLen=sizeof(Vec_t)/sizeof(Real_t);
  1062. //// Number of newton iterations
  1063. size_t NWTN_ITER=0;
  1064. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin0<Vec_t,Real_t>) NWTN_ITER=0;
  1065. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin1<Vec_t,Real_t>) NWTN_ITER=1;
  1066. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin2<Vec_t,Real_t>) NWTN_ITER=2;
  1067. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin3<Vec_t,Real_t>) NWTN_ITER=3;
  1068. Real_t nwtn_scal=1; // scaling factor for newton iterations
  1069. for(int i=0;i<NWTN_ITER;i++){
  1070. nwtn_scal=2*nwtn_scal*nwtn_scal*nwtn_scal;
  1071. }
  1072. const Real_t OOFP = -1.0/(4*nwtn_scal*nwtn_scal*nwtn_scal*const_pi<Real_t>());
  1073. size_t src_cnt_=src_coord.Dim(1);
  1074. size_t trg_cnt_=trg_coord.Dim(1);
  1075. for(size_t sblk=0;sblk<src_cnt_;sblk+=SRC_BLK){
  1076. size_t src_cnt=src_cnt_-sblk;
  1077. if(src_cnt>SRC_BLK) src_cnt=SRC_BLK;
  1078. for(size_t t=0;t<trg_cnt_;t+=VecLen){
  1079. Vec_t tx=load_intrin<Vec_t>(&trg_coord[0][t]);
  1080. Vec_t ty=load_intrin<Vec_t>(&trg_coord[1][t]);
  1081. Vec_t tz=load_intrin<Vec_t>(&trg_coord[2][t]);
  1082. Vec_t tv=zero_intrin<Vec_t>();
  1083. for(size_t s=sblk;s<sblk+src_cnt;s++){
  1084. Vec_t dx=sub_intrin(tx,bcast_intrin<Vec_t>(&src_coord[0][s]));
  1085. Vec_t dy=sub_intrin(ty,bcast_intrin<Vec_t>(&src_coord[1][s]));
  1086. Vec_t dz=sub_intrin(tz,bcast_intrin<Vec_t>(&src_coord[2][s]));
  1087. Vec_t sn0= bcast_intrin<Vec_t>(&src_value[0][s]) ;
  1088. Vec_t sn1= bcast_intrin<Vec_t>(&src_value[1][s]) ;
  1089. Vec_t sn2= bcast_intrin<Vec_t>(&src_value[2][s]) ;
  1090. Vec_t sv= bcast_intrin<Vec_t>(&src_value[3][s]) ;
  1091. Vec_t r2= mul_intrin(dx,dx) ;
  1092. r2=add_intrin(r2,mul_intrin(dy,dy));
  1093. r2=add_intrin(r2,mul_intrin(dz,dz));
  1094. Vec_t rinv=RSQRT_INTRIN(r2);
  1095. Vec_t r3inv=mul_intrin(mul_intrin(rinv,rinv),rinv);
  1096. Vec_t rdotn= mul_intrin(sn0,dx);
  1097. rdotn=add_intrin(rdotn, mul_intrin(sn1,dy));
  1098. rdotn=add_intrin(rdotn, mul_intrin(sn2,dz));
  1099. sv=mul_intrin(sv,rdotn);
  1100. tv=add_intrin(tv,mul_intrin(r3inv,sv));
  1101. }
  1102. Vec_t oofp=set_intrin<Vec_t,Real_t>(OOFP);
  1103. tv=add_intrin(mul_intrin(tv,oofp),load_intrin<Vec_t>(&trg_value[0][t]));
  1104. store_intrin(&trg_value[0][t],tv);
  1105. }
  1106. }
  1107. { // Add FLOPS
  1108. #ifndef __MIC__
  1109. Profile::Add_FLOP((long long)trg_cnt_*(long long)src_cnt_*(20+4*(NWTN_ITER)));
  1110. #endif
  1111. }
  1112. #undef SRC_BLK
  1113. }
  1114. template <class T, int newton_iter=0>
  1115. void laplace_dbl_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
  1116. #define LAP_KER_NWTN(nwtn) if(newton_iter==nwtn) \
  1117. generic_kernel<Real_t, 4, 1, laplace_dbl_uKernel<Real_t,Vec_t, rsqrt_intrin##nwtn<Vec_t,Real_t> > > \
  1118. ((Real_t*)r_src, src_cnt, (Real_t*)v_src, dof, (Real_t*)r_trg, trg_cnt, (Real_t*)v_trg, mem_mgr)
  1119. #define LAPLACE_KERNEL LAP_KER_NWTN(0); LAP_KER_NWTN(1); LAP_KER_NWTN(2); LAP_KER_NWTN(3);
  1120. if(mem::TypeTraits<T>::ID()==mem::TypeTraits<float>::ID()){
  1121. typedef float Real_t;
  1122. #if defined __MIC__
  1123. #define Vec_t Real_t
  1124. #elif defined __AVX__
  1125. #define Vec_t __m256
  1126. #elif defined __SSE3__
  1127. #define Vec_t __m128
  1128. #else
  1129. #define Vec_t Real_t
  1130. #endif
  1131. LAPLACE_KERNEL;
  1132. #undef Vec_t
  1133. }else if(mem::TypeTraits<T>::ID()==mem::TypeTraits<double>::ID()){
  1134. typedef double Real_t;
  1135. #if defined __MIC__
  1136. #define Vec_t Real_t
  1137. #elif defined __AVX__
  1138. #define Vec_t __m256d
  1139. #elif defined __SSE3__
  1140. #define Vec_t __m128d
  1141. #else
  1142. #define Vec_t Real_t
  1143. #endif
  1144. LAPLACE_KERNEL;
  1145. #undef Vec_t
  1146. }else{
  1147. typedef T Real_t;
  1148. #define Vec_t Real_t
  1149. LAPLACE_KERNEL;
  1150. #undef Vec_t
  1151. }
  1152. #undef LAP_KER_NWTN
  1153. #undef LAPLACE_KERNEL
  1154. }
  1155. // Laplace grdient kernel.
  1156. template <class Real_t, class Vec_t=Real_t, Vec_t (*RSQRT_INTRIN)(Vec_t)=rsqrt_intrin0<Vec_t> >
  1157. void laplace_grad_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value, Matrix<Real_t>& trg_coord, Matrix<Real_t>& trg_value){
  1158. #define SRC_BLK 500
  1159. size_t VecLen=sizeof(Vec_t)/sizeof(Real_t);
  1160. //// Number of newton iterations
  1161. size_t NWTN_ITER=0;
  1162. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin0<Vec_t,Real_t>) NWTN_ITER=0;
  1163. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin1<Vec_t,Real_t>) NWTN_ITER=1;
  1164. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin2<Vec_t,Real_t>) NWTN_ITER=2;
  1165. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin3<Vec_t,Real_t>) NWTN_ITER=3;
  1166. Real_t nwtn_scal=1; // scaling factor for newton iterations
  1167. for(int i=0;i<NWTN_ITER;i++){
  1168. nwtn_scal=2*nwtn_scal*nwtn_scal*nwtn_scal;
  1169. }
  1170. const Real_t OOFP = -1.0/(4*nwtn_scal*nwtn_scal*nwtn_scal*const_pi<Real_t>());
  1171. size_t src_cnt_=src_coord.Dim(1);
  1172. size_t trg_cnt_=trg_coord.Dim(1);
  1173. for(size_t sblk=0;sblk<src_cnt_;sblk+=SRC_BLK){
  1174. size_t src_cnt=src_cnt_-sblk;
  1175. if(src_cnt>SRC_BLK) src_cnt=SRC_BLK;
  1176. for(size_t t=0;t<trg_cnt_;t+=VecLen){
  1177. Vec_t tx=load_intrin<Vec_t>(&trg_coord[0][t]);
  1178. Vec_t ty=load_intrin<Vec_t>(&trg_coord[1][t]);
  1179. Vec_t tz=load_intrin<Vec_t>(&trg_coord[2][t]);
  1180. Vec_t tv0=zero_intrin<Vec_t>();
  1181. Vec_t tv1=zero_intrin<Vec_t>();
  1182. Vec_t tv2=zero_intrin<Vec_t>();
  1183. for(size_t s=sblk;s<sblk+src_cnt;s++){
  1184. Vec_t dx=sub_intrin(tx,bcast_intrin<Vec_t>(&src_coord[0][s]));
  1185. Vec_t dy=sub_intrin(ty,bcast_intrin<Vec_t>(&src_coord[1][s]));
  1186. Vec_t dz=sub_intrin(tz,bcast_intrin<Vec_t>(&src_coord[2][s]));
  1187. Vec_t sv= bcast_intrin<Vec_t>(&src_value[0][s]) ;
  1188. Vec_t r2= mul_intrin(dx,dx) ;
  1189. r2=add_intrin(r2,mul_intrin(dy,dy));
  1190. r2=add_intrin(r2,mul_intrin(dz,dz));
  1191. Vec_t rinv=RSQRT_INTRIN(r2);
  1192. Vec_t r3inv=mul_intrin(mul_intrin(rinv,rinv),rinv);
  1193. sv=mul_intrin(sv,r3inv);
  1194. tv0=add_intrin(tv0,mul_intrin(sv,dx));
  1195. tv1=add_intrin(tv1,mul_intrin(sv,dy));
  1196. tv2=add_intrin(tv2,mul_intrin(sv,dz));
  1197. }
  1198. Vec_t oofp=set_intrin<Vec_t,Real_t>(OOFP);
  1199. tv0=add_intrin(mul_intrin(tv0,oofp),load_intrin<Vec_t>(&trg_value[0][t]));
  1200. tv1=add_intrin(mul_intrin(tv1,oofp),load_intrin<Vec_t>(&trg_value[1][t]));
  1201. tv2=add_intrin(mul_intrin(tv2,oofp),load_intrin<Vec_t>(&trg_value[2][t]));
  1202. store_intrin(&trg_value[0][t],tv0);
  1203. store_intrin(&trg_value[1][t],tv1);
  1204. store_intrin(&trg_value[2][t],tv2);
  1205. }
  1206. }
  1207. { // Add FLOPS
  1208. #ifndef __MIC__
  1209. Profile::Add_FLOP((long long)trg_cnt_*(long long)src_cnt_*(19+4*(NWTN_ITER)));
  1210. #endif
  1211. }
  1212. #undef SRC_BLK
  1213. }
  1214. template <class T, int newton_iter=0>
  1215. void laplace_grad(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
  1216. #define LAP_KER_NWTN(nwtn) if(newton_iter==nwtn) \
  1217. generic_kernel<Real_t, 1, 3, laplace_grad_uKernel<Real_t,Vec_t, rsqrt_intrin##nwtn<Vec_t,Real_t> > > \
  1218. ((Real_t*)r_src, src_cnt, (Real_t*)v_src, dof, (Real_t*)r_trg, trg_cnt, (Real_t*)v_trg, mem_mgr)
  1219. #define LAPLACE_KERNEL LAP_KER_NWTN(0); LAP_KER_NWTN(1); LAP_KER_NWTN(2); LAP_KER_NWTN(3);
  1220. if(mem::TypeTraits<T>::ID()==mem::TypeTraits<float>::ID()){
  1221. typedef float Real_t;
  1222. #if defined __MIC__
  1223. #define Vec_t Real_t
  1224. #elif defined __AVX__
  1225. #define Vec_t __m256
  1226. #elif defined __SSE3__
  1227. #define Vec_t __m128
  1228. #else
  1229. #define Vec_t Real_t
  1230. #endif
  1231. LAPLACE_KERNEL;
  1232. #undef Vec_t
  1233. }else if(mem::TypeTraits<T>::ID()==mem::TypeTraits<double>::ID()){
  1234. typedef double Real_t;
  1235. #if defined __MIC__
  1236. #define Vec_t Real_t
  1237. #elif defined __AVX__
  1238. #define Vec_t __m256d
  1239. #elif defined __SSE3__
  1240. #define Vec_t __m128d
  1241. #else
  1242. #define Vec_t Real_t
  1243. #endif
  1244. LAPLACE_KERNEL;
  1245. #undef Vec_t
  1246. }else{
  1247. typedef T Real_t;
  1248. #define Vec_t Real_t
  1249. LAPLACE_KERNEL;
  1250. #undef Vec_t
  1251. }
  1252. #undef LAP_KER_NWTN
  1253. #undef LAPLACE_KERNEL
  1254. }
  1255. template<class T> const Kernel<T>& LaplaceKernel<T>::potential(){
  1256. static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,1>, laplace_dbl_poten<T,1> >("laplace" , 3, std::pair<int,int>(1,1),
  1257. NULL,NULL,NULL, NULL,NULL,NULL, NULL,NULL, &laplace_vol_poten<T>);
  1258. return potn_ker;
  1259. }
  1260. template<class T> const Kernel<T>& LaplaceKernel<T>::gradient(){
  1261. static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,1>, laplace_dbl_poten<T,1> >("laplace" , 3, std::pair<int,int>(1,1));
  1262. static Kernel<T> grad_ker=BuildKernel<T, laplace_grad <T,1> >("laplace_grad", 3, std::pair<int,int>(1,3),
  1263. &potn_ker, &potn_ker, NULL, &potn_ker, &potn_ker, NULL, &potn_ker, NULL);
  1264. return grad_ker;
  1265. }
  1266. template<> inline const Kernel<double>& LaplaceKernel<double>::potential(){
  1267. typedef double T;
  1268. static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,2>, laplace_dbl_poten<T,2> >("laplace" , 3, std::pair<int,int>(1,1),
  1269. NULL,NULL,NULL, NULL,NULL,NULL, NULL,NULL, &laplace_vol_poten<double>);
  1270. return potn_ker;
  1271. }
  1272. template<> inline const Kernel<double>& LaplaceKernel<double>::gradient(){
  1273. typedef double T;
  1274. static Kernel<T> potn_ker=BuildKernel<T, laplace_poten<T,2>, laplace_dbl_poten<T,2> >("laplace" , 3, std::pair<int,int>(1,1));
  1275. static Kernel<T> grad_ker=BuildKernel<T, laplace_grad <T,2> >("laplace_grad", 3, std::pair<int,int>(1,3),
  1276. &potn_ker, &potn_ker, NULL, &potn_ker, &potn_ker, NULL, &potn_ker, NULL);
  1277. return grad_ker;
  1278. }
  1279. ////////////////////////////////////////////////////////////////////////////////
  1280. //////// STOKES KERNEL ////////
  1281. ////////////////////////////////////////////////////////////////////////////////
  1282. /**
  1283. * \brief Green's function for the Stokes's equation. Kernel tensor
  1284. * dimension = 3x3.
  1285. */
  1286. template <class Real_t, class Vec_t=Real_t, Vec_t (*RSQRT_INTRIN)(Vec_t)=rsqrt_intrin0<Vec_t> >
  1287. void stokes_vel_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value, Matrix<Real_t>& trg_coord, Matrix<Real_t>& trg_value){
  1288. #define SRC_BLK 500
  1289. size_t VecLen=sizeof(Vec_t)/sizeof(Real_t);
  1290. //// Number of newton iterations
  1291. size_t NWTN_ITER=0;
  1292. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin0<Vec_t,Real_t>) NWTN_ITER=0;
  1293. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin1<Vec_t,Real_t>) NWTN_ITER=1;
  1294. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin2<Vec_t,Real_t>) NWTN_ITER=2;
  1295. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin3<Vec_t,Real_t>) NWTN_ITER=3;
  1296. Real_t nwtn_scal=1; // scaling factor for newton iterations
  1297. for(int i=0;i<NWTN_ITER;i++){
  1298. nwtn_scal=2*nwtn_scal*nwtn_scal*nwtn_scal;
  1299. }
  1300. const Real_t OOEP = 1.0/(8*nwtn_scal*const_pi<Real_t>());
  1301. Vec_t inv_nwtn_scal2=set_intrin<Vec_t,Real_t>(1.0/(nwtn_scal*nwtn_scal));
  1302. size_t src_cnt_=src_coord.Dim(1);
  1303. size_t trg_cnt_=trg_coord.Dim(1);
  1304. for(size_t sblk=0;sblk<src_cnt_;sblk+=SRC_BLK){
  1305. size_t src_cnt=src_cnt_-sblk;
  1306. if(src_cnt>SRC_BLK) src_cnt=SRC_BLK;
  1307. for(size_t t=0;t<trg_cnt_;t+=VecLen){
  1308. Vec_t tx=load_intrin<Vec_t>(&trg_coord[0][t]);
  1309. Vec_t ty=load_intrin<Vec_t>(&trg_coord[1][t]);
  1310. Vec_t tz=load_intrin<Vec_t>(&trg_coord[2][t]);
  1311. Vec_t tvx=zero_intrin<Vec_t>();
  1312. Vec_t tvy=zero_intrin<Vec_t>();
  1313. Vec_t tvz=zero_intrin<Vec_t>();
  1314. for(size_t s=sblk;s<sblk+src_cnt;s++){
  1315. Vec_t dx=sub_intrin(tx,bcast_intrin<Vec_t>(&src_coord[0][s]));
  1316. Vec_t dy=sub_intrin(ty,bcast_intrin<Vec_t>(&src_coord[1][s]));
  1317. Vec_t dz=sub_intrin(tz,bcast_intrin<Vec_t>(&src_coord[2][s]));
  1318. Vec_t svx= bcast_intrin<Vec_t>(&src_value[0][s]) ;
  1319. Vec_t svy= bcast_intrin<Vec_t>(&src_value[1][s]) ;
  1320. Vec_t svz= bcast_intrin<Vec_t>(&src_value[2][s]) ;
  1321. Vec_t r2= mul_intrin(dx,dx) ;
  1322. r2=add_intrin(r2,mul_intrin(dy,dy));
  1323. r2=add_intrin(r2,mul_intrin(dz,dz));
  1324. Vec_t rinv=RSQRT_INTRIN(r2);
  1325. Vec_t rinv2=mul_intrin(mul_intrin(rinv,rinv),inv_nwtn_scal2);
  1326. Vec_t inner_prod= mul_intrin(svx,dx) ;
  1327. inner_prod=add_intrin(inner_prod,mul_intrin(svy,dy));
  1328. inner_prod=add_intrin(inner_prod,mul_intrin(svz,dz));
  1329. inner_prod=mul_intrin(inner_prod,rinv2);
  1330. tvx=add_intrin(tvx,mul_intrin(rinv,add_intrin(svx,mul_intrin(dx,inner_prod))));
  1331. tvy=add_intrin(tvy,mul_intrin(rinv,add_intrin(svy,mul_intrin(dy,inner_prod))));
  1332. tvz=add_intrin(tvz,mul_intrin(rinv,add_intrin(svz,mul_intrin(dz,inner_prod))));
  1333. }
  1334. Vec_t ooep=set_intrin<Vec_t,Real_t>(OOEP);
  1335. tvx=add_intrin(mul_intrin(tvx,ooep),load_intrin<Vec_t>(&trg_value[0][t]));
  1336. tvy=add_intrin(mul_intrin(tvy,ooep),load_intrin<Vec_t>(&trg_value[1][t]));
  1337. tvz=add_intrin(mul_intrin(tvz,ooep),load_intrin<Vec_t>(&trg_value[2][t]));
  1338. store_intrin(&trg_value[0][t],tvx);
  1339. store_intrin(&trg_value[1][t],tvy);
  1340. store_intrin(&trg_value[2][t],tvz);
  1341. }
  1342. }
  1343. { // Add FLOPS
  1344. #ifndef __MIC__
  1345. Profile::Add_FLOP((long long)trg_cnt_*(long long)src_cnt_*(29+4*(NWTN_ITER)));
  1346. #endif
  1347. }
  1348. #undef SRC_BLK
  1349. }
  1350. template <class T, int newton_iter=0>
  1351. void stokes_vel(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
  1352. #define STK_KER_NWTN(nwtn) if(newton_iter==nwtn) \
  1353. generic_kernel<Real_t, 3, 3, stokes_vel_uKernel<Real_t,Vec_t, rsqrt_intrin##nwtn<Vec_t,Real_t> > > \
  1354. ((Real_t*)r_src, src_cnt, (Real_t*)v_src, dof, (Real_t*)r_trg, trg_cnt, (Real_t*)v_trg, mem_mgr)
  1355. #define STOKES_KERNEL STK_KER_NWTN(0); STK_KER_NWTN(1); STK_KER_NWTN(2); STK_KER_NWTN(3);
  1356. if(mem::TypeTraits<T>::ID()==mem::TypeTraits<float>::ID()){
  1357. typedef float Real_t;
  1358. #if defined __MIC__
  1359. #define Vec_t Real_t
  1360. #elif defined __AVX__
  1361. #define Vec_t __m256
  1362. #elif defined __SSE3__
  1363. #define Vec_t __m128
  1364. #else
  1365. #define Vec_t Real_t
  1366. #endif
  1367. STOKES_KERNEL;
  1368. #undef Vec_t
  1369. }else if(mem::TypeTraits<T>::ID()==mem::TypeTraits<double>::ID()){
  1370. typedef double Real_t;
  1371. #if defined __MIC__
  1372. #define Vec_t Real_t
  1373. #elif defined __AVX__
  1374. #define Vec_t __m256d
  1375. #elif defined __SSE3__
  1376. #define Vec_t __m128d
  1377. #else
  1378. #define Vec_t Real_t
  1379. #endif
  1380. STOKES_KERNEL;
  1381. #undef Vec_t
  1382. }else{
  1383. typedef T Real_t;
  1384. #define Vec_t Real_t
  1385. STOKES_KERNEL;
  1386. #undef Vec_t
  1387. }
  1388. #undef STK_KER_NWTN
  1389. #undef STOKES_KERNEL
  1390. }
  1391. template <class Real_t>
  1392. void stokes_vol_poten(const Real_t* coord, int n, Real_t* out){
  1393. for(int i=0;i<n;i++){
  1394. const Real_t* c=&coord[i*COORD_DIM];
  1395. Real_t rx_2=c[1]*c[1]+c[2]*c[2];
  1396. Real_t ry_2=c[0]*c[0]+c[2]*c[2];
  1397. Real_t rz_2=c[0]*c[0]+c[1]*c[1];
  1398. out[n*3*0+i*3+0]=-rx_2/6; out[n*3*0+i*3+1]= 0; out[n*3*0+i*3+2]= 0;
  1399. out[n*3*1+i*3+0]= 0; out[n*3*1+i*3+1]=-ry_2/6; out[n*3*1+i*3+2]= 0;
  1400. out[n*3*2+i*3+0]= 0; out[n*3*2+i*3+1]= 0; out[n*3*2+i*3+2]=-rz_2/6;
  1401. }
  1402. }
  1403. template <class T>
  1404. void stokes_sym_dip(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr){
  1405. #ifndef __MIC__
  1406. Profile::Add_FLOP((long long)trg_cnt*(long long)src_cnt*(47*dof));
  1407. #endif
  1408. const T mu=1.0;
  1409. const T OOEPMU = -1.0/(8.0*const_pi<T>()*mu);
  1410. for(int t=0;t<trg_cnt;t++){
  1411. for(int i=0;i<dof;i++){
  1412. T p[3]={0,0,0};
  1413. for(int s=0;s<src_cnt;s++){
  1414. T dR[3]={r_trg[3*t ]-r_src[3*s ],
  1415. r_trg[3*t+1]-r_src[3*s+1],
  1416. r_trg[3*t+2]-r_src[3*s+2]};
  1417. T R = (dR[0]*dR[0]+dR[1]*dR[1]+dR[2]*dR[2]);
  1418. if (R!=0){
  1419. T invR2=1.0/R;
  1420. T invR=pvfmm::sqrt<T>(invR2);
  1421. T invR3=invR2*invR;
  1422. T* f=&v_src[(s*dof+i)*6+0];
  1423. T* n=&v_src[(s*dof+i)*6+3];
  1424. T r_dot_n=(n[0]*dR[0]+n[1]*dR[1]+n[2]*dR[2]);
  1425. T r_dot_f=(f[0]*dR[0]+f[1]*dR[1]+f[2]*dR[2]);
  1426. T n_dot_f=(f[0]* n[0]+f[1]* n[1]+f[2]* n[2]);
  1427. p[0] += dR[0]*(n_dot_f - 3*r_dot_n*r_dot_f*invR2)*invR3;
  1428. p[1] += dR[1]*(n_dot_f - 3*r_dot_n*r_dot_f*invR2)*invR3;
  1429. p[2] += dR[2]*(n_dot_f - 3*r_dot_n*r_dot_f*invR2)*invR3;
  1430. }
  1431. }
  1432. k_out[(t*dof+i)*3+0] += p[0]*OOEPMU;
  1433. k_out[(t*dof+i)*3+1] += p[1]*OOEPMU;
  1434. k_out[(t*dof+i)*3+2] += p[2]*OOEPMU;
  1435. }
  1436. }
  1437. }
  1438. template <class T>
  1439. void stokes_press(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr){
  1440. #ifndef __MIC__
  1441. Profile::Add_FLOP((long long)trg_cnt*(long long)src_cnt*(17*dof));
  1442. #endif
  1443. const T OOFP = 1.0/(4.0*const_pi<T>());
  1444. for(int t=0;t<trg_cnt;t++){
  1445. for(int i=0;i<dof;i++){
  1446. T p=0;
  1447. for(int s=0;s<src_cnt;s++){
  1448. T dR[3]={r_trg[3*t ]-r_src[3*s ],
  1449. r_trg[3*t+1]-r_src[3*s+1],
  1450. r_trg[3*t+2]-r_src[3*s+2]};
  1451. T R = (dR[0]*dR[0]+dR[1]*dR[1]+dR[2]*dR[2]);
  1452. if (R!=0){
  1453. T invR2=1.0/R;
  1454. T invR=pvfmm::sqrt<T>(invR2);
  1455. T invR3=invR2*invR;
  1456. T v_src[3]={v_src_[(s*dof+i)*3 ],
  1457. v_src_[(s*dof+i)*3+1],
  1458. v_src_[(s*dof+i)*3+2]};
  1459. T inner_prod=(v_src[0]*dR[0] +
  1460. v_src[1]*dR[1] +
  1461. v_src[2]*dR[2])* invR3;
  1462. p += inner_prod;
  1463. }
  1464. }
  1465. k_out[t*dof+i] += p*OOFP;
  1466. }
  1467. }
  1468. }
  1469. template <class T>
  1470. void stokes_stress(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr){
  1471. #ifndef __MIC__
  1472. Profile::Add_FLOP((long long)trg_cnt*(long long)src_cnt*(45*dof));
  1473. #endif
  1474. const T TOFP = -3.0/(4.0*const_pi<T>());
  1475. for(int t=0;t<trg_cnt;t++){
  1476. for(int i=0;i<dof;i++){
  1477. T p[9]={0,0,0,
  1478. 0,0,0,
  1479. 0,0,0};
  1480. for(int s=0;s<src_cnt;s++){
  1481. T dR[3]={r_trg[3*t ]-r_src[3*s ],
  1482. r_trg[3*t+1]-r_src[3*s+1],
  1483. r_trg[3*t+2]-r_src[3*s+2]};
  1484. T R = (dR[0]*dR[0]+dR[1]*dR[1]+dR[2]*dR[2]);
  1485. if (R!=0){
  1486. T invR2=1.0/R;
  1487. T invR=pvfmm::sqrt<T>(invR2);
  1488. T invR3=invR2*invR;
  1489. T invR5=invR3*invR2;
  1490. T v_src[3]={v_src_[(s*dof+i)*3 ],
  1491. v_src_[(s*dof+i)*3+1],
  1492. v_src_[(s*dof+i)*3+2]};
  1493. T inner_prod=(v_src[0]*dR[0] +
  1494. v_src[1]*dR[1] +
  1495. v_src[2]*dR[2])* invR5;
  1496. p[0] += inner_prod*dR[0]*dR[0]; p[1] += inner_prod*dR[1]*dR[0]; p[2] += inner_prod*dR[2]*dR[0];
  1497. p[3] += inner_prod*dR[0]*dR[1]; p[4] += inner_prod*dR[1]*dR[1]; p[5] += inner_prod*dR[2]*dR[1];
  1498. p[6] += inner_prod*dR[0]*dR[2]; p[7] += inner_prod*dR[1]*dR[2]; p[8] += inner_prod*dR[2]*dR[2];
  1499. }
  1500. }
  1501. k_out[(t*dof+i)*9+0] += p[0]*TOFP;
  1502. k_out[(t*dof+i)*9+1] += p[1]*TOFP;
  1503. k_out[(t*dof+i)*9+2] += p[2]*TOFP;
  1504. k_out[(t*dof+i)*9+3] += p[3]*TOFP;
  1505. k_out[(t*dof+i)*9+4] += p[4]*TOFP;
  1506. k_out[(t*dof+i)*9+5] += p[5]*TOFP;
  1507. k_out[(t*dof+i)*9+6] += p[6]*TOFP;
  1508. k_out[(t*dof+i)*9+7] += p[7]*TOFP;
  1509. k_out[(t*dof+i)*9+8] += p[8]*TOFP;
  1510. }
  1511. }
  1512. }
  1513. template <class T>
  1514. void stokes_grad(T* r_src, int src_cnt, T* v_src_, int dof, T* r_trg, int trg_cnt, T* k_out, mem::MemoryManager* mem_mgr){
  1515. #ifndef __MIC__
  1516. Profile::Add_FLOP((long long)trg_cnt*(long long)src_cnt*(89*dof));
  1517. #endif
  1518. const T mu=1.0;
  1519. const T OOEPMU = 1.0/(8.0*const_pi<T>()*mu);
  1520. for(int t=0;t<trg_cnt;t++){
  1521. for(int i=0;i<dof;i++){
  1522. T p[9]={0,0,0,
  1523. 0,0,0,
  1524. 0,0,0};
  1525. for(int s=0;s<src_cnt;s++){
  1526. T dR[3]={r_trg[3*t ]-r_src[3*s ],
  1527. r_trg[3*t+1]-r_src[3*s+1],
  1528. r_trg[3*t+2]-r_src[3*s+2]};
  1529. T R = (dR[0]*dR[0]+dR[1]*dR[1]+dR[2]*dR[2]);
  1530. if (R!=0){
  1531. T invR2=1.0/R;
  1532. T invR=pvfmm::sqrt<T>(invR2);
  1533. T invR3=invR2*invR;
  1534. T v_src[3]={v_src_[(s*dof+i)*3 ],
  1535. v_src_[(s*dof+i)*3+1],
  1536. v_src_[(s*dof+i)*3+2]};
  1537. T inner_prod=(v_src[0]*dR[0] +
  1538. v_src[1]*dR[1] +
  1539. v_src[2]*dR[2]);
  1540. p[0] += ( inner_prod*(1-3*dR[0]*dR[0]*invR2))*invR3; //6
  1541. p[1] += (dR[1]*v_src[0]-v_src[1]*dR[0]+inner_prod*( -3*dR[1]*dR[0]*invR2))*invR3; //9
  1542. p[2] += (dR[2]*v_src[0]-v_src[2]*dR[0]+inner_prod*( -3*dR[2]*dR[0]*invR2))*invR3;
  1543. p[3] += (dR[0]*v_src[1]-v_src[0]*dR[1]+inner_prod*( -3*dR[0]*dR[1]*invR2))*invR3;
  1544. p[4] += ( inner_prod*(1-3*dR[1]*dR[1]*invR2))*invR3;
  1545. p[5] += (dR[2]*v_src[1]-v_src[2]*dR[1]+inner_prod*( -3*dR[2]*dR[1]*invR2))*invR3;
  1546. p[6] += (dR[0]*v_src[2]-v_src[0]*dR[2]+inner_prod*( -3*dR[0]*dR[2]*invR2))*invR3;
  1547. p[7] += (dR[1]*v_src[2]-v_src[1]*dR[2]+inner_prod*( -3*dR[1]*dR[2]*invR2))*invR3;
  1548. p[8] += ( inner_prod*(1-3*dR[2]*dR[2]*invR2))*invR3;
  1549. }
  1550. }
  1551. k_out[(t*dof+i)*9+0] += p[0]*OOEPMU;
  1552. k_out[(t*dof+i)*9+1] += p[1]*OOEPMU;
  1553. k_out[(t*dof+i)*9+2] += p[2]*OOEPMU;
  1554. k_out[(t*dof+i)*9+3] += p[3]*OOEPMU;
  1555. k_out[(t*dof+i)*9+4] += p[4]*OOEPMU;
  1556. k_out[(t*dof+i)*9+5] += p[5]*OOEPMU;
  1557. k_out[(t*dof+i)*9+6] += p[6]*OOEPMU;
  1558. k_out[(t*dof+i)*9+7] += p[7]*OOEPMU;
  1559. k_out[(t*dof+i)*9+8] += p[8]*OOEPMU;
  1560. }
  1561. }
  1562. }
  1563. #ifndef __MIC__
  1564. #ifdef USE_SSE
  1565. namespace
  1566. {
  1567. #define IDEAL_ALIGNMENT 16
  1568. #define SIMD_LEN (int)(IDEAL_ALIGNMENT / sizeof(double))
  1569. #define DECL_SIMD_ALIGNED __declspec(align(IDEAL_ALIGNMENT))
  1570. void stokesPressureSSE(
  1571. const int ns,
  1572. const int nt,
  1573. const double *sx,
  1574. const double *sy,
  1575. const double *sz,
  1576. const double *tx,
  1577. const double *ty,
  1578. const double *tz,
  1579. const double *srcDen,
  1580. double *trgVal)
  1581. {
  1582. if ( size_t(sx)%IDEAL_ALIGNMENT || size_t(sy)%IDEAL_ALIGNMENT || size_t(sz)%IDEAL_ALIGNMENT )
  1583. abort();
  1584. double OOFP = 1.0/(4.0*const_pi<double>());
  1585. __m128d temp_press;
  1586. double aux_arr[SIMD_LEN+1];
  1587. double *tempval_press;
  1588. if (size_t(aux_arr)%IDEAL_ALIGNMENT) // if aux_arr is misaligned
  1589. {
  1590. tempval_press = aux_arr + 1;
  1591. if (size_t(tempval_press)%IDEAL_ALIGNMENT)
  1592. abort();
  1593. }
  1594. else
  1595. tempval_press = aux_arr;
  1596. /*! One over eight pi */
  1597. __m128d oofp = _mm_set1_pd (OOFP);
  1598. __m128d half = _mm_set1_pd (0.5);
  1599. __m128d opf = _mm_set1_pd (1.5);
  1600. __m128d zero = _mm_setzero_pd ();
  1601. // loop over sources
  1602. int i = 0;
  1603. for (; i < nt; i++) {
  1604. temp_press = _mm_setzero_pd();
  1605. __m128d txi = _mm_load1_pd (&tx[i]);
  1606. __m128d tyi = _mm_load1_pd (&ty[i]);
  1607. __m128d tzi = _mm_load1_pd (&tz[i]);
  1608. int j = 0;
  1609. // Load and calculate in groups of SIMD_LEN
  1610. for (; j + SIMD_LEN <= ns; j+=SIMD_LEN) {
  1611. __m128d sxj = _mm_load_pd (&sx[j]);
  1612. __m128d syj = _mm_load_pd (&sy[j]);
  1613. __m128d szj = _mm_load_pd (&sz[j]);
  1614. __m128d sdenx = _mm_set_pd (srcDen[(j+1)*3], srcDen[j*3]);
  1615. __m128d sdeny = _mm_set_pd (srcDen[(j+1)*3+1], srcDen[j*3+1]);
  1616. __m128d sdenz = _mm_set_pd (srcDen[(j+1)*3+2], srcDen[j*3+2]);
  1617. __m128d dX, dY, dZ;
  1618. __m128d dR2;
  1619. __m128d S;
  1620. dX = _mm_sub_pd(txi , sxj);
  1621. dY = _mm_sub_pd(tyi , syj);
  1622. dZ = _mm_sub_pd(tzi , szj);
  1623. sxj = _mm_mul_pd(dX, dX);
  1624. syj = _mm_mul_pd(dY, dY);
  1625. szj = _mm_mul_pd(dZ, dZ);
  1626. dR2 = _mm_add_pd(sxj, syj);
  1627. dR2 = _mm_add_pd(szj, dR2);
  1628. __m128d temp = _mm_cmpeq_pd (dR2, zero);
  1629. __m128d xhalf = _mm_mul_pd (half, dR2);
  1630. __m128 dR2_s = _mm_cvtpd_ps(dR2);
  1631. __m128 S_s = _mm_rsqrt_ps(dR2_s);
  1632. __m128d S_d = _mm_cvtps_pd(S_s);
  1633. // To handle the condition when src and trg coincide
  1634. S_d = _mm_andnot_pd (temp, S_d);
  1635. S = _mm_mul_pd (S_d, S_d);
  1636. S = _mm_mul_pd (S, xhalf);
  1637. S = _mm_sub_pd (opf, S);
  1638. S = _mm_mul_pd (S, S_d);
  1639. __m128d dotx = _mm_mul_pd (dX, sdenx);
  1640. __m128d doty = _mm_mul_pd (dY, sdeny);
  1641. __m128d dotz = _mm_mul_pd (dZ, sdenz);
  1642. __m128d dot_sum = _mm_add_pd (dotx, doty);
  1643. dot_sum = _mm_add_pd (dot_sum, dotz);
  1644. dot_sum = _mm_mul_pd (dot_sum, S);
  1645. dot_sum = _mm_mul_pd (dot_sum, S);
  1646. dot_sum = _mm_mul_pd (dot_sum, S);
  1647. temp_press = _mm_add_pd (dot_sum, temp_press);
  1648. }
  1649. temp_press = _mm_mul_pd (temp_press, oofp);
  1650. _mm_store_pd(tempval_press, temp_press);
  1651. for (int k = 0; k < SIMD_LEN; k++) {
  1652. trgVal[i] += tempval_press[k];
  1653. }
  1654. for (; j < ns; j++) {
  1655. double x = tx[i] - sx[j];
  1656. double y = ty[i] - sy[j];
  1657. double z = tz[i] - sz[j];
  1658. double r2 = x*x + y*y + z*z;
  1659. double r = pvfmm::sqrt<T>(r2);
  1660. double invdr;
  1661. if (r == 0)
  1662. invdr = 0;
  1663. else
  1664. invdr = 1/r;
  1665. double dot = (x*srcDen[j*3] + y*srcDen[j*3+1] + z*srcDen[j*3+2]) * invdr * invdr * invdr;
  1666. trgVal[i] += dot*OOFP;
  1667. }
  1668. }
  1669. return;
  1670. }
  1671. void stokesStressSSE(
  1672. const int ns,
  1673. const int nt,
  1674. const double *sx,
  1675. const double *sy,
  1676. const double *sz,
  1677. const double *tx,
  1678. const double *ty,
  1679. const double *tz,
  1680. const double *srcDen,
  1681. double *trgVal)
  1682. {
  1683. if ( size_t(sx)%IDEAL_ALIGNMENT || size_t(sy)%IDEAL_ALIGNMENT || size_t(sz)%IDEAL_ALIGNMENT )
  1684. abort();
  1685. double TOFP = -3.0/(4.0*const_pi<double>());
  1686. __m128d tempxx; __m128d tempxy; __m128d tempxz;
  1687. __m128d tempyx; __m128d tempyy; __m128d tempyz;
  1688. __m128d tempzx; __m128d tempzy; __m128d tempzz;
  1689. double aux_arr[9*SIMD_LEN+1];
  1690. double *tempvalxx, *tempvalxy, *tempvalxz;
  1691. double *tempvalyx, *tempvalyy, *tempvalyz;
  1692. double *tempvalzx, *tempvalzy, *tempvalzz;
  1693. if (size_t(aux_arr)%IDEAL_ALIGNMENT) // if aux_arr is misaligned
  1694. {
  1695. tempvalxx = aux_arr + 1;
  1696. if (size_t(tempvalxx)%IDEAL_ALIGNMENT)
  1697. abort();
  1698. }
  1699. else
  1700. tempvalxx = aux_arr;
  1701. tempvalxy=tempvalxx+SIMD_LEN;
  1702. tempvalxz=tempvalxy+SIMD_LEN;
  1703. tempvalyx=tempvalxz+SIMD_LEN;
  1704. tempvalyy=tempvalyx+SIMD_LEN;
  1705. tempvalyz=tempvalyy+SIMD_LEN;
  1706. tempvalzx=tempvalyz+SIMD_LEN;
  1707. tempvalzy=tempvalzx+SIMD_LEN;
  1708. tempvalzz=tempvalzy+SIMD_LEN;
  1709. /*! One over eight pi */
  1710. __m128d tofp = _mm_set1_pd (TOFP);
  1711. __m128d half = _mm_set1_pd (0.5);
  1712. __m128d opf = _mm_set1_pd (1.5);
  1713. __m128d zero = _mm_setzero_pd ();
  1714. // loop over sources
  1715. int i = 0;
  1716. for (; i < nt; i++) {
  1717. tempxx = _mm_setzero_pd(); tempxy = _mm_setzero_pd(); tempxz = _mm_setzero_pd();
  1718. tempyx = _mm_setzero_pd(); tempyy = _mm_setzero_pd(); tempyz = _mm_setzero_pd();
  1719. tempzx = _mm_setzero_pd(); tempzy = _mm_setzero_pd(); tempzz = _mm_setzero_pd();
  1720. __m128d txi = _mm_load1_pd (&tx[i]);
  1721. __m128d tyi = _mm_load1_pd (&ty[i]);
  1722. __m128d tzi = _mm_load1_pd (&tz[i]);
  1723. int j = 0;
  1724. // Load and calculate in groups of SIMD_LEN
  1725. for (; j + SIMD_LEN <= ns; j+=SIMD_LEN) {
  1726. __m128d sxj = _mm_load_pd (&sx[j]);
  1727. __m128d syj = _mm_load_pd (&sy[j]);
  1728. __m128d szj = _mm_load_pd (&sz[j]);
  1729. __m128d sdenx = _mm_set_pd (srcDen[(j+1)*3], srcDen[j*3]);
  1730. __m128d sdeny = _mm_set_pd (srcDen[(j+1)*3+1], srcDen[j*3+1]);
  1731. __m128d sdenz = _mm_set_pd (srcDen[(j+1)*3+2], srcDen[j*3+2]);
  1732. __m128d dX, dY, dZ;
  1733. __m128d dR2;
  1734. __m128d S;
  1735. __m128d S2;
  1736. dX = _mm_sub_pd(txi , sxj);
  1737. dY = _mm_sub_pd(tyi , syj);
  1738. dZ = _mm_sub_pd(tzi , szj);
  1739. sxj = _mm_mul_pd(dX, dX);
  1740. syj = _mm_mul_pd(dY, dY);
  1741. szj = _mm_mul_pd(dZ, dZ);
  1742. dR2 = _mm_add_pd(sxj, syj);
  1743. dR2 = _mm_add_pd(szj, dR2);
  1744. __m128d temp = _mm_cmpeq_pd (dR2, zero);
  1745. __m128d xhalf = _mm_mul_pd (half, dR2);
  1746. __m128 dR2_s = _mm_cvtpd_ps(dR2);
  1747. __m128 S_s = _mm_rsqrt_ps(dR2_s);
  1748. __m128d S_d = _mm_cvtps_pd(S_s);
  1749. // To handle the condition when src and trg coincide
  1750. S_d = _mm_andnot_pd (temp, S_d);
  1751. S = _mm_mul_pd (S_d, S_d);
  1752. S = _mm_mul_pd (S, xhalf);
  1753. S = _mm_sub_pd (opf, S);
  1754. S = _mm_mul_pd (S, S_d);
  1755. S2 = _mm_mul_pd (S, S);
  1756. __m128d dotx = _mm_mul_pd (dX, sdenx);
  1757. __m128d doty = _mm_mul_pd (dY, sdeny);
  1758. __m128d dotz = _mm_mul_pd (dZ, sdenz);
  1759. __m128d dot_sum = _mm_add_pd (dotx, doty);
  1760. dot_sum = _mm_add_pd (dot_sum, dotz);
  1761. dot_sum = _mm_mul_pd (dot_sum, S);
  1762. dot_sum = _mm_mul_pd (dot_sum, S2);
  1763. dot_sum = _mm_mul_pd (dot_sum, S2);
  1764. dotx = _mm_mul_pd (dot_sum, dX);
  1765. doty = _mm_mul_pd (dot_sum, dY);
  1766. dotz = _mm_mul_pd (dot_sum, dZ);
  1767. tempxx = _mm_add_pd (_mm_mul_pd(dotx,dX), tempxx);
  1768. tempxy = _mm_add_pd (_mm_mul_pd(dotx,dY), tempxy);
  1769. tempxz = _mm_add_pd (_mm_mul_pd(dotx,dZ), tempxz);
  1770. tempyx = _mm_add_pd (_mm_mul_pd(doty,dX), tempyx);
  1771. tempyy = _mm_add_pd (_mm_mul_pd(doty,dY), tempyy);
  1772. tempyz = _mm_add_pd (_mm_mul_pd(doty,dZ), tempyz);
  1773. tempzx = _mm_add_pd (_mm_mul_pd(dotz,dX), tempzx);
  1774. tempzy = _mm_add_pd (_mm_mul_pd(dotz,dY), tempzy);
  1775. tempzz = _mm_add_pd (_mm_mul_pd(dotz,dZ), tempzz);
  1776. }
  1777. tempxx = _mm_mul_pd (tempxx, tofp);
  1778. tempxy = _mm_mul_pd (tempxy, tofp);
  1779. tempxz = _mm_mul_pd (tempxz, tofp);
  1780. tempyx = _mm_mul_pd (tempyx, tofp);
  1781. tempyy = _mm_mul_pd (tempyy, tofp);
  1782. tempyz = _mm_mul_pd (tempyz, tofp);
  1783. tempzx = _mm_mul_pd (tempzx, tofp);
  1784. tempzy = _mm_mul_pd (tempzy, tofp);
  1785. tempzz = _mm_mul_pd (tempzz, tofp);
  1786. _mm_store_pd(tempvalxx, tempxx); _mm_store_pd(tempvalxy, tempxy); _mm_store_pd(tempvalxz, tempxz);
  1787. _mm_store_pd(tempvalyx, tempyx); _mm_store_pd(tempvalyy, tempyy); _mm_store_pd(tempvalyz, tempyz);
  1788. _mm_store_pd(tempvalzx, tempzx); _mm_store_pd(tempvalzy, tempzy); _mm_store_pd(tempvalzz, tempzz);
  1789. for (int k = 0; k < SIMD_LEN; k++) {
  1790. trgVal[i*9 ] += tempvalxx[k];
  1791. trgVal[i*9+1] += tempvalxy[k];
  1792. trgVal[i*9+2] += tempvalxz[k];
  1793. trgVal[i*9+3] += tempvalyx[k];
  1794. trgVal[i*9+4] += tempvalyy[k];
  1795. trgVal[i*9+5] += tempvalyz[k];
  1796. trgVal[i*9+6] += tempvalzx[k];
  1797. trgVal[i*9+7] += tempvalzy[k];
  1798. trgVal[i*9+8] += tempvalzz[k];
  1799. }
  1800. for (; j < ns; j++) {
  1801. double x = tx[i] - sx[j];
  1802. double y = ty[i] - sy[j];
  1803. double z = tz[i] - sz[j];
  1804. double r2 = x*x + y*y + z*z;
  1805. double r = pvfmm::sqrt<T>(r2);
  1806. double invdr;
  1807. if (r == 0)
  1808. invdr = 0;
  1809. else
  1810. invdr = 1/r;
  1811. double invdr2=invdr*invdr;
  1812. double dot = (x*srcDen[j*3] + y*srcDen[j*3+1] + z*srcDen[j*3+2]) * invdr2 * invdr2 * invdr;
  1813. double denx = dot*x;
  1814. double deny = dot*y;
  1815. double denz = dot*z;
  1816. trgVal[i*9 ] += denx*x*TOFP;
  1817. trgVal[i*9+1] += denx*y*TOFP;
  1818. trgVal[i*9+2] += denx*z*TOFP;
  1819. trgVal[i*9+3] += deny*x*TOFP;
  1820. trgVal[i*9+4] += deny*y*TOFP;
  1821. trgVal[i*9+5] += deny*z*TOFP;
  1822. trgVal[i*9+6] += denz*x*TOFP;
  1823. trgVal[i*9+7] += denz*y*TOFP;
  1824. trgVal[i*9+8] += denz*z*TOFP;
  1825. }
  1826. }
  1827. return;
  1828. }
  1829. void stokesGradSSE(
  1830. const int ns,
  1831. const int nt,
  1832. const double *sx,
  1833. const double *sy,
  1834. const double *sz,
  1835. const double *tx,
  1836. const double *ty,
  1837. const double *tz,
  1838. const double *srcDen,
  1839. double *trgVal,
  1840. const double cof )
  1841. {
  1842. if ( size_t(sx)%IDEAL_ALIGNMENT || size_t(sy)%IDEAL_ALIGNMENT || size_t(sz)%IDEAL_ALIGNMENT )
  1843. abort();
  1844. double mu = cof;
  1845. double OOEP = 1.0/(8.0*const_pi<double>());
  1846. __m128d tempxx; __m128d tempxy; __m128d tempxz;
  1847. __m128d tempyx; __m128d tempyy; __m128d tempyz;
  1848. __m128d tempzx; __m128d tempzy; __m128d tempzz;
  1849. double oomeu = 1/mu;
  1850. double aux_arr[9*SIMD_LEN+1];
  1851. double *tempvalxx, *tempvalxy, *tempvalxz;
  1852. double *tempvalyx, *tempvalyy, *tempvalyz;
  1853. double *tempvalzx, *tempvalzy, *tempvalzz;
  1854. if (size_t(aux_arr)%IDEAL_ALIGNMENT) // if aux_arr is misaligned
  1855. {
  1856. tempvalxx = aux_arr + 1;
  1857. if (size_t(tempvalxx)%IDEAL_ALIGNMENT)
  1858. abort();
  1859. }
  1860. else
  1861. tempvalxx = aux_arr;
  1862. tempvalxy=tempvalxx+SIMD_LEN;
  1863. tempvalxz=tempvalxy+SIMD_LEN;
  1864. tempvalyx=tempvalxz+SIMD_LEN;
  1865. tempvalyy=tempvalyx+SIMD_LEN;
  1866. tempvalyz=tempvalyy+SIMD_LEN;
  1867. tempvalzx=tempvalyz+SIMD_LEN;
  1868. tempvalzy=tempvalzx+SIMD_LEN;
  1869. tempvalzz=tempvalzy+SIMD_LEN;
  1870. /*! One over eight pi */
  1871. __m128d ooep = _mm_set1_pd (OOEP);
  1872. __m128d half = _mm_set1_pd (0.5);
  1873. __m128d opf = _mm_set1_pd (1.5);
  1874. __m128d three = _mm_set1_pd (3.0);
  1875. __m128d zero = _mm_setzero_pd ();
  1876. __m128d oomu = _mm_set1_pd (1/mu);
  1877. __m128d ooepmu = _mm_mul_pd(ooep,oomu);
  1878. // loop over sources
  1879. int i = 0;
  1880. for (; i < nt; i++) {
  1881. tempxx = _mm_setzero_pd(); tempxy = _mm_setzero_pd(); tempxz = _mm_setzero_pd();
  1882. tempyx = _mm_setzero_pd(); tempyy = _mm_setzero_pd(); tempyz = _mm_setzero_pd();
  1883. tempzx = _mm_setzero_pd(); tempzy = _mm_setzero_pd(); tempzz = _mm_setzero_pd();
  1884. __m128d txi = _mm_load1_pd (&tx[i]);
  1885. __m128d tyi = _mm_load1_pd (&ty[i]);
  1886. __m128d tzi = _mm_load1_pd (&tz[i]);
  1887. int j = 0;
  1888. // Load and calculate in groups of SIMD_LEN
  1889. for (; j + SIMD_LEN <= ns; j+=SIMD_LEN) {
  1890. __m128d sxj = _mm_load_pd (&sx[j]);
  1891. __m128d syj = _mm_load_pd (&sy[j]);
  1892. __m128d szj = _mm_load_pd (&sz[j]);
  1893. __m128d sdenx = _mm_set_pd (srcDen[(j+1)*3], srcDen[j*3]);
  1894. __m128d sdeny = _mm_set_pd (srcDen[(j+1)*3+1], srcDen[j*3+1]);
  1895. __m128d sdenz = _mm_set_pd (srcDen[(j+1)*3+2], srcDen[j*3+2]);
  1896. __m128d dX, dY, dZ;
  1897. __m128d dR2;
  1898. __m128d S;
  1899. __m128d S2;
  1900. __m128d S3;
  1901. dX = _mm_sub_pd(txi , sxj);
  1902. dY = _mm_sub_pd(tyi , syj);
  1903. dZ = _mm_sub_pd(tzi , szj);
  1904. sxj = _mm_mul_pd(dX, dX);
  1905. syj = _mm_mul_pd(dY, dY);
  1906. szj = _mm_mul_pd(dZ, dZ);
  1907. dR2 = _mm_add_pd(sxj, syj);
  1908. dR2 = _mm_add_pd(szj, dR2);
  1909. __m128d temp = _mm_cmpeq_pd (dR2, zero);
  1910. __m128d xhalf = _mm_mul_pd (half, dR2);
  1911. __m128 dR2_s = _mm_cvtpd_ps(dR2);
  1912. __m128 S_s = _mm_rsqrt_ps(dR2_s);
  1913. __m128d S_d = _mm_cvtps_pd(S_s);
  1914. // To handle the condition when src and trg coincide
  1915. S_d = _mm_andnot_pd (temp, S_d);
  1916. S = _mm_mul_pd (S_d, S_d);
  1917. S = _mm_mul_pd (S, xhalf);
  1918. S = _mm_sub_pd (opf, S);
  1919. S = _mm_mul_pd (S, S_d);
  1920. S2 = _mm_mul_pd (S, S);
  1921. S3 = _mm_mul_pd (S2, S);
  1922. __m128d dotx = _mm_mul_pd (dX, sdenx);
  1923. __m128d doty = _mm_mul_pd (dY, sdeny);
  1924. __m128d dotz = _mm_mul_pd (dZ, sdenz);
  1925. __m128d dot_sum = _mm_add_pd (dotx, doty);
  1926. dot_sum = _mm_add_pd (dot_sum, dotz);
  1927. dot_sum = _mm_mul_pd (dot_sum, S2);
  1928. tempxx = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dX, sdenx), _mm_mul_pd(sdenx, dX)), _mm_mul_pd(dot_sum, _mm_sub_pd(dR2 , _mm_mul_pd(three, _mm_mul_pd(dX, dX)))))),tempxx);
  1929. tempxy = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dY, sdenx), _mm_mul_pd(sdeny, dX)), _mm_mul_pd(dot_sum, _mm_sub_pd(zero, _mm_mul_pd(three, _mm_mul_pd(dY, dX)))))),tempxy);
  1930. tempxz = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dZ, sdenx), _mm_mul_pd(sdenz, dX)), _mm_mul_pd(dot_sum, _mm_sub_pd(zero, _mm_mul_pd(three, _mm_mul_pd(dZ, dX)))))),tempxz);
  1931. tempyx = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dX, sdeny), _mm_mul_pd(sdenx, dY)), _mm_mul_pd(dot_sum, _mm_sub_pd(zero, _mm_mul_pd(three, _mm_mul_pd(dX, dY)))))),tempyx);
  1932. tempyy = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dY, sdeny), _mm_mul_pd(sdeny, dY)), _mm_mul_pd(dot_sum, _mm_sub_pd(dR2 , _mm_mul_pd(three, _mm_mul_pd(dY, dY)))))),tempyy);
  1933. tempyz = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dZ, sdeny), _mm_mul_pd(sdenz, dY)), _mm_mul_pd(dot_sum, _mm_sub_pd(zero, _mm_mul_pd(three, _mm_mul_pd(dZ, dY)))))),tempyz);
  1934. tempzx = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dX, sdenz), _mm_mul_pd(sdenx, dZ)), _mm_mul_pd(dot_sum, _mm_sub_pd(zero, _mm_mul_pd(three, _mm_mul_pd(dX, dZ)))))),tempzx);
  1935. tempzy = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dY, sdenz), _mm_mul_pd(sdeny, dZ)), _mm_mul_pd(dot_sum, _mm_sub_pd(zero, _mm_mul_pd(three, _mm_mul_pd(dY, dZ)))))),tempzy);
  1936. tempzz = _mm_add_pd(_mm_mul_pd(S3,_mm_add_pd(_mm_sub_pd(_mm_mul_pd(dZ, sdenz), _mm_mul_pd(sdenz, dZ)), _mm_mul_pd(dot_sum, _mm_sub_pd(dR2 , _mm_mul_pd(three, _mm_mul_pd(dZ, dZ)))))),tempzz);
  1937. }
  1938. tempxx = _mm_mul_pd (tempxx, ooepmu);
  1939. tempxy = _mm_mul_pd (tempxy, ooepmu);
  1940. tempxz = _mm_mul_pd (tempxz, ooepmu);
  1941. tempyx = _mm_mul_pd (tempyx, ooepmu);
  1942. tempyy = _mm_mul_pd (tempyy, ooepmu);
  1943. tempyz = _mm_mul_pd (tempyz, ooepmu);
  1944. tempzx = _mm_mul_pd (tempzx, ooepmu);
  1945. tempzy = _mm_mul_pd (tempzy, ooepmu);
  1946. tempzz = _mm_mul_pd (tempzz, ooepmu);
  1947. _mm_store_pd(tempvalxx, tempxx); _mm_store_pd(tempvalxy, tempxy); _mm_store_pd(tempvalxz, tempxz);
  1948. _mm_store_pd(tempvalyx, tempyx); _mm_store_pd(tempvalyy, tempyy); _mm_store_pd(tempvalyz, tempyz);
  1949. _mm_store_pd(tempvalzx, tempzx); _mm_store_pd(tempvalzy, tempzy); _mm_store_pd(tempvalzz, tempzz);
  1950. for (int k = 0; k < SIMD_LEN; k++) {
  1951. trgVal[i*9 ] += tempvalxx[k];
  1952. trgVal[i*9+1] += tempvalxy[k];
  1953. trgVal[i*9+2] += tempvalxz[k];
  1954. trgVal[i*9+3] += tempvalyx[k];
  1955. trgVal[i*9+4] += tempvalyy[k];
  1956. trgVal[i*9+5] += tempvalyz[k];
  1957. trgVal[i*9+6] += tempvalzx[k];
  1958. trgVal[i*9+7] += tempvalzy[k];
  1959. trgVal[i*9+8] += tempvalzz[k];
  1960. }
  1961. for (; j < ns; j++) {
  1962. double x = tx[i] - sx[j];
  1963. double y = ty[i] - sy[j];
  1964. double z = tz[i] - sz[j];
  1965. double r2 = x*x + y*y + z*z;
  1966. double r = pvfmm::sqrt<T>(r2);
  1967. double invdr;
  1968. if (r == 0)
  1969. invdr = 0;
  1970. else
  1971. invdr = 1/r;
  1972. double invdr2=invdr*invdr;
  1973. double invdr3=invdr2*invdr;
  1974. double dot = (x*srcDen[j*3] + y*srcDen[j*3+1] + z*srcDen[j*3+2]);
  1975. trgVal[i*9 ] += OOEP*oomeu*invdr3*( x*srcDen[j*3 ] - srcDen[j*3 ]*x + dot*(1-3*x*x*invdr2) );
  1976. trgVal[i*9+1] += OOEP*oomeu*invdr3*( y*srcDen[j*3 ] - srcDen[j*3+1]*x + dot*(0-3*y*x*invdr2) );
  1977. trgVal[i*9+2] += OOEP*oomeu*invdr3*( z*srcDen[j*3 ] - srcDen[j*3+2]*x + dot*(0-3*z*x*invdr2) );
  1978. trgVal[i*9+3] += OOEP*oomeu*invdr3*( x*srcDen[j*3+1] - srcDen[j*3 ]*y + dot*(0-3*x*y*invdr2) );
  1979. trgVal[i*9+4] += OOEP*oomeu*invdr3*( y*srcDen[j*3+1] - srcDen[j*3+1]*y + dot*(1-3*y*y*invdr2) );
  1980. trgVal[i*9+5] += OOEP*oomeu*invdr3*( z*srcDen[j*3+1] - srcDen[j*3+2]*y + dot*(0-3*z*y*invdr2) );
  1981. trgVal[i*9+6] += OOEP*oomeu*invdr3*( x*srcDen[j*3+2] - srcDen[j*3 ]*z + dot*(0-3*x*z*invdr2) );
  1982. trgVal[i*9+7] += OOEP*oomeu*invdr3*( y*srcDen[j*3+2] - srcDen[j*3+1]*z + dot*(0-3*y*z*invdr2) );
  1983. trgVal[i*9+8] += OOEP*oomeu*invdr3*( z*srcDen[j*3+2] - srcDen[j*3+2]*z + dot*(1-3*z*z*invdr2) );
  1984. }
  1985. }
  1986. return;
  1987. }
  1988. #undef SIMD_LEN
  1989. #define X(s,k) (s)[(k)*COORD_DIM]
  1990. #define Y(s,k) (s)[(k)*COORD_DIM+1]
  1991. #define Z(s,k) (s)[(k)*COORD_DIM+2]
  1992. void stokesPressureSSEShuffle(const int ns, const int nt, double const src[], double const trg[], double const den[], double pot[], mem::MemoryManager* mem_mgr=NULL)
  1993. {
  1994. std::vector<double> xs(ns+1); std::vector<double> xt(nt);
  1995. std::vector<double> ys(ns+1); std::vector<double> yt(nt);
  1996. std::vector<double> zs(ns+1); std::vector<double> zt(nt);
  1997. int x_shift = size_t(&xs[0]) % IDEAL_ALIGNMENT ? 1:0;
  1998. int y_shift = size_t(&ys[0]) % IDEAL_ALIGNMENT ? 1:0;
  1999. int z_shift = size_t(&zs[0]) % IDEAL_ALIGNMENT ? 1:0;
  2000. //1. reshuffle memory
  2001. for (int k =0;k<ns;k++){
  2002. xs[k+x_shift]=X(src,k);
  2003. ys[k+y_shift]=Y(src,k);
  2004. zs[k+z_shift]=Z(src,k);
  2005. }
  2006. for (int k=0;k<nt;k++){
  2007. xt[k]=X(trg,k);
  2008. yt[k]=Y(trg,k);
  2009. zt[k]=Z(trg,k);
  2010. }
  2011. //2. perform caclulation
  2012. stokesPressureSSE(ns,nt,&xs[x_shift],&ys[y_shift],&zs[z_shift],&xt[0],&yt[0],&zt[0],den,pot);
  2013. return;
  2014. }
  2015. void stokesStressSSEShuffle(const int ns, const int nt, double const src[], double const trg[], double const den[], double pot[], mem::MemoryManager* mem_mgr=NULL)
  2016. {
  2017. std::vector<double> xs(ns+1); std::vector<double> xt(nt);
  2018. std::vector<double> ys(ns+1); std::vector<double> yt(nt);
  2019. std::vector<double> zs(ns+1); std::vector<double> zt(nt);
  2020. int x_shift = size_t(&xs[0]) % IDEAL_ALIGNMENT ? 1:0;
  2021. int y_shift = size_t(&ys[0]) % IDEAL_ALIGNMENT ? 1:0;
  2022. int z_shift = size_t(&zs[0]) % IDEAL_ALIGNMENT ? 1:0;
  2023. //1. reshuffle memory
  2024. for (int k =0;k<ns;k++){
  2025. xs[k+x_shift]=X(src,k);
  2026. ys[k+y_shift]=Y(src,k);
  2027. zs[k+z_shift]=Z(src,k);
  2028. }
  2029. for (int k=0;k<nt;k++){
  2030. xt[k]=X(trg,k);
  2031. yt[k]=Y(trg,k);
  2032. zt[k]=Z(trg,k);
  2033. }
  2034. //2. perform caclulation
  2035. stokesStressSSE(ns,nt,&xs[x_shift],&ys[y_shift],&zs[z_shift],&xt[0],&yt[0],&zt[0],den,pot);
  2036. return;
  2037. }
  2038. void stokesGradSSEShuffle(const int ns, const int nt, double const src[], double const trg[], double const den[], double pot[], const double kernel_coef, mem::MemoryManager* mem_mgr=NULL)
  2039. {
  2040. std::vector<double> xs(ns+1); std::vector<double> xt(nt);
  2041. std::vector<double> ys(ns+1); std::vector<double> yt(nt);
  2042. std::vector<double> zs(ns+1); std::vector<double> zt(nt);
  2043. int x_shift = size_t(&xs[0]) % IDEAL_ALIGNMENT ? 1:0;
  2044. int y_shift = size_t(&ys[0]) % IDEAL_ALIGNMENT ? 1:0;
  2045. int z_shift = size_t(&zs[0]) % IDEAL_ALIGNMENT ? 1:0;
  2046. //1. reshuffle memory
  2047. for (int k =0;k<ns;k++){
  2048. xs[k+x_shift]=X(src,k);
  2049. ys[k+y_shift]=Y(src,k);
  2050. zs[k+z_shift]=Z(src,k);
  2051. }
  2052. for (int k=0;k<nt;k++){
  2053. xt[k]=X(trg,k);
  2054. yt[k]=Y(trg,k);
  2055. zt[k]=Z(trg,k);
  2056. }
  2057. //2. perform caclulation
  2058. stokesGradSSE(ns,nt,&xs[x_shift],&ys[y_shift],&zs[z_shift],&xt[0],&yt[0],&zt[0],den,pot,kernel_coef);
  2059. return;
  2060. }
  2061. #undef X
  2062. #undef Y
  2063. #undef Z
  2064. #undef IDEAL_ALIGNMENT
  2065. #undef DECL_SIMD_ALIGNED
  2066. }
  2067. template <>
  2068. inline void stokes_press<double>(double* r_src, int src_cnt, double* v_src_, int dof, double* r_trg, int trg_cnt, double* k_out, mem::MemoryManager* mem_mgr){
  2069. Profile::Add_FLOP((long long)trg_cnt*(long long)src_cnt*(17*dof));
  2070. stokesPressureSSEShuffle(src_cnt, trg_cnt, r_src, r_trg, v_src_, k_out, mem_mgr);
  2071. return;
  2072. }
  2073. template <>
  2074. inline void stokes_stress<double>(double* r_src, int src_cnt, double* v_src_, int dof, double* r_trg, int trg_cnt, double* k_out, mem::MemoryManager* mem_mgr){
  2075. Profile::Add_FLOP((long long)trg_cnt*(long long)src_cnt*(45*dof));
  2076. stokesStressSSEShuffle(src_cnt, trg_cnt, r_src, r_trg, v_src_, k_out, mem_mgr);
  2077. }
  2078. template <>
  2079. inline void stokes_grad<double>(double* r_src, int src_cnt, double* v_src_, int dof, double* r_trg, int trg_cnt, double* k_out, mem::MemoryManager* mem_mgr){
  2080. Profile::Add_FLOP((long long)trg_cnt*(long long)src_cnt*(89*dof));
  2081. const double mu=1.0;
  2082. stokesGradSSEShuffle(src_cnt, trg_cnt, r_src, r_trg, v_src_, k_out, mu, mem_mgr);
  2083. }
  2084. #endif
  2085. #endif
  2086. template<class T> const Kernel<T>& StokesKernel<T>::velocity(){
  2087. static Kernel<T> ker=BuildKernel<T, stokes_vel<T,1>, stokes_sym_dip>("stokes_vel" , 3, std::pair<int,int>(3,3),
  2088. NULL,NULL,NULL, NULL,NULL,NULL, NULL,NULL, &stokes_vol_poten<T>);
  2089. return ker;
  2090. }
  2091. template<class T> const Kernel<T>& StokesKernel<T>::pressure(){
  2092. static Kernel<T> ker=BuildKernel<T, stokes_press >("stokes_press" , 3, std::pair<int,int>(3,1));
  2093. return ker;
  2094. }
  2095. template<class T> const Kernel<T>& StokesKernel<T>::stress(){
  2096. static Kernel<T> ker=BuildKernel<T, stokes_stress >("stokes_stress", 3, std::pair<int,int>(3,9));
  2097. return ker;
  2098. }
  2099. template<class T> const Kernel<T>& StokesKernel<T>::vel_grad(){
  2100. static Kernel<T> ker=BuildKernel<T, stokes_grad >("stokes_grad" , 3, std::pair<int,int>(3,9));
  2101. return ker;
  2102. }
  2103. template<> inline const Kernel<double>& StokesKernel<double>::velocity(){
  2104. typedef double T;
  2105. static Kernel<T> ker=BuildKernel<T, stokes_vel<T,2>, stokes_sym_dip>("stokes_vel" , 3, std::pair<int,int>(3,3),
  2106. NULL,NULL,NULL, NULL,NULL,NULL, NULL,NULL, &stokes_vol_poten<double>);
  2107. return ker;
  2108. }
  2109. ////////////////////////////////////////////////////////////////////////////////
  2110. //////// BIOT-SAVART KERNEL ////////
  2111. ////////////////////////////////////////////////////////////////////////////////
  2112. template <class Real_t, class Vec_t=Real_t, Vec_t (*RSQRT_INTRIN)(Vec_t)=rsqrt_intrin0<Vec_t> >
  2113. void biot_savart_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value, Matrix<Real_t>& trg_coord, Matrix<Real_t>& trg_value){
  2114. #define SRC_BLK 500
  2115. size_t VecLen=sizeof(Vec_t)/sizeof(Real_t);
  2116. //// Number of newton iterations
  2117. size_t NWTN_ITER=0;
  2118. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin0<Vec_t,Real_t>) NWTN_ITER=0;
  2119. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin1<Vec_t,Real_t>) NWTN_ITER=1;
  2120. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin2<Vec_t,Real_t>) NWTN_ITER=2;
  2121. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin3<Vec_t,Real_t>) NWTN_ITER=3;
  2122. Real_t nwtn_scal=1; // scaling factor for newton iterations
  2123. for(int i=0;i<NWTN_ITER;i++){
  2124. nwtn_scal=2*nwtn_scal*nwtn_scal*nwtn_scal;
  2125. }
  2126. const Real_t OOFP = 1.0/(4*nwtn_scal*nwtn_scal*nwtn_scal*const_pi<Real_t>());
  2127. size_t src_cnt_=src_coord.Dim(1);
  2128. size_t trg_cnt_=trg_coord.Dim(1);
  2129. for(size_t sblk=0;sblk<src_cnt_;sblk+=SRC_BLK){
  2130. size_t src_cnt=src_cnt_-sblk;
  2131. if(src_cnt>SRC_BLK) src_cnt=SRC_BLK;
  2132. for(size_t t=0;t<trg_cnt_;t+=VecLen){
  2133. Vec_t tx=load_intrin<Vec_t>(&trg_coord[0][t]);
  2134. Vec_t ty=load_intrin<Vec_t>(&trg_coord[1][t]);
  2135. Vec_t tz=load_intrin<Vec_t>(&trg_coord[2][t]);
  2136. Vec_t tvx=zero_intrin<Vec_t>();
  2137. Vec_t tvy=zero_intrin<Vec_t>();
  2138. Vec_t tvz=zero_intrin<Vec_t>();
  2139. for(size_t s=sblk;s<sblk+src_cnt;s++){
  2140. Vec_t dx=sub_intrin(tx,bcast_intrin<Vec_t>(&src_coord[0][s]));
  2141. Vec_t dy=sub_intrin(ty,bcast_intrin<Vec_t>(&src_coord[1][s]));
  2142. Vec_t dz=sub_intrin(tz,bcast_intrin<Vec_t>(&src_coord[2][s]));
  2143. Vec_t svx= bcast_intrin<Vec_t>(&src_value[0][s]) ;
  2144. Vec_t svy= bcast_intrin<Vec_t>(&src_value[1][s]) ;
  2145. Vec_t svz= bcast_intrin<Vec_t>(&src_value[2][s]) ;
  2146. Vec_t r2= mul_intrin(dx,dx) ;
  2147. r2=add_intrin(r2,mul_intrin(dy,dy));
  2148. r2=add_intrin(r2,mul_intrin(dz,dz));
  2149. Vec_t rinv=RSQRT_INTRIN(r2);
  2150. Vec_t rinv3=mul_intrin(mul_intrin(rinv,rinv),rinv);
  2151. tvx=sub_intrin(tvx,mul_intrin(rinv3,sub_intrin(mul_intrin(svy,dz),mul_intrin(svz,dy))));
  2152. tvy=sub_intrin(tvy,mul_intrin(rinv3,sub_intrin(mul_intrin(svz,dx),mul_intrin(svx,dz))));
  2153. tvz=sub_intrin(tvz,mul_intrin(rinv3,sub_intrin(mul_intrin(svx,dy),mul_intrin(svy,dx))));
  2154. }
  2155. Vec_t oofp=set_intrin<Vec_t,Real_t>(OOFP);
  2156. tvx=add_intrin(mul_intrin(tvx,oofp),load_intrin<Vec_t>(&trg_value[0][t]));
  2157. tvy=add_intrin(mul_intrin(tvy,oofp),load_intrin<Vec_t>(&trg_value[1][t]));
  2158. tvz=add_intrin(mul_intrin(tvz,oofp),load_intrin<Vec_t>(&trg_value[2][t]));
  2159. store_intrin(&trg_value[0][t],tvx);
  2160. store_intrin(&trg_value[1][t],tvy);
  2161. store_intrin(&trg_value[2][t],tvz);
  2162. }
  2163. }
  2164. { // Add FLOPS
  2165. #ifndef __MIC__
  2166. Profile::Add_FLOP((long long)trg_cnt_*(long long)src_cnt_*(29+4*(NWTN_ITER)));
  2167. #endif
  2168. }
  2169. #undef SRC_BLK
  2170. }
  2171. template <class T, int newton_iter=0>
  2172. void biot_savart(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
  2173. #define BS_KER_NWTN(nwtn) if(newton_iter==nwtn) \
  2174. generic_kernel<Real_t, 3, 3, biot_savart_uKernel<Real_t,Vec_t, rsqrt_intrin##nwtn<Vec_t,Real_t> > > \
  2175. ((Real_t*)r_src, src_cnt, (Real_t*)v_src, dof, (Real_t*)r_trg, trg_cnt, (Real_t*)v_trg, mem_mgr)
  2176. #define BIOTSAVART_KERNEL BS_KER_NWTN(0); BS_KER_NWTN(1); BS_KER_NWTN(2); BS_KER_NWTN(3);
  2177. if(mem::TypeTraits<T>::ID()==mem::TypeTraits<float>::ID()){
  2178. typedef float Real_t;
  2179. #if defined __MIC__
  2180. #define Vec_t Real_t
  2181. #elif defined __AVX__
  2182. #define Vec_t __m256
  2183. #elif defined __SSE3__
  2184. #define Vec_t __m128
  2185. #else
  2186. #define Vec_t Real_t
  2187. #endif
  2188. BIOTSAVART_KERNEL;
  2189. #undef Vec_t
  2190. }else if(mem::TypeTraits<T>::ID()==mem::TypeTraits<double>::ID()){
  2191. typedef double Real_t;
  2192. #if defined __MIC__
  2193. #define Vec_t Real_t
  2194. #elif defined __AVX__
  2195. #define Vec_t __m256d
  2196. #elif defined __SSE3__
  2197. #define Vec_t __m128d
  2198. #else
  2199. #define Vec_t Real_t
  2200. #endif
  2201. BIOTSAVART_KERNEL;
  2202. #undef Vec_t
  2203. }else{
  2204. typedef T Real_t;
  2205. #define Vec_t Real_t
  2206. BIOTSAVART_KERNEL;
  2207. #undef Vec_t
  2208. }
  2209. #undef BS_KER_NWTN
  2210. #undef BIOTSAVART_KERNEL
  2211. }
  2212. template<class T> const Kernel<T>& BiotSavartKernel<T>::potential(){
  2213. static Kernel<T> ker=BuildKernel<T, biot_savart<T,1> >("biot_savart", 3, std::pair<int,int>(3,3));
  2214. return ker;
  2215. }
  2216. template<> inline const Kernel<double>& BiotSavartKernel<double>::potential(){
  2217. typedef double T;
  2218. static Kernel<T> ker=BuildKernel<T, biot_savart<T,2> >("biot_savart", 3, std::pair<int,int>(3,3));
  2219. return ker;
  2220. }
  2221. ////////////////////////////////////////////////////////////////////////////////
  2222. //////// HELMHOLTZ KERNEL ////////
  2223. ////////////////////////////////////////////////////////////////////////////////
  2224. /**
  2225. * \brief Green's function for the Helmholtz's equation. Kernel tensor
  2226. * dimension = 2x2.
  2227. */
  2228. template <class Real_t, class Vec_t=Real_t, Vec_t (*RSQRT_INTRIN)(Vec_t)=rsqrt_intrin0<Vec_t> >
  2229. void helmholtz_poten_uKernel(Matrix<Real_t>& src_coord, Matrix<Real_t>& src_value, Matrix<Real_t>& trg_coord, Matrix<Real_t>& trg_value){
  2230. #define SRC_BLK 500
  2231. size_t VecLen=sizeof(Vec_t)/sizeof(Real_t);
  2232. //// Number of newton iterations
  2233. size_t NWTN_ITER=0;
  2234. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin0<Vec_t,Real_t>) NWTN_ITER=0;
  2235. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin1<Vec_t,Real_t>) NWTN_ITER=1;
  2236. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin2<Vec_t,Real_t>) NWTN_ITER=2;
  2237. if(RSQRT_INTRIN==(Vec_t (*)(Vec_t))rsqrt_intrin3<Vec_t,Real_t>) NWTN_ITER=3;
  2238. Real_t nwtn_scal=1; // scaling factor for newton iterations
  2239. for(int i=0;i<NWTN_ITER;i++){
  2240. nwtn_scal=2*nwtn_scal*nwtn_scal*nwtn_scal;
  2241. }
  2242. const Real_t OOFP = 1.0/(4*nwtn_scal*const_pi<Real_t>());
  2243. const Vec_t mu = set_intrin<Vec_t,Real_t>(20.0*const_pi<Real_t>()/nwtn_scal);
  2244. size_t src_cnt_=src_coord.Dim(1);
  2245. size_t trg_cnt_=trg_coord.Dim(1);
  2246. for(size_t sblk=0;sblk<src_cnt_;sblk+=SRC_BLK){
  2247. size_t src_cnt=src_cnt_-sblk;
  2248. if(src_cnt>SRC_BLK) src_cnt=SRC_BLK;
  2249. for(size_t t=0;t<trg_cnt_;t+=VecLen){
  2250. Vec_t tx=load_intrin<Vec_t>(&trg_coord[0][t]);
  2251. Vec_t ty=load_intrin<Vec_t>(&trg_coord[1][t]);
  2252. Vec_t tz=load_intrin<Vec_t>(&trg_coord[2][t]);
  2253. Vec_t tvx=zero_intrin<Vec_t>();
  2254. Vec_t tvy=zero_intrin<Vec_t>();
  2255. for(size_t s=sblk;s<sblk+src_cnt;s++){
  2256. Vec_t dx=sub_intrin(tx,bcast_intrin<Vec_t>(&src_coord[0][s]));
  2257. Vec_t dy=sub_intrin(ty,bcast_intrin<Vec_t>(&src_coord[1][s]));
  2258. Vec_t dz=sub_intrin(tz,bcast_intrin<Vec_t>(&src_coord[2][s]));
  2259. Vec_t svx= bcast_intrin<Vec_t>(&src_value[0][s]) ;
  2260. Vec_t svy= bcast_intrin<Vec_t>(&src_value[1][s]) ;
  2261. Vec_t r2= mul_intrin(dx,dx) ;
  2262. r2=add_intrin(r2,mul_intrin(dy,dy));
  2263. r2=add_intrin(r2,mul_intrin(dz,dz));
  2264. Vec_t rinv=RSQRT_INTRIN(r2);
  2265. Vec_t mu_r=mul_intrin(mu,mul_intrin(r2,rinv));
  2266. Vec_t G0=mul_intrin(cos_intrin(mu_r),rinv);
  2267. Vec_t G1=mul_intrin(sin_intrin(mu_r),rinv);
  2268. tvx=add_intrin(tvx,sub_intrin(mul_intrin(svx,G0),mul_intrin(svy,G1)));
  2269. tvy=add_intrin(tvy,add_intrin(mul_intrin(svx,G1),mul_intrin(svy,G0)));
  2270. }
  2271. Vec_t oofp=set_intrin<Vec_t,Real_t>(OOFP);
  2272. tvx=add_intrin(mul_intrin(tvx,oofp),load_intrin<Vec_t>(&trg_value[0][t]));
  2273. tvy=add_intrin(mul_intrin(tvy,oofp),load_intrin<Vec_t>(&trg_value[1][t]));
  2274. store_intrin(&trg_value[0][t],tvx);
  2275. store_intrin(&trg_value[1][t],tvy);
  2276. }
  2277. }
  2278. { // Add FLOPS
  2279. #ifndef __MIC__
  2280. Profile::Add_FLOP((long long)trg_cnt_*(long long)src_cnt_*(24+4*(NWTN_ITER)));
  2281. #endif
  2282. }
  2283. #undef SRC_BLK
  2284. }
  2285. template <class T, int newton_iter=0>
  2286. void helmholtz_poten(T* r_src, int src_cnt, T* v_src, int dof, T* r_trg, int trg_cnt, T* v_trg, mem::MemoryManager* mem_mgr){
  2287. #define HELM_KER_NWTN(nwtn) if(newton_iter==nwtn) \
  2288. generic_kernel<Real_t, 2, 2, helmholtz_poten_uKernel<Real_t,Vec_t, rsqrt_intrin##nwtn<Vec_t,Real_t> > > \
  2289. ((Real_t*)r_src, src_cnt, (Real_t*)v_src, dof, (Real_t*)r_trg, trg_cnt, (Real_t*)v_trg, mem_mgr)
  2290. #define HELMHOLTZ_KERNEL HELM_KER_NWTN(0); HELM_KER_NWTN(1); HELM_KER_NWTN(2); HELM_KER_NWTN(3);
  2291. if(mem::TypeTraits<T>::ID()==mem::TypeTraits<float>::ID()){
  2292. typedef float Real_t;
  2293. #if defined __MIC__
  2294. #define Vec_t Real_t
  2295. #elif defined __AVX__
  2296. #define Vec_t __m256
  2297. #elif defined __SSE3__
  2298. #define Vec_t __m128
  2299. #else
  2300. #define Vec_t Real_t
  2301. #endif
  2302. HELMHOLTZ_KERNEL;
  2303. #undef Vec_t
  2304. }else if(mem::TypeTraits<T>::ID()==mem::TypeTraits<double>::ID()){
  2305. typedef double Real_t;
  2306. #if defined __MIC__
  2307. #define Vec_t Real_t
  2308. #elif defined __AVX__
  2309. #define Vec_t __m256d
  2310. #elif defined __SSE3__
  2311. #define Vec_t __m128d
  2312. #else
  2313. #define Vec_t Real_t
  2314. #endif
  2315. HELMHOLTZ_KERNEL;
  2316. #undef Vec_t
  2317. }else{
  2318. typedef T Real_t;
  2319. #define Vec_t Real_t
  2320. HELMHOLTZ_KERNEL;
  2321. #undef Vec_t
  2322. }
  2323. #undef HELM_KER_NWTN
  2324. #undef HELMHOLTZ_KERNEL
  2325. }
  2326. template<class T> const Kernel<T>& HelmholtzKernel<T>::potential(){
  2327. static Kernel<T> ker=BuildKernel<T, helmholtz_poten<T,1> >("helmholtz" , 3, std::pair<int,int>(2,2));
  2328. return ker;
  2329. }
  2330. template<> inline const Kernel<double>& HelmholtzKernel<double>::potential(){
  2331. typedef double T;
  2332. static Kernel<T> ker=BuildKernel<T, helmholtz_poten<T,3> >("helmholtz" , 3, std::pair<int,int>(2,2));
  2333. return ker;
  2334. }
  2335. }//end namespace