OpenAtom  Version1.5a
CP_State_RealParticlePlane.C
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 ///////////////////////////////////////////////////////////////////////////
3 ///////////////////////////////////////////////////////////////////////////
4 /** \file CP_State_RealParticlePlane.C
5  * Life-cycle of a CP_State_RealParticlePlane:
6  *
7  * Driven by the CP_State_GSpaceParticlePlane, this implements the
8  * Real space side of the fft based Euler Exponential Spline method
9  * for computing the non-local quantities.
10  */
11 ///////////////////////////////////////////////////////////////////////////
12 
13 #include "utility/util.h" ///@note: Putting this declaration further down seems to screw things. Whats the issue with fft malloc & free declarations?
14 
15 #include "CP_State_ParticlePlane.h"
16 #include "CP_State_Plane.h"
17 #include "structure_factor/StructFactorCache.h"
18 #include "fft_slab_ctrl/fftCacheSlab.h"
19 #include "main/AtomsCache.h"
20 #include "main/eesCache.h"
21 #include "main/cpaimd.h"
22 #include "utility/matrix2file.h"
23 
24 #include "ckmulticast.h"
25 #include "charm++.h"
26 #include "../../src_piny_physics_v1.0/include/class_defs/CP_OPERATIONS/class_cpnonlocal.h"
27 #include "../../src_piny_physics_v1.0/include/class_defs/CP_OPERATIONS/class_cplocal.h"
28 
29 ///////////////////////////////////////////////////////////////////////////
30 extern CProxy_InstanceController instControllerProxy;
31 extern CProxy_main mainProxy;
32 extern CkVec <CProxy_CP_State_GSpacePlane> UgSpacePlaneProxy;
33 extern CkVec <CProxy_AtomsCache> UatomsCacheProxy;
34 extern CkVec <CProxy_CP_State_ParticlePlane> UparticlePlaneProxy;
35 extern CkVec <CProxy_StructFactCache> UsfCacheProxy;
36 extern CkVec <CProxy_eesCache> UeesCacheProxy;
37 extern CkVec <CProxy_FFTcache> UfftCacheProxy;
38 extern CkVec <CProxy_CP_State_RealParticlePlane> UrealParticlePlaneProxy;
39 extern CkVec <MapType2> RPPImaptable;
40 extern CkGroupID mCastGrpId;
41 extern ComlibInstanceHandle mssPInstance;
42 extern CkReduction::reducerType sumFastDoubleType;
43 extern int nstates;
44 extern int nchareG;
45 extern Config config;
46 
47 //#define _CP_DEBUG_STATE_RPP_VERBOSE_
48 
49 ///////////////////////////////////////////////////////////////////////////
50 
51 /** @addtogroup Particle
52  @{
53 */
54 
55 //////////////////////////////////////////////////////////////////////////////
56 /// Energy reduction client for ees method!
57 //////////////////////////////////////////////////////////////////////////////
58 ///////////////////////////////////////////////////////////////////////////cc
59 //////////////////////////////////////////////////////////////////////////////
60 void CP_State_RealParticlePlane::printEnlR(CkReductionMsg *m){
61 
62  //unpack
63  double d = ((double *)m->getData())[0];
64  itimeRed = m->getUserFlag();
65  delete m;
66  //output and save the data
67  FILE *temperScreenFile = UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch()->temperScreenFile;
68  int iteration= UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch()->iteration;
69  fprintf(temperScreenFile,"Iter [%d] ENL(EES) = %5.8lf\n", iteration,d);
70  UgSpacePlaneProxy[thisInstance.proxyOffset](0,0).computeEnergies(ENERGY_ENL, d);
71 }
72 //////////////////////////////////////////////////////////////////////////////
73 
74 
75 //////////////////////////////////////////////////////////////////////////////
76 /// Energy reduction client for ees method!
77 //////////////////////////////////////////////////////////////////////////////
78 ///////////////////////////////////////////////////////////////////////////cc
79 //////////////////////////////////////////////////////////////////////////////
80 void CP_State_RealParticlePlane::printEnlRSimp(double cp_enl_loc,int index,int itime_in){
81 
82  // upack
83  countEnl++;
84  if(countEnl==1){itimeRed=itime_in;}
85  cp_enlTot += cp_enl_loc;
86 
87  if(itimeRed!=itime_in){
88  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
89  CkPrintf("Yo dawg, ees nonlocal is racin' into the enl reduction\n");
90  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
91  CkExit();
92  }//endif
93 
94 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
95  CkPrintf("Enl contrib arrived from %d : %g %d %d %d\n",index,cp_enl_loc,countEnl,
96  itime_in,itimeRed);
97 #endif
98 
99  // output and save the data
100  if(countEnl==nstates){
101  CkPrintf("{%d} ENL(EES) = %5.8lf\n", thisInstance.proxyOffset, cp_enlTot);
102  UgSpacePlaneProxy[thisInstance.proxyOffset](0,0).computeEnergies(ENERGY_ENL,cp_enlTot);
103  countEnl = 0;
104  cp_enlTot = 0.0;
105  }//endif
106 
107 }
108 //////////////////////////////////////////////////////////////////////////////
109 
110 
111 
112 //////////////////////////////////////////////////////////////////////////////
113 /// The constructor : Called only once
114 //////////////////////////////////////////////////////////////////////////////
115 //////////////////////////////////////////////////////////////////////////////
116 //////////////////////////////////////////////////////////////////////////////
117 CP_State_RealParticlePlane::CP_State_RealParticlePlane(
118  int ngridA_in, int ngridB_in, int ngridC_in,
119  int numIterNL_in, int zmatSizeMax_in,
120  int Rstates_per_pe_in,int nChareG_in,
121  int ees_nonlocal_in, UberCollection _instance):
122  thisInstance(_instance)
123 //////////////////////////////////////////////////////////////////////////////
124  {//begin routine
125 //////////////////////////////////////////////////////////////////////////////
126 /// The non-local guts variables
127 
128  ngridA = ngridA_in; // grid size along a
129  ngridB = ngridB_in; // grid size along b
130  ngridC = ngridC_in; // grid size along c
131 
132  nChareR = ngridC; // Real Space chares=# C-planes
133  nChareG = nChareG_in; // G Space chares
134  Rstates_per_pe = Rstates_per_pe_in; // Real Space topomap
135  myPlane = thisIndex.y; // Real space plane number
136  ibead_ind = thisInstance.idxU.x;
137  kpoint_ind = thisInstance.idxU.y;
138  itemper_ind = thisInstance.idxU.z;
139 
140  if(config.doublePack){
141  planeSize = (ngridA+2)*ngridB; // expanded plane size for FFTing
142  planeSizeT = ngridA*ngridB; // true plane size
143  csize = (ngridA+2)*ngridB/2; // complex variable size
144  }else{
145  planeSize = 2*ngridA*ngridB; // expanded plane size for FFTing
146  planeSizeT = 2*ngridA*ngridB; // true plane size
147  csize = ngridA*ngridB; // complex variable size
148  }//endif
149  CkAssert(csize*2==planeSize);
150  numIterNl = numIterNL_in; // # of non-local iterations per time step
151  zmatSizeMax = zmatSizeMax_in; // zmatrix size
152  ees_nonlocal = ees_nonlocal_in;
153 
154  rhoRTime = 0;
155  cp_enl = 0.0; // non-local energy
156  cp_enlTot = 0.0;
157  count = 0; // fft communication counter
158  countEnl = 0; // Energy counter
159  iterNL = 0; // Nl iteration counter
160  itimeRed = 0; // time of reduction
161  itime = 0; // time step;
162  registrationFlag = 0;
163  recvBlock = 0;
164  fftDataDone = false;
165  launchFFT = false;
166  countZ = 0; // Zmat communication counter
167 
168 //////////////////////////////////////////////////////////////////////////////
169 /// No migration: No atSync load-balancing act
170 
171  setMigratable(false);
172  usesAtSync = false;
173 #ifdef _CP_GS_DEBUG_COMPARE_VKS_
174  savedprojpsiC=NULL;
175  savedProjpsiCScr=NULL;
176  savedProjpsiRScr=NULL;
177  savedzmat=NULL;
178  savedmn=NULL;
179  saveddmn_x=NULL;
180  saveddmn_y=NULL;
181  saveddmn_z=NULL;
182  savedigrid=NULL;
183 #endif
184 //----------------------------------------------------------------------------
185  }//end routine
186 //////////////////////////////////////////////////////////////////////////////
187 
188 //////////////////////////////////////////////////////////////////////////////
189 /// Post construction initialization.
190 /// you can't make sections until the multicast manager is done
191 //////////////////////////////////////////////////////////////////////////////
192 //////////////////////////////////////////////////////////////////////////////
193 //////////////////////////////////////////////////////////////////////////////
194 void CP_State_RealParticlePlane::init(){
195 
196 //////////////////////////////////////////////////////////////////////////////
197 /// Malloc the projector memory, non-local matrix and register with your cache
198 
199  if(ees_nonlocal==1){
200  //--------
201  // malloc
202  projPsiC = (complex*) fftw_malloc(csize*sizeof(complex));
203  projPsiR = reinterpret_cast<double*> (projPsiC);
204  if(config.doublePack){
205  zmat = (double*) fftw_malloc(zmatSizeMax*sizeof(double));
206  zmatScr = (double*) fftw_malloc(zmatSizeMax*sizeof(double));
207  }else{
208  zmatC = (complex*) fftw_malloc(zmatSizeMax*sizeof(complex));
209  zmatScrC = (complex*) fftw_malloc(zmatSizeMax*sizeof(complex));
210  zmat = reinterpret_cast<double*> (zmatC);
211  zmatScr = reinterpret_cast<double*> (zmatScrC);
212  }//endif
213  //--------
214  // Register
215  eesCache *eesData = UeesCacheProxy[thisInstance.proxyOffset].ckLocalBranch ();
216  eesData->registerCacheRPP(thisIndex.y);
217  }//endif
218 
219 
220 //////////////////////////////////////////////////////////////////////////////
221 /// Choose reduction plane reasonably intelligently
222 
223  int *red_pl = new int[nstates];
224  //foreach state, try to find a new proc for each state's reduction plane
225  int *usedProc= new int[CkNumPes()];
226  memset(usedProc,0,sizeof(int)*CkNumPes());
227  int charperpe=nstates/CkNumPes();
228  if(nstates%CkNumPes()!=0) charperpe++;
229  if(charperpe<1) charperpe=1;
230  for(int state=0; state<nstates;state++){
231  int plane=0;
232  while(plane<nChareR){
233  bool used=false;
234  int thisstateplaneproc=RPPImaptable[thisInstance.proxyOffset].get(state,plane)%CkNumPes();
235  if(usedProc[thisstateplaneproc]>=charperpe);
236  {
237  used=true;
238  }
239  if(!used || (plane+1==nChareR)){
240  usedProc[thisstateplaneproc]++;
241  red_pl[state]=plane;
242  plane=nChareR;
243  }//endif
244  plane++;
245  }//end while : plane < nchareR
246  }//end for : i
247  reductionPlaneNum = red_pl[thisIndex.x];
248 
249  delete [] usedProc;
250 
251 
252 //////////////////////////////////////////////////////////////////////////////
253 /// Build section reductions
254 
255  //-----------------------------------------------------------
256  // The plane section reduction to the reduction plane for zmat
257  // Only the root of the reduction does setup dance.
258  //-----------------------------------------------------------
259  enlSectionComplete=false;
260  planeRedSectionComplete=false;
261  if(thisIndex.y==reductionPlaneNum){
262  rPlaneSectProxy =
263  CProxySection_CP_State_RealParticlePlane::ckNew(thisProxy.ckGetArrayID(),
264  thisIndex.x,thisIndex.x,1,
265  0,nChareR-1,1);
266  CkMulticastMgr *mcastGrp = CProxy_CkMulticastMgr(mCastGrpId).ckLocalBranch();
267  rPlaneSectProxy.ckDelegate(mcastGrp);
268  mcastGrp->setSection(rPlaneSectProxy);
269  EnlCookieMsg *emsg= new EnlCookieMsg;
270  rPlaneSectProxy.setPlaneRedCookie(emsg);
271  }//endif
272  else
273  { // if you aren't in the reduction plane you are setup
274  planeRedSectionComplete=true;
275  }
276 
277  //-----------------------------------------------------------
278  // The reduction over states involving the reduction planes for cp_enl
279  // Only the root of the reduction does setup dance, state=0
280  //-----------------------------------------------------------
281 
282  if(thisIndex.y==reductionPlaneNum && thisIndex.x==0){
283  CkArrayIndexMax *elems = new CkArrayIndexMax[nstates];
284  CkArrayIndex2D idx(0, reductionPlaneNum); // cheesy constructor
285  for (int j = 0; j < nstates; j++) {
286  idx.index[0] = j; idx.index[1] = red_pl[j];
287  elems[j] = idx;
288  }//endfor
289  rPlaneENLProxy =
290  CProxySection_CP_State_RealParticlePlane::ckNew(thisProxy.ckGetArrayID(),
291  elems, nstates);
292  CkMulticastMgr *mcastGrp = CProxy_CkMulticastMgr(mCastGrpId).ckLocalBranch();
293  rPlaneENLProxy.ckDelegate(mcastGrp);
294  mcastGrp->setSection(rPlaneENLProxy);
295  EnlCookieMsg *emsg= new EnlCookieMsg;
296  rPlaneENLProxy.setEnlCookie(emsg);
297  delete [] elems;
298  }//endif
299  else
300  { // if not part of enl, setup complete
301  bool inEnl=false;
302  for (int j = 0; j < nstates; j++) {
303  if(thisIndex.x ==j && thisIndex.y== red_pl[j])
304  inEnl=true;
305  }//endfor
306  if(!inEnl)
307  enlSectionComplete=true;
308  }
309  initComplete();
310  delete [] red_pl;
311 
312 //////////////////////////////////////////////////////////////////////////////
313 /// Setup the comlib to talk to GPP
314 
315  gPP_proxy = UparticlePlaneProxy[thisInstance.proxyOffset];
316 #ifdef USE_COMLIB
317  if (config.useMssInsGPP){
318  ComlibAssociateProxy(mssPInstance,gPP_proxy);
319  }//endif
320 #endif
321 
322 }
323 
324 //////////////////////////////////////////////////////////////////////////////
325 /// All cookies initialized for enl section, only reduction root receives this
326 //////////////////////////////////////////////////////////////////////////////
327 //////////////////////////////////////////////////////////////////////////////
328 //////////////////////////////////////////////////////////////////////////////
329 void CP_State_RealParticlePlane::enlSectDone(CkReductionMsg *m)
330 {
331  // CkPrintf("RPP[%d %d] gets enlSectDone\n",thisIndex.x, thisIndex.y);
332  delete m;
333  enlSectionComplete=true;
334  initComplete();
335 }
336 
337 //////////////////////////////////////////////////////////////////////////////
338 /// All cookies initialized for zmat plane reduction
339 //////////////////////////////////////////////////////////////////////////////
340 //////////////////////////////////////////////////////////////////////////////
341 //////////////////////////////////////////////////////////////////////////////
342 void CP_State_RealParticlePlane::planeRedSectDone(CkReductionMsg *m)
343 {
344  // CkPrintf("RPP[%d %d] gets rPlaneRedSectDone\n",thisIndex.x, thisIndex.y);
345  delete m;
346  planeRedSectionComplete=true;
347  initComplete();
348 }
349 
350 //////////////////////////////////////////////////////////////////////////////
351 /// Initialization and registration done for this chare
352 //////////////////////////////////////////////////////////////////////////////
353 //////////////////////////////////////////////////////////////////////////////
354 //////////////////////////////////////////////////////////////////////////////
355 void CP_State_RealParticlePlane::initComplete()
356 {
357  // If the sections are done for enl and zmat, then we are also done
358  // registering. The coalesced completion state is reported in to the
359  // global initialization phase process.
360  if( enlSectionComplete && planeRedSectionComplete)
361  {
362  // CkPrintf("RPP[%d %d] initComplete\n",thisIndex.x, thisIndex.y);
363  int constructed=1;
364  contribute(sizeof(int), &constructed, CkReduction::sum_int,
365  CkCallback(CkIndex_InstanceController::doneInit(NULL),
366  CkArrayIndex1D(thisInstance.proxyOffset),
367  instControllerProxy), thisInstance.proxyOffset);
368  }
369 }
370 
371 //////////////////////////////////////////////////////////////////////////////
372 /// The destructor : never called directly but I guess migration needs it
373 //////////////////////////////////////////////////////////////////////////////
374 //////////////////////////////////////////////////////////////////////////////
375 //////////////////////////////////////////////////////////////////////////////
376 CP_State_RealParticlePlane::~CP_State_RealParticlePlane(){
377 
378  if(ees_nonlocal==1){
379  fftw_free(projPsiC);
380  if(config.doublePack){
381  fftw_free(zmat);
382  fftw_free(zmatScr);
383  }else{
384  fftw_free(zmatC); // free the master
385  fftw_free(zmatScrC); // free the master
386  }//endif
387  }//endif
388 
389  projPsiC = NULL;
390  projPsiR = NULL;
391  zmat = NULL;
392  zmatScr = NULL;
393  zmatC = NULL;
394  zmatScrC = NULL;
395 
396 }
397 //////////////////////////////////////////////////////////////////////////////
398 
399 
400 //////////////////////////////////////////////////////////////////////////////
401 /// Recv FFT data/psi-projector from CP_StateParticlePlane
402 /// kicking things off for this Chare.
403 //////////////////////////////////////////////////////////////////////////////
404 //////////////////////////////////////////////////////////////////////////////
405 //////////////////////////////////////////////////////////////////////////////
406 void CP_State_RealParticlePlane::recvFromEesGPP(NLFFTMsg *msg){
407 //////////////////////////////////////////////////////////////////////////////
408 /// Make sure you belong here.
409 
410  if(ees_nonlocal==0){
411  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
412  CkPrintf("Yo dawg, ees nonlocal is off. RPP can't recvFromEesGPP\n");
413  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
414  CkExit();
415  }//endif
416 
417 //////////////////////////////////////////////////////////////////////////////
418 #ifdef _NAN_CHECK_
419  for(int i=0;i<msg->size ;i++)
420  {
421  CkAssert(finite(msg->data[i].re));
422  CkAssert(finite(msg->data[i].im));
423  }
424 #endif
425 //////////////////////////////////////////////////////////////////////////////
426 /// Unpack the message, get some local variables and increment counters
427 
428  int size = msg->size;
429  int iterNLNow = msg->step;
430  int Index = msg->senderIndex; // which g-space chare sent the data
431  complex *partiallyFFTd = msg->data;
432 
433  CPcharmParaInfo *sim = CPcharmParaInfo::get();
434  int nchareG = sim->nchareG;
435  int **tranUnpack = sim->index_tran_upackNL;
436  int *nline_per_chareG = sim->nlines_per_chareG;
437 
438  count++;
439  // if we are starting a totally new time step, increment the time; zero the nonlocal energy
440  if(iterNL==0 && count==1){itime++; cp_enl=0.0;}
441 
442  // If we doing a new iteration, refresh the memory, increment iterNL
443  // In real space, we must zero out here because not every value
444  // is initialized. e.g. some zero elements are not set.
445  if(count==1){
446  iterNL++;
447  bzero(projPsiR,planeSize*sizeof(double)); // OK when !doublepack too
448  }//endif
449 
450 //////////////////////////////////////////////////////////////////////////////
451 /// Consistency Checking
452 
453  if (count > nchareG) {
454  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
455  CkPrintf(
456  "Mismatch in allowed # of NL-gspace chare arrays sending to NLR-chare: %d %d %d %d\n",
457  count,nchareG,thisIndex.x,thisIndex.y);
458  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
459  CkExit();
460  }//endif
461 
462  if(size!=nline_per_chareG[Index]){
463  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
464  CkPrintf("Dude, size %d != %d for NL-Rchare %d %d from G-chare %d\n",
465  size,nline_per_chareG[Index],thisIndex.y,Index,Index);
466  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
467  CkExit();
468  }//endif
469 
470  if(iterNLNow!=iterNL || recvBlock==1){
471  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
472  CkPrintf("Dude, iter count %d != %d for NL-Rchare %d %d %d\n",iterNLNow,iterNL,
473  thisIndex.y,Index,recvBlock);
474  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
475  CkExit();
476  }//endif
477 
478 //////////////////////////////////////////////////////////////////////////////
479 /// Copy out the data set and delete the message
480 
481 /// You have received packed data (x,y) from processor sendIndex
482 /// Every real space chare receives the same x,y indicies.
483 /// For double pack, x=0,1,2,3,4 ... y= {-K ... K}
484 /// The x increase with processor number. The y are split.
485 /// The rundescriptor contains all we need to unpack the data.
486 /// For doublepack : nffty*run[i][j].x + run[i][j].y
487 /// we store this stuff in the convenient package
488 /// Pictorially a half cylinder is sent which is unpacked into
489 /// a half cube for easy FFTing. Y is the inner index.
490 
491  for(int i=0;i< size;i++){projPsiC[tranUnpack[Index][i]] = partiallyFFTd[i];}
492  delete msg;
493 
494 //////////////////////////////////////////////////////////////////////////////
495 /// When you have collected the full data set, continue
496 
497  if(count == nchareG) {
498  count = 0;
499  recvBlock = 1;
500 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
501  if(thisIndex.x==0)
502  CkPrintf("HI, I am rPP %d %d in recvfromGpp : %d\n",thisIndex.x,thisIndex.y,iterNL);
503 #endif
504  fftDataDone=true;
505  if(launchFFT){
506  thisProxy(thisIndex.x,thisIndex.y).FFTNLEesFwdR();
507  if(rhoRTime!=itime){CkPrintf("Badddd launchFFT.1\n");CkExit();}
508  }else{
509  if(iterNL!=1){CkPrintf("Badddd launchFFT.2\n");CkExit();}
510  }//endif
511  }//endif
512 
513 //----------------------------------------------------------------------------
514  }// end routine
515 //////////////////////////////////////////////////////////////////////////////
516 
517 
518 //////////////////////////////////////////////////////////////////////////////
519 /// Complete the Forward FFT of psi-projector sent from g-space chares
520 //////////////////////////////////////////////////////////////////////////////
521 //////////////////////////////////////////////////////////////////////////////
522 //////////////////////////////////////////////////////////////////////////////
523 void CP_State_RealParticlePlane::FFTNLEesFwdR(){
524 //////////////////////////////////////////////////////////////////////////////
525 
526  fftDataDone=false;
527  CPcharmParaInfo *sim = CPcharmParaInfo::get();
528  int nplane_x = sim->nplane_x;
529  if(!config.doublePack){nplane_x = (nplane_x+1)/2;} // convert to the doublePack definition
530 
531 //////////////////////////////////////////////////////////////////////////////
532 
533 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
534  if(thisIndex.x==0)
535  CkPrintf("HI, I am rPP %d %d in FFTNL : %d %d %d %d\n",
536  thisIndex.x,thisIndex.y,iterNL,ngridA,ngridB,nplane_x);
537 #endif
538  // This is actually in place, projPsiR and projPsiC use the same memory location
539 #if CMK_TRACE_ENABLED
540  double StartTime=CmiWallTimer();
541 #endif
542 
543 #ifdef _Glenn_DBG_KPT_
544  char junk[1000];
545  FILE *fp;
546  if(thisIndex.x==0 && thisIndex.y==0){
547  sprintf(junk,"rpp%d.%d_beforeFFT.out",thisIndex.x,thisIndex.y);
548  fp = fopen(junk,"w");
549  for(int i=0;i<ngridB;i++){
550  for(int j=i*ngridA,k=0;j<(i+1)*ngridA;j++,k++){
551  fprintf(fp,"%d %d %.10g %.10g\n",i,k,projPsiC[j].re,projPsiC[j].im);
552  }//endfor
553  }//endfor
554  fclose(fp);
555  }//endif
556 #endif
557 
558  UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch()->doNlFFTGtoR_Rchare(projPsiC,projPsiR,
559  nplane_x,ngridA,ngridB,myPlane);
560 #ifdef _Glenn_DBG_KPT_
561  if(thisIndex.x==0 && thisIndex.y==0){
562  sprintf(junk,"rpp%d.%d_afterFFT.out",thisIndex.x,thisIndex.y);
563  fp = fopen(junk,"w");
564  for(int i=0;i<ngridB;i++){
565  for(int j=i*ngridA,k=0;j<(i+1)*ngridA;j++,k++){
566  fprintf(fp,"%d %d %.10g %.10g\n",i,k,projPsiC[j].re,projPsiC[j].im);
567  }//endfor
568  }//endfor
569  fclose(fp);
570  }//endif
571 #endif
572 
573 #if CMK_TRACE_ENABLED
574  traceUserBracketEvent(doNlFFTGtoR_, StartTime, CmiWallTimer());
575 #endif
576 
577 //////////////////////////////////////////////////////////////////////////////
578 /// Lots of Ugly debugging stuff
579 
580 #ifdef _CP_GS_DUMP_VKS_
581  dumpMatrix("projPsiC",(double *)projPsiC, 1, csize*2,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
582 #endif
583 
584 #ifdef _CP_GS_DEBUG_COMPARE_VKS_
585  // at this point we're checking that we did the fft correctly
586  // which implicitly checks that we collated the input correctly
587 
588  if(savedprojpsiC==NULL){ // load it
589  savedprojpsiC= new complex[csize];
590  loadMatrix("projPsiC",(double *)savedprojpsiC, 1, csize*2,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
591  }//endif
592  for(int i=0;i<csize;i++){
593  if(fabs(projPsiC[i].re-savedprojpsiC[i].re)>0.0001){
594  fprintf(stderr, "RPP [%d,%d] %d element projpsi %.10g not %.10g\n",thisIndex.x, thisIndex.y,i,
595  projPsiC[i].re, savedprojpsiC[i].re);
596  }//endif
597  CkAssert(fabs(projPsiC[i].re-savedprojpsiC[i].re)<0.0001);
598  CkAssert(fabs(projPsiC[i].im-savedprojpsiC[i].im)<0.0001);
599  }//endfor
600 #endif
601 
602 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
603  if(thisIndex.x==0)
604  CkPrintf("HI, I am rPP %d %d in FFTNL.2 : %d %d %d %d\n",
605  thisIndex.x,thisIndex.y,iterNL,ngridA,ngridB,nplane_x);
606 #endif
607 
608 //////////////////////////////////////////////////////////////////////////////
609 /// You can't start computing until all the cache stuff has completed.
610 
611  if(registrationFlag==1){computeZmatEes();}
612 
613  if(registrationFlag==0 && iterNL!=1 || iterNL==0){
614  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
615  CkPrintf("Dude, iter count %d > 1 for NL-Rchare %d %d and no reg?\n",iterNL,
616  thisIndex.x,thisIndex.y);
617  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
618  CkExit();
619  }//endif
620 
621 //----------------------------------------------------------------------------
622  }//end routine
623 //////////////////////////////////////////////////////////////////////////////
624 
625 
626 //////////////////////////////////////////////////////////////////////////////
627 /// Compute the Zmat elements of this iteration : Spawn the section reduction
628 //////////////////////////////////////////////////////////////////////////////
629 //////////////////////////////////////////////////////////////////////////////
630 //////////////////////////////////////////////////////////////////////////////
631 void CP_State_RealParticlePlane::computeZmatEes(){
632 
633  CPcharmParaInfo *sim = CPcharmParaInfo::get();
634  int iterNL1 = iterNL-1; // silly C++ convention
635  int *nmem_zmat = sim->nmem_zmat; // zmat size now
636  int nZmat = nmem_zmat[iterNL1];
637  if(!config.doublePack){nZmat *=2;} //size in doubles of complex beasty
638 
639 //////////////////////////////////////////////////////////////////////////////
640 /// Check out your B-splines from the cache
641 
642  eesCache *eesData = UeesCacheProxy[thisInstance.proxyOffset].ckLocalBranch ();
643  int *allowedRppChares = eesData->allowedRppChares;
644  CkAssert(allowedRppChares[myPlane]==1);
645 
646  if(iterNL==1){eesData->queryCacheRPP(myPlane,itime,iterNL);}// query once a t-step
647 
648  int *plane_index = eesData->RppData[myPlane]->plane_index;
649  int **igrid = eesData->RppData[myPlane]->igrid;
650  double **mn = eesData->RppData[myPlane]->mn;
651 
652 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
653  if(thisIndex.x==0)
654  CkPrintf("HI, I am rPP %d %d in computeZmat : %d\n",thisIndex.x,thisIndex.y,iterNL);
655 #endif
656 
657 //////////////////////////////////////////////////////////////////////////////
658 /// Compute Zmat
659 
660 #if CMK_TRACE_ENABLED
661  double StartTime=CmiWallTimer();
662 #endif
663 
664  if(config.doublePack){
665  CPNONLOCAL::eesZmatRchare(projPsiR,iterNL,zmat,igrid,mn,
666  plane_index,thisIndex.x,myPlane);
667  }else{
668  CPNONLOCAL::eesZmatRchareC(projPsiC,iterNL,zmatC,igrid,mn,
669  plane_index,thisIndex.x,myPlane);
670  }//endif
671 
672 #if CMK_TRACE_ENABLED
673  traceUserBracketEvent(eesZmatR_, StartTime, CmiWallTimer());
674 #endif
675 
676 //////////////////////////////////////////////////////////////////////////////
677 /// Launch section reduction : If !doublePack reduction is still done as doubles
678 
679 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
680  if(thisIndex.x==0)
681  CkPrintf("HI, I am rPP %d %d in send to zmat-red : %d %d %d d\n",
682  thisIndex.x,thisIndex.y,iterNL,reductionPlaneNum,nZmat,zmatSizeMax);
683 #endif
684 
685 #define _FANCY_RED_METHOD_
686 #ifdef _FANCY_RED_METHOD_
687  CkMulticastMgr *mcastGrp = CProxy_CkMulticastMgr(mCastGrpId).ckLocalBranch();
688  CkCallback cb(CkIndex_CP_State_RealParticlePlane::recvZMatEes(NULL),
689  CkArrayIndex2D(thisIndex.x,reductionPlaneNum),
690  UrealParticlePlaneProxy[thisInstance.proxyOffset].ckGetArrayID());
691  mcastGrp->contribute((nZmat*sizeof(double)),zmat,sumFastDoubleType,
692  rPlaneRedCookie,cb, iterNL);
693 #else
694  thisProxy(thisIndex.x,reductionPlaneNum).recvZMatEesSimp(nZmat,zmat,
695  thisIndex.x,thisIndex.y,iterNL);
696 #endif
697 
698 
699 //----------------------------------------------------------------------------
700  }//end routine
701 //////////////////////////////////////////////////////////////////////////////
702 
703 
704 //////////////////////////////////////////////////////////////////////////////
705 /// A chare can be behind by 1 iteration only. This message can arrive before
706 /// this chare has sent its zmat. This chare must therefore recv in zmatScr
707 //////////////////////////////////////////////////////////////////////////////
708 void CP_State_RealParticlePlane::recvZMatEesSimp(int size, double *_zmat,
709  int state, int index,int iterNL_in){
710 //////////////////////////////////////////////////////////////////////////////
711 
712  CPcharmParaInfo *sim = CPcharmParaInfo::get();
713  int iterNL1 = iterNL_in-1; // silly C++ convention
714  int *nmem_zmat = sim->nmem_zmat; // zmat size now
715  int nZmat = nmem_zmat[iterNL1];
716  if(!config.doublePack){nZmat *= 2;}
717 
718 //////////////////////////////////////////////////////////////////////////////
719 /// check for errors
720 
721  if(size != nZmat){
722  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
723  CkPrintf("Dude, size != nzmat : %d : %d %d : %d %d\n",
724  iterNL,thisIndex.x,thisIndex.y,state,index);
725  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
726  CkExit();
727  }//endif
728  if(iterNL!=iterNL_in && iterNL!=iterNL_in-1){
729  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
730  CkPrintf("Dude, wrong iterNls : %d %d : %d %d : %d %d\n",
731  iterNL,iterNL_in,thisIndex.x,thisIndex.y,state,index);
732  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
733  CkExit();
734  }//endif
735  if(state != thisIndex.x){
736  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
737  CkPrintf("Dude, you sent to the wrong state : %d %d : %d %d : %d %d\n",
738  iterNL,iterNL_in,thisIndex.x,thisIndex.y,state,index);
739  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
740  CkExit();
741  }//endif
742 
743 //////////////////////////////////////////////////////////////////////////////
744 /// Recv the reduced zmatrix and chuck the message. Recv into scratch because
745 /// it could be this chare has not yet contributed and you would overwrite
746 /// critical data if you used zmat[]
747 
748  countZ++;
749  if(countZ==1){bzero(zmatScr,nZmat*sizeof(double));}
750  for(int i=0;i<nZmat;i++){zmatScr[i]+=_zmat[i];}
751 
752 //////////////////////////////////////////////////////////////////////////////
753 /// Bcast the puppy dog back out to your friends : Bow-Wow-Wow
754 
755  if(countZ==nChareR){
756  if(iterNL != iterNL){
757  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
758  CkPrintf("Wrong iterNls strange racing: %d %d : %d %d : %d %d\n",
759  iterNL,iterNL_in,thisIndex.x,thisIndex.y,state,index);
760  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
761  CkExit();
762  }//endif
763 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
764  CkPrintf("HI, I am rPP %d %d in recvZmatSimp blasting off to forces: %d %d\n",
765  thisIndex.x,thisIndex.y,iterNL,index);
766 #endif
767  countZ=0;
768  for(int i=0;i<nChareR;i++){
769  CompAtmForcMsg *rmsg = new (nZmat, 8*sizeof(int)) CompAtmForcMsg;
770  if(config.prioNLFFTMsg){
771  CkSetQueueing(rmsg, CK_QUEUEING_IFIFO);
772  *(int*)CkPriorityPtr(rmsg) = config.gsNLfftpriority+thisIndex.x+iterNL_in;
773  }//endif
774  rmsg->init(nZmat,zmatScr,iterNL_in);
775  thisProxy(thisIndex.x,i).computeAtmForcEes(rmsg);
776  }//endfor
777  }//endif
778 
779 //----------------------------------------------------------------------------
780  }//end routine
781 //////////////////////////////////////////////////////////////////////////////
782 
783 
784 //////////////////////////////////////////////////////////////////////////////
785 /// Reduction client of zmat : Here everyone must have the same iteration
786 //////////////////////////////////////////////////////////////////////////////
787 //////////////////////////////////////////////////////////////////////////////
788 //////////////////////////////////////////////////////////////////////////////
789 void CP_State_RealParticlePlane::recvZMatEes(CkReductionMsg *msg){
790 //////////////////////////////////////////////////////////////////////////////
791 
792  CPcharmParaInfo *sim = CPcharmParaInfo::get();
793  int iterNL_in = msg->getUserFlag();
794  // int iterNL1 = iterNL_in-1; // silly C++ convention
795  int iterNL1 = iterNL-1; // silly C++ convention
796  int *nmem_zmat = sim->nmem_zmat; // zmat size now
797  int nZmat = nmem_zmat[iterNL1];
798  if(!config.doublePack){nZmat *=2;}
799 
800 //////////////////////////////////////////////////////////////////////////////
801 /// Some debuggin
802 
803 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
804  if(thisIndex.x==0)
805  CkPrintf("HI, I am rPP %d %d in recvZmat : %d\n",thisIndex.x,thisIndex.y,iterNL);
806 #endif
807 #ifdef _NAN_CHECK_
808  for(int i=0;i<msg->getSize()/sizeof(double) ;i++){
809  CkAssert(finite(((double*) msg->getData())[i]));
810  }//endfor
811 #endif
812 
813 //////////////////////////////////////////////////////////////////////////////
814 /// Recv the reduced zmatrix into zmat and chuck the message
815 
816  CkAssert(msg->getSize() == nZmat* sizeof(double));
817  double *realValues = (double *) msg->getData();
818  for(int i=0;i<nZmat;i++){zmat[i]=realValues[i];}
819 
820  delete msg;
821 
822 //////////////////////////////////////////////////////////////////////////////
823 /// Bcast the puppy dog back out to your friends : Bow-Wow-Wow
824 
825 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
826  if(thisIndex.x==0)
827  CkPrintf("HI, I am rPP %d %d in recvZmat sending to compute: %d\n",
828  thisIndex.x,thisIndex.y,iterNL);
829 #endif
830  CompAtmForcMsg *rmsg=new (nZmat, 8*sizeof(int)) CompAtmForcMsg;
831  if(config.prioNLFFTMsg){
832  CkSetQueueing(rmsg, CK_QUEUEING_IFIFO);
833  *(int*)CkPriorityPtr(rmsg) = config.gsNLfftpriority+thisIndex.x+iterNL_in;
834  }//endif
835 
836  rmsg->init(nZmat,zmat,iterNL);
837  CkAssert(rmsg->iterNL>0);
838  rPlaneSectProxy.computeAtmForcEes(rmsg);
839 
840 //----------------------------------------------------------------------------
841  }//end routine
842 //////////////////////////////////////////////////////////////////////////////
843 
844 
845 //////////////////////////////////////////////////////////////////////////////
846 /// Use Zmat and ProjPsi to get atmForces, Energy and psiforces
847 //////////////////////////////////////////////////////////////////////////////
848 //////////////////////////////////////////////////////////////////////////////
849 //////////////////////////////////////////////////////////////////////////////
850 void CP_State_RealParticlePlane::computeAtmForcEes(CompAtmForcMsg *msg){
851 //////////////////////////////////////////////////////////////////////////////
852 /// Unpack the message : Get some sizes : copy out message.
853 /// : Do not delete msg. Its a nokeep!!!!!!!!!!!!
854 //////////////////////////////////////////////////////////////////////////////
855 
856 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
857  if(thisIndex.x==0)
858  CkPrintf("HI, I am rPP %d %d in compteAtmforc : %d\n",thisIndex.x,thisIndex.y,iterNL);
859 #endif
860 
861  int nZmat_in=msg->nZmat;
862 
863  double *zmat_loc=msg->zmat;
864  int iterNL_in=msg->iterNL;
865  CkAssert(msg->iterNL>0);
866  CPcharmParaInfo *sim = CPcharmParaInfo::get();
867  int iterNL1 = iterNL-1; // silly C++ convention
868  int *nmem_zmat = sim->nmem_zmat; // zmat memory size
869  int nZmat = nmem_zmat[iterNL1];
870  if(!config.doublePack){nZmat *=2;} //size in doubles of complex beasty
871 
872  /* In the reduction case each iteration arrives together. This coordination scheme is redundant */
873  if(iterNL_in != iterNL || nZmat_in != nZmat){
874  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
875  CkPrintf("RPP [%d,%d] Dude, iteration mismatch : %d %d z %d %d\n",
876  thisIndex.x,thisIndex.y,
877  iterNL,iterNL_in,nZmat_in, nZmat );
878  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
879  CkExit();
880  }//endif
881 
882  CmiMemcpy(zmat,zmat_loc,sizeof(double)*nZmat_in);
883 
884 //////////////////////////////////////////////////////////////////////////////
885 /// Check out your B-splines from the cache and get some parameters
886 
887  eesCache *eesData = UeesCacheProxy[thisInstance.proxyOffset].ckLocalBranch ();
888  FFTcache *fftcache = UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
889 
890  int *allowedRppChares = eesData->allowedRppChares;
891  CkAssert(allowedRppChares[myPlane]==1);
892 
893  int *plane_index = eesData->RppData[myPlane]->plane_index;
894  int **igrid = eesData->RppData[myPlane]->igrid;
895  int *nBreakJ = eesData->RppData[myPlane]->nBreakJ;
896  int **sBreakJ = eesData->RppData[myPlane]->sBreakJ;
897  double **dmn_x = eesData->RppData[myPlane]->dmn_x;
898  double **dmn_y = eesData->RppData[myPlane]->dmn_y;
899  double **dmn_z = eesData->RppData[myPlane]->dmn_z;
900  double **mn = eesData->RppData[myPlane]->mn;
901  double *projPsiRScr = fftcache->tmpDataR;
902  complex *projPsiCScr = fftcache->tmpData;
903  fftcache->getCacheMem("CP_State_RealParticlePlane::computeAtmForcEes");
904 
905  AtomsCache *ag = UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
906  FastAtoms *fastAtoms = &(ag->fastAtoms);
907 
908  // projPsiR comes in with info for atoms
909  // projPsiRSsr leaves with info for psiforces to (ngrid_a+2)*ngrid_b elements
910  int n_a, n_b, n_c, n_interp, nAtm;
911  CPNONLOCAL::getEesPrms(&n_a,&n_b,&n_c,&n_interp,&nAtm);
912  int n_interp2=n_interp*n_interp;
913 
914 //////////////////////////////////////////////////////////////////////////////
915 /// Debug yourself
916 
917 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
918  if(thisIndex.x==0)
919  CkPrintf("HI, I am rPP %d %d in compteAtmforc : %d\n",thisIndex.x,thisIndex.y,iterNL);
920 #endif
921 
922 #ifdef _CP_GS_DUMP_VKS_
923 #ifdef _INSANEO_PARANOID_COMPARE_THAT_EATS_HUGE_MEMORY_
924  dumpMatrix2DDouble("mn",mn, nAtm, n_interp2, thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
925  dumpMatrix2DDouble("dmn_x",dmn_x, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
926  dumpMatrix2DDouble("dmn_y",dmn_y, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
927 
928  dumpMatrix2DDouble("dmn_z",dmn_z, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
929  dumpMatrix2DInt("igrid",igrid, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
930 #endif
931 #endif
932 
933 #ifdef _CP_GS_DEBUG_COMPARE_VKS_
934  // at this point we're checking that nobody poisoned projPsiR, igrid dmn_[xyz] mn
935  // projPsiR uses the same location as projpsiC, so we hold our nose
936  // and cut and paste
937  if(savedprojpsiC==NULL)
938  { // load it
939  savedprojpsiC= new complex[csize];
940  loadMatrix("projPsiC",(double *)savedprojpsiC, 1, csize*2,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
941  }
942 #ifdef _INSANEO_PARANOID_COMPARE_THAT_EATS_HUGE_MEMORY_
943  if(savedmn==NULL){ // load it
944  savedigrid = (int **)fftw_malloc(nAtm*sizeof(int*));
945 
946  savedmn = (double **)fftw_malloc(nAtm*sizeof(double*));
947  saveddmn_x = (double **)fftw_malloc(nAtm*sizeof(double*));
948  saveddmn_y = (double **)fftw_malloc(nAtm*sizeof(double*));
949  saveddmn_z = (double **)fftw_malloc(nAtm*sizeof(double*));
950  // Argh! this crazy iterate from 1 stuff escaped from
951  // the piny box of evil fortranisms
952  for(int i=0;i<nAtm;i++){
953  savedigrid[i] = (int *)fftw_malloc(n_interp2*sizeof(int))-1;
954  double *tmp = (double *)fftw_malloc(4*n_interp2*sizeof(double));
955  // by dint of more ugly trickery we store 4 array rows in 1
956  // which provides a nice way for us to quietly violate logical
957  // array boundaries without tripping over memory we don't own
958  int ioff = 0;
959  savedmn[i] = &tmp[ioff]-1; ioff+=n_interp2;
960  saveddmn_x[i] = &tmp[ioff]-1; ioff+=n_interp2;
961  saveddmn_y[i] = &tmp[ioff]-1; ioff+=n_interp2;
962  saveddmn_z[i] = &tmp[ioff]-1;
963  }//endfor
964 
965  loadMatrix2DDouble("mn",savedmn, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
966  loadMatrix2DDouble("dmn_x",saveddmn_x, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
967  loadMatrix2DDouble("dmn_y",saveddmn_y, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
968  loadMatrix2DDouble("dmn_z",saveddmn_z, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
969  loadMatrix2DInt("igrid",savedigrid, nAtm, n_interp2,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
970  }//endif : loading
971 
972  for(int i=0;i<nAtm;i++){
973  if(plane_index[i]!=0){
974  for(int j=1;j<=n_interp2;j++){
975  if(fabs(mn[i][j]-savedmn[i][j])>0.0001){
976  fprintf(stderr,"RPP [%d,%d] %d %d element mn %.10g not %.10g\n",thisIndex.x,thisIndex.y,i,j,mn[i][j],savedmn[i][j]);
977  }//endif
978  CkAssert(fabs(mn[i][j]-savedmn[i][j])<0.0001);
979  if(fabs(dmn_x[i][j]-saveddmn_x[i][j])>0.0001){
980  fprintf(stderr,"RPP [%d,%d] %d %d element dmn_x %.10g not %.10g\n",thisIndex.x,thisIndex.y,i,j,dmn_x[i][j],saveddmn_x[i][j]);
981  }//endif
982  CkAssert(fabs(dmn_x[i][j]-saveddmn_x[i][j])<0.0001);
983  if(fabs(dmn_y[i][j]-saveddmn_y[i][j])>0.0001){
984  fprintf(stderr,"RPP [%d,%d] %d %d element dmn_y %.10g not %.10g\n",thisIndex.x,thisIndex.y,i,j,dmn_y[i][j],saveddmn_y[i][j]);
985  }//endif
986  CkAssert(fabs(dmn_y[i][j]-saveddmn_y[i][j])<0.0001);
987  if(fabs(dmn_z[i][j]-saveddmn_z[i][j])>0.0001){
988  fprintf(stderr,"RPP [%d,%d] %d %d element dmn_z %.10g not %.10g\n",thisIndex.x,thisIndex.y,i,j,dmn_z[i][j],saveddmn_z[i][j]);
989  }//endif
990  CkAssert(fabs(dmn_z[i][j]-saveddmn_z[i][j])<0.0001);
991  if(igrid[i][j]!=savedigrid[i][j]){
992  fprintf(stderr,"RPP [%d,%d] %d %d element igrid %d not %d\n",thisIndex.x,thisIndex.y,i,j,igrid[i][j],savedigrid[i][j]);
993  }//endif
994  CkAssert(igrid[i][j]==savedigrid[i][j]);
995  }//end for
996  }//end if
997  }//endfor
998 #endif
999 #endif
1000 
1001 //////////////////////////////////////////////////////////////////////////////
1002 /// Perform the computation of the atom forces and energy finally
1003 
1004 #if CMK_TRACE_ENABLED
1005  double StartTime=CmiWallTimer();
1006 #endif
1007 
1008  double cp_enl_now = 0.0;
1009  if(config.doublePack){
1010  CPNONLOCAL::eesEnergyAtmForcRchare(iterNL,&cp_enl_now,zmat,igrid,mn,dmn_x,dmn_y,dmn_z,
1011  projPsiR,projPsiRScr,plane_index,nBreakJ,sBreakJ,
1012  myPlane,thisIndex.x,kpoint_ind,fastAtoms);
1013  }else{
1014  CPNONLOCAL::eesEnergyAtmForcRchareC(iterNL,&cp_enl_now,zmatC,igrid,mn,dmn_x,dmn_y,dmn_z,
1015  projPsiC,projPsiCScr,plane_index,nBreakJ,sBreakJ,
1016  myPlane,thisIndex.x,kpoint_ind,fastAtoms);
1017  }//endif
1018  cp_enl += cp_enl_now;
1019 
1020 #if CMK_TRACE_ENABLED
1021  traceUserBracketEvent(eesEnergyAtmForcR_, StartTime, CmiWallTimer());
1022 #endif
1023 
1024 //////////////////////////////////////////////////////////////////////////////
1025 /// If we are done, send out the energy : HELP HELP Evil Section Multicast
1026 
1027  if(thisIndex.y==reductionPlaneNum && iterNL==numIterNl){
1028 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
1029  if(thisIndex.x==0)
1030  CkPrintf("HI, I am rPP %d %d sending Enl : %d\n",thisIndex.x,thisIndex.y,iterNL);
1031 #endif
1032 
1033 #ifdef _FANCY_RED_METHOD_
1034  CkMulticastMgr *mcastGrp = CProxy_CkMulticastMgr(mCastGrpId).ckLocalBranch();
1035  // CkCallback cb=CkCallback(printEnlR, NULL);
1036  CkCallback cb(CkIndex_CP_State_RealParticlePlane::printEnlR(NULL),
1037  CkArrayIndex2D(0,reductionPlaneNum),
1038  UrealParticlePlaneProxy[thisInstance.proxyOffset].ckGetArrayID());
1039  mcastGrp->contribute(sizeof(double),(void*) &cp_enl,
1040  CkReduction::sum_double,rEnlCookie, cb, itime);
1041 #else
1042  thisProxy(0,0).printEnlRSimp(cp_enl,thisIndex.x,itime);
1043 #endif
1044  }//endif : we are done and can reduce the non-local energy
1045 
1046 //////////////////////////////////////////////////////////////////////////////
1047 /// More nasty debuggin
1048 
1049 #ifdef _CP_GS_DUMP_VKS_
1050  dumpMatrix("zmat",(double *)zmat, 1, nZmat_in,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
1051  dumpMatrix("projPsiRScr",(double *)projPsiRScr, 1, planeSize,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
1052 #endif
1053 
1054 #ifdef _CP_GS_DEBUG_COMPARE_VKS_
1055  // check that we did the zmat gather/scatter correctly and
1056 
1057  if(savedzmat==NULL)
1058  { // load it
1059  savedzmat= new double[nZmat_in];
1060  loadMatrix("zmat",(double *)savedzmat, 1, nZmat_in,thisIndex.y,thisIndex.x,thisIndex.x,iterNL,false);
1061  }
1062  for(int i=0;i<nZmat_in; i++)
1063  {
1064  if(fabs(zmat[i]-savedzmat[i])>0.0001)
1065  {
1066  fprintf(stderr, "RPP [%d,%d] %d element projpsi %.10g not %.10g\n",thisIndex.x, thisIndex.y,i, zmat[i], savedzmat[i]);
1067  }
1068  CkAssert(fabs(zmat[i]-savedzmat[i])<0.0001);
1069 
1070  }
1071  if(savedProjpsiRScr==NULL)
1072  { // load it
1073  savedProjpsiRScr= new double[planeSize];
1074  loadMatrix("projPsiRScr",(double *)savedProjpsiRScr, 1, planeSize,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
1075  }
1076  for(int i=0;i<planeSize;i++)
1077  {
1078  if(fabs(projPsiRScr[i]-savedProjpsiRScr[i])>0.0001)
1079  {
1080  fprintf(stderr, "RPP [%d,%d] %d element projpsi %.10g not %.10g\n",thisIndex.x, thisIndex.y,i, projPsiRScr[i], savedProjpsiRScr[i]);
1081  dumpMatrix("badprojPsiRScr",(double *)projPsiRScr, 1, planeSize,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
1082  }
1083  CkAssert(fabs(projPsiRScr[i]-savedProjpsiRScr[i])<0.0001);
1084 
1085  }
1086 
1087 #endif
1088 
1089 //////////////////////////////////////////////////////////////////////////////
1090 /// zero the total enl energy if we are done.
1091 
1092  if(iterNL==numIterNl){
1093  cp_enl = 0.0;
1094  launchFFT=false;
1095  }//endif
1096 
1097 //////////////////////////////////////////////////////////////////////////////
1098 /// Time to make the Psiforces (donuts!)
1099 
1100  FFTNLEesBckR();
1101 
1102 //////////////////////////////////////////////////////////////////////////////
1103 /// Do not delete the nokeep message!!!
1104 //----------------------------------------------------------------------------
1105  }//end routine
1106 //////////////////////////////////////////////////////////////////////////////
1107 
1108 //////////////////////////////////////////////////////////////////////////////
1109 /// Control launch of FFT based on Rho having more of its act together
1110 //////////////////////////////////////////////////////////////////////////////
1111 //////////////////////////////////////////////////////////////////////////////
1112 //////////////////////////////////////////////////////////////////////////////
1113 void CP_State_RealParticlePlane::launchFFTControl(int time_in){
1114  rhoRTime = time_in;
1115  launchFFT=true;
1116  if(fftDataDone){
1117  if(iterNL!=1){CkPrintf("Badddd launchFFT.3\n");CkExit();}
1118  if(rhoRTime!=itime){CkPrintf("Badddd launchFFT.4\n");CkExit();}
1119  thisProxy(thisIndex.x,thisIndex.y).FFTNLEesFwdR();
1120  }//endif
1121 //----------------------------------------------------------------------------
1122 }
1123 //////////////////////////////////////////////////////////////////////////////
1124 
1125 
1126 //////////////////////////////////////////////////////////////////////////////
1127 /// Do the FFT of projPsi(x,y,z) to get projPsi(gx,gy,z)
1128 //////////////////////////////////////////////////////////////////////////////
1129 //////////////////////////////////////////////////////////////////////////////
1130 //////////////////////////////////////////////////////////////////////////////
1131 void CP_State_RealParticlePlane::FFTNLEesBckR(){
1132 //////////////////////////////////////////////////////////////////////////////
1133 
1134  CPcharmParaInfo *sim = CPcharmParaInfo::get();
1135  int nplane_x = sim->nplane_x;
1136  if(!config.doublePack){nplane_x = (nplane_x+1)/2;} // convert to the doublePack definition
1137 
1138 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
1139  if(thisIndex.x==0)
1140  CkPrintf("HI, I am rPP %d %d in FFTNLbck : %d %d %d %d\n",
1141  thisIndex.x,thisIndex.y,iterNL,ngridA,ngridB,nplane_x);
1142 #endif
1143 
1144  FFTcache *fftcache = UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
1145  double *projPsiRScr = fftcache->tmpDataR;
1146  complex *projPsiCScr = fftcache->tmpData;
1147 
1148 //////////////////////////////////////////////////////////////////////////////
1149 /// Debugging stuff
1150 
1151 #ifdef _CP_GS_DEBUG_COMPARE_VKS_
1152  // double check that projPsiRScr is still ok
1153  if(savedProjpsiRScr==NULL){ // load it
1154  savedProjpsiRScr= new double[planeSize];
1155  loadMatrix("projPsiRScr",(double *)savedProjpsiRScr, 1, planeSize,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
1156  }//endif
1157 
1158  for(int i=0;i<planeSize;i++){
1159  if(fabs(projPsiRScr[i]-savedProjpsiRScr[i])>0.0001){
1160  fprintf(stderr,"RPP[%d,%d] %d element projpsi %.10g not %.10g\n",thisIndex.x,thisIndex.y,i,projPsiRScr[i],savedProjpsiRScr[i]);
1161  }//endif
1162  CkAssert(fabs(projPsiRScr[i]-savedProjpsiRScr[i])<0.0001);
1163  CkAssert(fabs(projPsiRScr[i]-savedProjpsiRScr[i])<0.0001);
1164  }//endfor
1165 #endif
1166 
1167 //////////////////////////////////////////////////////////////////////////////
1168 /// here is the computation of the fft
1169 
1170 #if CMK_TRACE_ENABLED
1171  double StartTime= CmiWallTimer();
1172 #endif
1173 
1174  fftcache->doNlFFTRtoG_Rchare(projPsiCScr,projPsiRScr,nplane_x,ngridA,ngridB,myPlane);
1175 
1176 #ifndef CMK_OPTIMIZE
1177  traceUserBracketEvent(doNlFFTRtoG_, StartTime, CmiWallTimer());
1178 #endif
1179 
1180 //////////////////////////////////////////////////////////////////////////////
1181 /// Here is the tedious debuggins
1182 
1183 #ifdef _CP_GS_DUMP_VKS_
1184  dumpMatrix("projPsiCScr",(double *)projPsiCScr, 1, (ngridA+2)*ngridB,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
1185 #endif
1186 
1187 #ifdef _CP_GS_DEBUG_COMPARE_VKS_
1188  // at this point we're checking that we did the fft correctly
1189  // which implicitly checks that we did the zmat gather/scatter correctly
1190 
1191  if(savedProjpsiCScr==NULL){ // load it
1192  savedProjpsiCScr= new complex[csize];
1193  loadMatrix("projPsiCScr",(double *)savedProjpsiCScr, 1, csize*2,thisIndex.y,thisIndex.x,thisIndex.x,0,false);
1194  }//endif
1195  for(int i=0;i<csize;i++){
1196  if(fabs(projPsiCScr[i].re-savedProjpsiCScr[i].re)>0.0001){
1197  fprintf(stderr,"RPP [%d,%d] %d element projpsi %.10g not %.10g\n",thisIndex.x,thisIndex.y,i,
1198  projPsiCScr[i].re,savedProjpsiCScr[i].re);
1199  fprintf(stderr, "HI, I am rPP %d %d in FFTNLbck : %d %d %d %d\n",
1200  thisIndex.x,thisIndex.y,iterNL,ngridA,ngridB,nplane_x);
1201  }//endif
1202  CkAssert(fabs(projPsiCScr[i].re-savedProjpsiCScr[i].re)<0.0001);
1203  CkAssert(fabs(projPsiCScr[i].im-savedProjpsiCScr[i].im)<0.0001);
1204  }//endfor
1205 #endif
1206 
1207 #if CMK_TRACE_ENABLED
1208  traceUserBracketEvent(doNlFFTRtoG_, StartTime, CmiWallTimer());
1209 #endif
1210 //////////////////////////////////////////////////////////////////////////////
1211 /// Now that the FFT is done, go back to Gspace
1212 
1213  sendToEesGPP();
1214 
1215 //////////////////////////////////////////////////////////////////////////////
1216  }//end routine
1217 //////////////////////////////////////////////////////////////////////////////
1218 
1219 
1220 //////////////////////////////////////////////////////////////////////////////
1221 /// Send the PsiForce NL FFT back to GparticlePlane.
1222 /// Until ALL RealPP chares send forces, PP cannot start the next iteration
1223 //////////////////////////////////////////////////////////////////////////////
1224 //////////////////////////////////////////////////////////////////////////////
1225 //////////////////////////////////////////////////////////////////////////////
1226 void CP_State_RealParticlePlane::sendToEesGPP(){
1227 //////////////////////////////////////////////////////////////////////////////
1228 
1229  CPcharmParaInfo *sim = CPcharmParaInfo::get();
1230  FFTcache *fftcache = UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
1231  int nchareG = sim->nchareG;
1232  int **tranpack = sim->index_tran_upackNL;
1233  int *nlines_per_chareG = sim->nlines_per_chareG;
1234  complex *projPsiCScr = fftcache->tmpData;
1235 
1236 ///////////////////////////////////////////////////////////////////==
1237 /// Perform the transpose and then the blast off the final 1D-FFT
1238 
1239 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
1240  if(thisIndex.x==0)
1241  CkPrintf("HI, I am rPP %d %d in sendtoGPP : %d\n",thisIndex.x,thisIndex.y,iterNL);
1242 #endif
1243 
1244 #ifdef USE_COMLIB
1245 #ifdef OLD_COMMLIB
1246  if(config.useMssInsGPP){mssPInstance.beginIteration();}
1247 #else
1248 /// if(config.useMssInsGPP){ComlibBegin(UparticlePlaneProxy[thisInstance.proxyOffset], iterNL);}
1249 #endif
1250 #endif
1251 
1252  for (int ic=0;ic<nchareG;ic++) { // chare arrays to which we will send
1253  int sendFFTDataSize = nlines_per_chareG[ic];
1254  GSPPIFFTMsg *msg = new (sendFFTDataSize, 8 * sizeof(int)) GSPPIFFTMsg;
1255  msg->iterNL = iterNL;
1256  msg->size = sendFFTDataSize;
1257  msg->offset = thisIndex.y; // c-plane-index
1258  complex *data = msg->data;
1259  if(config.prioNLFFTMsg){
1260  CkSetQueueing(msg, CK_QUEUEING_IFIFO);
1261  *(int*)CkPriorityPtr(msg) = config.gsNLfftpriority+thisIndex.x;
1262  }//endif
1263  for(int i=0;i<sendFFTDataSize;i++){data[i] = projPsiCScr[tranpack[ic][i]];}
1264  UparticlePlaneProxy[thisInstance.proxyOffset](thisIndex.x, ic).recvFromEesRPP(msg); // send the message
1265  }//end for : chare sending
1266 
1267 #ifdef USE_COMLIB
1268 #ifdef OLD_COMMLIB
1269  if(config.useMssInsGPP){mssPInstance.endIteration();}
1270 #else
1271 /// if(config.useMssInsGPP){ComlibEnd(UparticlePlaneProxy[thisInstance.proxyOffset], iterNL);}
1272 #endif
1273 #endif
1274 
1275 //////////////////////////////////////////////////////////////////////////////
1276 /// If it looks like this is the end, reset my counters baby.
1277 /// If its not the end, GParticlePlane will invoke the entry methods again
1278 
1279  recvBlock = 0; // I am done, I can recv messages now from PP
1280  // If I recv with recvblock==1, something is terribly wrong
1281  if(iterNL==numIterNl){
1282  iterNL = 0;
1283  cp_enl = 0;
1284  launchFFT=false;
1285  }//endif
1286 
1287  fftcache->freeCacheMem("CP_State_RealParticlePlane::sendToEesGPP");
1288 
1289 //----------------------------------------------------------------------------
1290  }//end routine
1291 //////////////////////////////////////////////////////////////////////////////
1292 
1293 
1294 //////////////////////////////////////////////////////////////////////////////
1295 /**
1296  * spread the reduction plane numbers around to minimize map collisions
1297  */
1298 //////////////////////////////////////////////////////////////////////////////
1299 //////////////////////////////////////////////////////////////////////////////
1300 //////////////////////////////////////////////////////////////////////////////
1301 int CP_State_RealParticlePlane::calcReductionPlaneNum(int state){
1302 //////////////////////////////////////////////////////////////////////////////
1303 
1304  int nstatemax=nstates-1;
1305  int ncharemax=nChareR-1;
1306  int planeNum= (state %ncharemax);
1307  if(planeNum<0){
1308  CkPrintf(" PP [%d %d] calc nstatemax %d ncharemax %d state %d planenum %d\n",
1309  thisIndex.x, thisIndex.y,nstatemax, ncharemax, state, planeNum);
1310  CkExit();
1311  }//endif
1312  return planeNum;
1313 
1314 }//end routine
1315 //////////////////////////////////////////////////////////////////////////////
1316 
1317 
1318 //////////////////////////////////////////////////////////////////////////////
1319 /// Recv a dummy message, 1 integer, and set the cookie monster
1320 //////////////////////////////////////////////////////////////////////////////
1321 //////////////////////////////////////////////////////////////////////////////
1322 //////////////////////////////////////////////////////////////////////////////
1323 void CP_State_RealParticlePlane::setPlaneRedCookie(EnlCookieMsg *m){
1324 //////////////////////////////////////////////////////////////////////////////
1325 /// Do not delete msg. Its a nokeep.
1326 //////////////////////////////////////////////////////////////////////////////
1327 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
1328  if(thisIndex.x==0)
1329  CkPrintf("RPP[%d %d] gets rPlaneRedCookie\n",thisIndex.x, thisIndex.y);
1330 #endif
1331 
1332  CkGetSectionInfo(rPlaneRedCookie,m);
1333  int cookiedone=1;
1334  CkMulticastMgr *mcastGrp = CProxy_CkMulticastMgr(mCastGrpId).ckLocalBranch();
1335  CkCallback cb(CkIndex_CP_State_RealParticlePlane::planeRedSectDone(NULL),
1336  CkArrayIndex2D(thisIndex.x,reductionPlaneNum),
1337  UrealParticlePlaneProxy[thisInstance.proxyOffset].ckGetArrayID());
1338  mcastGrp->contribute(sizeof(int),&cookiedone,CkReduction::sum_int,
1339  rPlaneRedCookie,cb);
1340  // if not root, we are set up
1341  if(thisIndex.y!=reductionPlaneNum)
1342  {
1343  planeRedSectionComplete=true;
1344  initComplete();
1345  }
1346 }
1347 //////////////////////////////////////////////////////////////////////////////
1348 
1349 
1350 //////////////////////////////////////////////////////////////////////////////
1351 /// Recv a dummy message, 1 integer, and set the cookie monster
1352 //////////////////////////////////////////////////////////////////////////////
1353 //////////////////////////////////////////////////////////////////////////////
1354 //////////////////////////////////////////////////////////////////////////////
1355 void CP_State_RealParticlePlane::setEnlCookie(EnlCookieMsg *m){
1356 //////////////////////////////////////////////////////////////////////////////
1357 /// Do not delete msg. Its a nokeep.
1358 //////////////////////////////////////////////////////////////////////////////
1359 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
1360  if(thisIndex.x==0)
1361  CkPrintf("RPP[%d %d] gets rEnlCookie\n",thisIndex.x, thisIndex.y);
1362 #endif
1363  CkGetSectionInfo(rEnlCookie,m);
1364  int cookiedone=1;
1365  CkMulticastMgr *mcastGrp = CProxy_CkMulticastMgr(mCastGrpId).ckLocalBranch();
1366  CkCallback cb(CkIndex_CP_State_RealParticlePlane::enlSectDone(NULL),
1367  CkArrayIndex2D(0,reductionPlaneNum),
1368  UrealParticlePlaneProxy[thisInstance.proxyOffset].ckGetArrayID());
1369  mcastGrp->contribute(sizeof(int), &cookiedone,
1370  CkReduction::sum_int,rEnlCookie, cb);
1371  // if not root, we are set up
1372  if(thisIndex.y!=reductionPlaneNum || thisIndex.x!=0){
1373  enlSectionComplete=true;
1374  initComplete();
1375  }
1376 }
1377 //////////////////////////////////////////////////////////////////////////////
1378 
1379 ///////////////////////////////////////////////////////////////////////////=
1380 /// Make sure everyone is registered on the 1st time step
1381 ///////////////////////////////////////////////////////////////////////////=
1382 ///////////////////////////////////////////////////////////////////////////c
1383 ///////////////////////////////////////////////////////////////////////////=
1384  void CP_State_RealParticlePlane::registrationDone() {
1385 ///////////////////////////////////////////////////////////////////////////=
1386 
1387 
1388 #ifdef _CP_DEBUG_STATE_RPP_VERBOSE_
1389  if(thisIndex.x==0)
1390  CkPrintf("HI, I am realPPa %d %d \n",thisIndex.x,thisIndex.y);
1391 #endif
1392 
1393  registrationFlag=1;
1394  if(iterNL==1){computeZmatEes();}
1395 
1396  if(iterNL>1){
1397  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1398  CkPrintf("Homes, registeration must occur before the 1st time step in RPP\n");
1399  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1400  CkExit();
1401  }//endif
1402 
1403  }//end routine
1404 ///////////////////////////////////////////////////////////////////////////=
1405 
1406 
1407 //////////////////////////////////////////////////////////////////////////////
1408 //////////////////////////////////////////////////////////////////////////////
1409 //////////////////////////////////////////////////////////////////////////////
1410 void CP_State_RealParticlePlane::pup(PUP::er &p) {
1411 //////////////////////////////////////////////////////////////////////////////
1412 
1413  p|rhoRTime;
1414  p|ees_nonlocal;
1415  p|nChareR;
1416  p|nChareG;
1417  p|Rstates_per_pe;
1418  p|myPlane;
1419  p|registrationFlag;
1420  p|recvBlock;
1421 
1422  p|ibead_ind; p|kpoint_ind; p|itemper_ind;
1423 
1424  p|numIterNl;
1425  p|countEnl;
1426  p|count;
1427  p|iterNL;
1428  p|itime;
1429  p|itimeRed;
1430  p|countZ;
1431 
1432  p|ngridA;
1433  p|ngridB;
1434  p|ngridC;
1435  p|planeSize;
1436  p|planeSizeT;
1437  p|csize;
1438  p|zmatSizeMax;
1439  p|reductionPlaneNum;
1440  p|cp_enl;
1441  p|cp_enlTot;
1442 
1443  p|rPlaneSectProxy;
1444  p|rPlaneENLProxy;
1445  p|rPlaneRedCookie;
1446  p|rEnlCookie;
1447  p|gPP_proxy;
1448 
1449  if(ees_nonlocal==1){
1450  int nZ = zmatSizeMax;
1451  if(config.doublePack){nZ *=2;}
1452  if (p.isUnpacking()) {
1453  projPsiC = (complex*) fftw_malloc(csize*sizeof(complex));
1454  projPsiR = reinterpret_cast<double*> (projPsiC);
1455  if(config.doublePack){
1456  zmat = (double*) fftw_malloc(zmatSizeMax*sizeof(double));
1457  zmatScr = (double*) fftw_malloc(zmatSizeMax*sizeof(double));
1458  }else{
1459  zmatC = (complex*) fftw_malloc(zmatSizeMax*sizeof(complex));
1460  zmatScrC = (complex*) fftw_malloc(zmatSizeMax*sizeof(complex));
1461  zmat = reinterpret_cast<double*> (zmatC);
1462  zmatScr = reinterpret_cast<double*> (zmatScrC);
1463  }//endif
1464  }//endif
1465  p(zmat,nZ); // pup as double even when zmat is complex
1466  p(zmatScr,nZ); // pup as double even when zmat is complex
1467  p((char*)projPsiC,csize *sizeof(complex)); // Why is there a cast to char ?
1468  }//endif
1469 
1470 //---------------------------------------------------------------------------
1471  }//end routine
1472 //////////////////////////////////////////////////////////////////////////////
1473 /*@}*/
CP_State_RealParticlePlane::recvZMatEes
void recvZMatEes(CkReductionMsg *)
Reduction client of zmat : Here everyone must have the same iteration ///////////////////////////////...
Definition: CP_State_RealParticlePlane.C:789
FFTcache::doNlFFTRtoG_Rchare
void doNlFFTRtoG_Rchare(complex *, double *, int, int, int, int)
non-local : Rchare : data(x,y,z) -> data(gx,gy,z) : backward ////////////////////////////////////////...
Definition: fftCache.C:696
CP_State_RealParticlePlane::planeRedSectDone
void planeRedSectDone(CkReductionMsg *m)
All cookies initialized for zmat plane reduction ////////////////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:342
UberCollection
holds the UberIndex and the offset for proxies
Definition: Uber.h:68
CP_State_RealParticlePlane::printEnlR
void printEnlR(CkReductionMsg *m)
Energy reduction client for ees method! /////////////////////////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:60
CP_State_RealParticlePlane::calcReductionPlaneNum
int calcReductionPlaneNum(int)
spread the reduction plane numbers around to minimize map collisions
Definition: CP_State_RealParticlePlane.C:1301
CP_State_RealParticlePlane::computeAtmForcEes
void computeAtmForcEes(CompAtmForcMsg *msg)
Use Zmat and ProjPsi to get atmForces, Energy and psiforces /////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:850
CP_State_RealParticlePlane::init
void init()
Post construction initialization.
Definition: CP_State_RealParticlePlane.C:194
instControllerProxy
CProxy_InstanceController instControllerProxy
Definition: cpaimd.ci:236
CP_State_RealParticlePlane::FFTNLEesBckR
void FFTNLEesBckR()
Do the FFT of projPsi(x,y,z) to get projPsi(gx,gy,z) ////////////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:1131
CP_State_RealParticlePlane::sendToEesGPP
void sendToEesGPP()
Send the PsiForce NL FFT back to GparticlePlane.
Definition: CP_State_RealParticlePlane.C:1226
CP_State_RealParticlePlane::initComplete
void initComplete()
Initialization and registration done for this chare /////////////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:355
CP_State_RealParticlePlane::printEnlRSimp
void printEnlRSimp(double, int, int)
Energy reduction client for ees method! /////////////////////////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:80
CP_State_RealParticlePlane::~CP_State_RealParticlePlane
~CP_State_RealParticlePlane()
The destructor : never called directly but I guess migration needs it ///////////////////////////////...
Definition: CP_State_RealParticlePlane.C:376
CP_State_RealParticlePlane::setPlaneRedCookie
void setPlaneRedCookie(EnlCookieMsg *)
Recv a dummy message, 1 integer, and set the cookie monster /////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:1323
CP_State_RealParticlePlane::FFTNLEesFwdR
void FFTNLEesFwdR()
Complete the Forward FFT of psi-projector sent from g-space chares //////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:523
mCastGrpId
CkGroupID mCastGrpId
Multicast manager group that handles many mcast/redns in the code. Grep for info. ...
Definition: cpaimd.C:216
fftCacheSlab.h
Add type declarations for simulationConstants class (readonly vars) and once class for each type of o...
CP_State_RealParticlePlane::recvFromEesGPP
void recvFromEesGPP(NLFFTMsg *)
Recv FFT data/psi-projector from CP_StateParticlePlane kicking things off for this Chare...
Definition: CP_State_RealParticlePlane.C:406
eesCache::queryCacheRPP
void queryCacheRPP(int, int, int)
= realParticlePlane Cache Management tool
Definition: eesCache.C:394
CP_State_RealParticlePlane::enlSectDone
void enlSectDone(CkReductionMsg *m)
All cookies initialized for enl section, only reduction root receives this //////////////////////////...
Definition: CP_State_RealParticlePlane.C:329
CP_State_RealParticlePlane::recvZMatEesSimp
void recvZMatEesSimp(int, double *, int, int, int)
A chare can be behind by 1 iteration only.
Definition: CP_State_RealParticlePlane.C:708
cpaimd.h
Some basic data structures and the array map classes are defined here.
CP_State_RealParticlePlane::launchFFTControl
void launchFFTControl(int)
Control launch of FFT based on Rho having more of its act together //////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:1113
CP_State_RealParticlePlane::registrationDone
void registrationDone()
= Make sure everyone is registered on the 1st time step
Definition: CP_State_RealParticlePlane.C:1384
CP_State_Plane.h
eesCache
Group Container class : Only allowed chare data classes have data.
Definition: eesCache.h:18
CP_State_RealParticlePlane::setEnlCookie
void setEnlCookie(EnlCookieMsg *)
Recv a dummy message, 1 integer, and set the cookie monster /////////////////////////////////////////...
Definition: CP_State_RealParticlePlane.C:1355
CP_State_RealParticlePlane::computeZmatEes
void computeZmatEes()
Compute the Zmat elements of this iteration : Spawn the section reduction ///////////////////////////...
Definition: CP_State_RealParticlePlane.C:631
eesCache::registerCacheRPP
void registerCacheRPP(int)
= realParticlePlane Cache Management tool
Definition: eesCache.C:118