OpenAtom  Version1.5a
fftCache.C
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////=
2 ///////////////////////////////////////////////////////////////////////////////c
3 ///////////////////////////////////////////////////////////////////////////////=
4 /** \file fftCache.C
5  * Add functions to allow application programmers to initialize these and
6  * the corresponding functions in Charm++ to call these functions with
7  * appropriate parameters
8  */
9 ///////////////////////////////////////////////////////////////////////////////=
10 
11 #include "utility/util.h"
12 #include "debug_flags.h"
13 #include <cmath>
14 #include "main/cpaimd.h"
15 #include "fft_slab_ctrl/fftCacheSlab.h"
16 #include "cp_state_ctrl/CP_State_Plane.h"
17 #include "src_mathlib/mathlib.h"
18 
19 //#define TEST_ALIGN
20 ///////////////////////////////////////////////////////////////////////////////=
21 
22 CPcharmParaInfo *sim = CPcharmParaInfo::get();
23 extern CkVec <CProxy_FFTcache> UfftCacheProxy;
24 extern Config config;
25 extern int nstates;
26 extern int sizeX;
27 CmiNodeLock FFTcache::fftw_plan_lock;
28 ///////////////////////////////////////////////////////////////////////////////=
29 
30 void initFFTLock(void) {
31  FFTcache::fftw_plan_lock = CmiCreateLock();
32 }
33 
34 ///////////////////////////////////////////////////////////////////////////////=
35 ///////////////////////////////////////////////////////////////////////////////c
36 ///////////////////////////////////////////////////////////////////////////////=
37 ///FFTcache - sets up a fftcache on each processor
38 ///////////////////////////////////////////////////////////////////////////////=
39 FFTcache::FFTcache(
40  int _ngrida, int _ngridb,int _ngridc,
41  int _ngridaEext, int _ngridbEext, int _ngridcEext,
42  int _ees_eext_on, int _ngridaNL, int _ngridbNL, int _ngridcNL,
43  int _ees_NL_on, int _nlines_max, int _nlines_max_rho,
44  int _nchareGState, int _nchareRState,
45  int _nchareGNL, int _nchareRNL,
46  int _nchareGRho, int _nchareRRho, int _nchareRRhoTot,
47  int _nchareGEext, int _nchareREext, int _nchareREextTot,
48  int *numGState, int *numRXState, int *numRYState, int *numRYStateLower,
49  int *numGNL, int *numRXNL, int *numRYNL, int *numRYNLLower,
50  int *numGRho, int *numRXRho, int *numRYRho ,
51  int *numGEext, int *numRXEext, int *numRYEext ,
52  int _fftopt, int _nsplitR, int _nsplitG,
53  int _rhoRsubPlanes, UberCollection _thisInstance): thisInstance(_thisInstance){
54 ///////////////////////////////////////////////////////////////////////////////=
55 /// Local Variables
56  if ( CmiMyRank() == 0 ) {
57  //fftw_plan_lock = CmiCreateLock();
58  }
59 
60  int size[3];
61  complex *cin; complex *cout;
62  double *din; double *dout;
63 
64 ///////////////////////////////////////////////////////////////////////////////=
65 /// Copy out
66 
67  cacheMemFlag = 0;
68 
69  ngrida = _ngrida;
70  ngridb = _ngridb;
71  ngridc = _ngridc;
72 
73  ngridaEext = _ngridaEext;
74  ngridbEext = _ngridbEext;
75  ngridcEext = _ngridcEext;
76  ees_eext_on = _ees_eext_on;
77 
78  ngridaNL = _ngridaNL;
79  ngridbNL = _ngridbNL;
80  ngridcNL = _ngridcNL;
81  ees_NL_on = _ees_NL_on;
82 
83  int nlines_max = _nlines_max;
84  int nlines_max_rho = _nlines_max_rho;
85 
86  int iopt = _fftopt;
87 
88  nchareGState = _nchareGState;
89  nchareRState = _nchareRState;
90  nchareGNL = _nchareGNL;
91  nchareRNL = _nchareRNL;
92  nchareGRho = _nchareGRho;
93  nchareRRho = _nchareRRho;
94  nchareRRhoTot = _nchareRRhoTot;
95  nchareGEext = _nchareGEext;
96  nchareREext = _nchareREext;
97  nchareREextTot = _nchareREextTot;
98  rhoRsubPlanes = _rhoRsubPlanes;
99 
100  nsplitR = _nsplitR;
101  nsplitG = _nsplitG;
102 
103 ///////////////////////////////////////////////////////////////////////////////=
104 /// Set the fft stuff generically
105 
106  int skipR = 1;
107  int skipC = 1;
108  int nf_max = 0;
109  int nwork1 = 0;
110  int nwork2 = 0;
111  int nwork2T = 0;
112  int plus = 1;
113  int mnus = -1;
114  int unit = 1;
115  double scale = 1.0;
116 
117  if(iopt<0 || iopt>1){
118  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
119  CkPrintf("Bad fft option, dude! %d\n",iopt);
120  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
121  CkExit();
122  }//endif
123  CmiLock(fftw_plan_lock);
124 ///////////////////////////////////////////////////////////////////////////////=
125 /// Density, State and EES Scratch
126 
127 
128  int pGsize = nlines_max*ngridc; // z-collections
129  int pGsizeHart = nlines_max_rho*ngridc; // z-collections
130  int pGsizeNL = nlines_max*ngridcNL; // z-collections
131  int pGsizeEext = nlines_max_rho*ngridcEext; // z-collections
132 
133  int pRsize = (ngrida/2+1)*ngridb; // x-y plane
134  int pRsizeNL = (ngridaNL/2+1)*ngridbNL; // x-y plane
135  if(!config.doublePack){
136  pRsize = 2*ngrida*ngridb; // x-y plane
137  pRsizeNL = 2*ngridaNL*ngridbNL; // x-y plane
138  }//endif
139 
140  int pmax = pGsize;
141  pmax = (pmax > pRsize ) ? pmax : pRsize;
142  pmax = (pmax > pGsizeHart) ? pmax : pGsizeHart;
143 
144  if(ees_eext_on==1){pmax = (pmax>pGsizeEext) ? pmax : pGsizeEext;}
145  if(ees_NL_on ==1){pmax = (pmax>pGsizeNL) ? pmax : pGsizeNL;}
146  if(ees_NL_on ==1){pmax = (pmax>pRsizeNL) ? pmax : pRsizeNL;}
147 
148  tmpData = (complex*) fftw_malloc(pmax*sizeof(complex));
149  tmpDataR = reinterpret_cast<double*> (tmpData);
150 
151 ///////////////////////////////////////////////////////////////////////////////=
152 /// State plans : funky names : non-cubic broken
153 
154  //------------------------------------------------
155  // Double Pack Plans
156 
157  nf_max = (ngrida > ngridb) ? ngrida : ngridb;
158  nf_max = (ngridc > nf_max) ? ngridc : nf_max;
159  nwork1 = 0; nwork2 = 0; nwork2T = 0;
160  if(iopt==1){
161  int iadd = (int) (2.3*(double)nf_max);
162  if(nf_max<=2048){iadd=0;}
163  nwork1 = 20000+iadd;
164  nwork2 = nwork1;
165  nwork2T = nwork1 + (2*nf_max+256)*64;
166  }//endif
167  if(config.doublePack){
168  skipR = ngrida+2;
169  skipC = ngrida/2+1;
170  }else{
171  skipR = 2*ngrida;
172  skipC = ngrida;
173  }//endif
174 
175  initFFTholder (&fwdYPlanState, &iopt,&nwork1,&nwork2T,&scale,&plus,&ngridb,
176  &skipC,&unit,nchareRState,&nsplitR,numRYState);
177  initFFTholder (&bwdYPlanState, &iopt,&nwork1,&nwork2T,&scale,&mnus,&ngridb,
178  &skipC,&unit,nchareRState,&nsplitR,numRYState);
179  if(!config.doublePack){
180  initFFTholder (&fwdYPlanStateLower, &iopt,&nwork1,&nwork2T,&scale,&plus,&ngridb,
181  &skipC,&unit,nchareRState,&nsplitR,numRYStateLower);
182  initFFTholder (&bwdYPlanStateLower, &iopt,&nwork1,&nwork2T,&scale,&mnus,&ngridb,
183  &skipC,&unit,nchareRState,&nsplitR,numRYStateLower);
184  }//endif
185 
186  if(config.doublePack){
187  initCRFFTholder(&fwdXPlanState, &iopt,&nwork1,&nwork2,&scale,&plus,&ngrida,
188  &skipR,&skipC,nchareRState,&nsplitR,numRXState);
189  initRCFFTholder(&bwdXPlanState, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngrida,
190  &skipR,&skipC,nchareRState,&nsplitR,numRXState);
191  }else{
192  initFFTholder (&fwdXPlanStateK, &iopt,&nwork1,&nwork2,&scale,&plus,&ngrida,
193  &unit,&skipC,nchareRState,&nsplitR,numRXState);
194  initFFTholder (&bwdXPlanStateK, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngrida,
195  &unit,&skipC,nchareRState,&nsplitR,numRXState);
196  }//endif
197 
198  initFFTholder (&fwdZPlanState, &iopt,&nwork1,&nwork2,&scale,&plus,&ngridc,
199  &unit, &ngridc,nchareGState,&nsplitG,numGState);
200  initFFTholder (&bwdZPlanState, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngridc,
201  &unit, &ngridc,nchareGState,&nsplitG,numGState);
202  char wisstring[1024];
203  snprintf(wisstring,1024,"fftwwisdom_%d.dat",CkNumPes());
204  if(iopt==0){
205  /**
206  * If there is a wisdom file already made for this decomp, load
207  * it. We assume here that the rho decomposition choice changes
208  * only with the number of processors. It should be noted that
209  * those engaging in serious rho parameter tuning would do well
210  * to delete their wisdom files. FFTW will make a new plan if
211  * the wisdom file doesn't address the size of a plan anyway so
212  * this probably isn't a big deal.
213  */
214 
215  FILE *wisdomFile=fopen(wisstring,"r");
216  if(wisdomFile ==NULL && CkMyPe()==0){
217  CkPrintf("@@@@@@@@@@@@@@@@@@@@_warning_@@@@@@@@@@@@@@@@@@@@\n");
218  CkPrintf("Can't open wisdom file %s plan creation will take a little longer\n",wisstring);
219  CkPrintf("@@@@@@@@@@@@@@@@@@@@_warning_@@@@@@@@@@@@@@@@@@@@\n");
220  }else{
221 
222  if(wisdomFile != NULL && fftw_import_wisdom_from_file(wisdomFile)==FFTW_SUCCESS){
223  // only print this once
224  if(CkMyPe()==0) CkPrintf("[%d] FFTCache loaded FFTW Wisdom %s\n",CkMyPe(),wisstring);
225  }else{
226  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
227  CkPrintf("Wisdom load failed on CkMyPe() %d !\n",CkMyPe());
228  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
229  }//endif
230  fclose(wisdomFile);
231  }//endif
232  size[0] = ngrida; size[1] = ngridb; size[2] = 1;
233  if(config.doublePack){
234  fwdXPlanState.rfftwPlan = rfftwnd_create_plan(1, (const int*)size, FFTW_COMPLEX_TO_REAL,
235  FFTW_MEASURE | FFTW_IN_PLACE|FFTW_USE_WISDOM);
236  bwdXPlanState.rfftwPlan = rfftwnd_create_plan(1, (const int*)size, FFTW_REAL_TO_COMPLEX,
237  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
238  }else{
239  fwdXPlanStateK.fftwPlan = fftw_create_plan(ngrida, FFTW_FORWARD,
240  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
241  bwdXPlanStateK.fftwPlan = fftw_create_plan(ngrida, FFTW_BACKWARD,
242  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
243  }//endif
244  fwdYPlanState.fftwPlan = fftw_create_plan(ngridb, FFTW_FORWARD,
245  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
246  bwdYPlanState.fftwPlan = fftw_create_plan(ngridb, FFTW_BACKWARD,
247  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
248  if(!config.doublePack){
249  fwdYPlanStateLower.fftwPlan = fftw_create_plan(ngridb, FFTW_FORWARD,
250  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
251  bwdYPlanStateLower.fftwPlan = fftw_create_plan(ngridb, FFTW_BACKWARD,
252  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
253  }//endif
254  fwdZPlanState.fftwPlan = fftw_create_plan(ngridc,FFTW_FORWARD,
255  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
256  bwdZPlanState.fftwPlan = fftw_create_plan(ngridc,FFTW_BACKWARD,
257  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
258  }//endif
259 
260 ///////////////////////////////////////////////////////////////////////////////=
261 /// Density plans : funky names : non-cubic broken
262 
263  //------------------------------------------------
264  // Double Pack Plans
265 
266  nf_max = (ngrida > ngridb) ? ngrida : ngridb;
267  nf_max = (ngridc > nf_max) ? ngridc : nf_max;
268  nwork1 = 0; nwork2 = 0; nwork2T = 0;
269  if(iopt==1){
270  int iadd = (int) (2.3*(double)nf_max);
271  if(nf_max<=2048){iadd=0;}
272  nwork1 = 20000+iadd;
273  nwork2 = nwork1;
274  nwork2T = nwork1 + (2*nf_max+256)*64;
275  }//endif
276 
277  if(config.doublePack){
278  skipR = ngrida+2;
279  skipC = ngrida/2+1;
280  }else{
281  skipR = 2*ngrida;
282  skipC = ngrida;
283  }//endif
284 
285  if(rhoRsubPlanes==1){
286  initFFTholder (&fwdYPlanRho, &iopt,&nwork1,&nwork2T,&scale,&plus,&ngridb,
287  &skipC,&unit,nchareRRhoTot,&nsplitR,numRYRho);
288  initFFTholder (&bwdYPlanRho, &iopt,&nwork1,&nwork2T,&scale,&mnus,&ngridb,
289  &skipC,&unit,nchareRRhoTot,&nsplitR,numRYRho);
290  }else{
291  initFFTholder (&fwdYPlanRhoS,&iopt,&nwork1,&nwork2,&scale,&plus,&ngridb,
292  &unit, &ngridb,nchareRRhoTot,&nsplitR,numRYRho);
293  initFFTholder (&bwdYPlanRhoS,&iopt,&nwork1,&nwork2,&scale,&mnus,&ngridb,
294  &unit, &ngridb,nchareRRhoTot,&nsplitR,numRYRho);
295  }//endif
296  initCRFFTholder(&fwdXPlanRho, &iopt,&nwork1,&nwork2,&scale,&plus,&ngrida,
297  &skipR,&skipC,nchareRRhoTot,&nsplitR,numRXRho);
298  initRCFFTholder(&bwdXPlanRho, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngrida,
299  &skipR,&skipC,nchareRRhoTot,&nsplitR,numRXRho);
300  initFFTholder (&fwdZPlanRho, &iopt,&nwork1,&nwork2,&scale,&plus,&ngridc,
301  &unit, &ngridc,nchareGRho,&nsplitG,numGRho);
302  initFFTholder (&bwdZPlanRho, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngridc,
303  &unit, &ngridc,nchareGRho,&nsplitG,numGRho);
304  // I'm special
305  initFFTholder (&fwdZPlanRhoHart,&iopt,&nwork1,&nwork2,&scale,&plus,&ngridc,
306  &unit, &ngridc,nchareGEext,&nsplitG,numGEext);
307 
308  if(iopt==0){
309  size[0] = ngrida; size[1] = ngridb; size[2] = 1;
310  fwdXPlanRho.rfftwPlan = rfftwnd_create_plan(1, (const int*)size, FFTW_COMPLEX_TO_REAL,
311  FFTW_MEASURE | FFTW_IN_PLACE|FFTW_USE_WISDOM);
312  bwdXPlanRho.rfftwPlan = rfftwnd_create_plan(1, (const int*)size, FFTW_REAL_TO_COMPLEX,
313  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
314  fwdYPlanRho.fftwPlan = fftw_create_plan(ngridb, FFTW_FORWARD,
315  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
316  bwdYPlanRho.fftwPlan = fftw_create_plan(ngridb, FFTW_BACKWARD,
317  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
318  fwdYPlanRhoS.fftwPlan = fftw_create_plan(ngridb, FFTW_FORWARD,
319  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
320  bwdYPlanRhoS.fftwPlan = fftw_create_plan(ngridb, FFTW_BACKWARD,
321  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
322  fwdZPlanRho.fftwPlan = fftw_create_plan(ngridc,FFTW_FORWARD,
323  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
324  bwdZPlanRho.fftwPlan = fftw_create_plan(ngridc,FFTW_BACKWARD,
325  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
326  fwdZPlanRhoHart.fftwPlan=fftw_create_plan(ngridc,FFTW_FORWARD,
327  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
328  }//endif
329 
330 ///////////////////////////////////////////////////////////////////////////////=
331 /// Euler spline plans : better names : non-cubic broken
332 
333  //-----------------------------------
334  // Double Pack Plan NL :
335 
336  nf_max = (ngridaNL> ngridbNL) ? ngridaNL : ngridbNL;
337  nf_max = (ngridcNL > nf_max ) ? ngridcNL : nf_max;
338  nwork1 = 0; nwork2 = 0; nwork2T = 0;
339  if(iopt==1){
340  int iadd = (int) (2.3*(double)nf_max);
341  if(nf_max<=2048){iadd=0;}
342  nwork1 = 20000+iadd;
343  nwork2 = nwork1;
344  nwork2T = nwork1 + (2*nf_max+256)*64;
345  }//endif
346  if(config.doublePack){
347  skipR = ngridaNL+2;
348  skipC = ngridaNL/2+1;
349  }else{
350  skipR = 2*ngridaNL;
351  skipC = ngridaNL;
352  }//endif
353 
354  if(ees_NL_on) {
355  initFFTholder (&fwdYPlanNL,&iopt,&nwork1,&nwork2T,&scale,&plus,&ngridbNL,
356  &skipC,&unit,nchareRNL,&nsplitR,numRYNL);
357  initFFTholder (&bwdYPlanNL,&iopt,&nwork1,&nwork2T,&scale,&mnus,&ngridbNL,
358  &skipC,&unit,nchareRNL,&nsplitR,numRYNL);
359  if(!config.doublePack){
360  initFFTholder (&fwdYPlanNLLower,&iopt,&nwork1,&nwork2T,&scale,&plus,&ngridbNL,
361  &skipC,&unit,nchareRNL,&nsplitR,numRYNLLower);
362  initFFTholder (&bwdYPlanNLLower, &iopt,&nwork1,&nwork2T,&scale,&mnus,&ngridbNL,
363  &skipC,&unit,nchareRNL,&nsplitR,numRYNLLower);
364  }//endif
365 
366  if(config.doublePack){
367  initCRFFTholder(&fwdXPlanNL,&iopt,&nwork1,&nwork2,&scale,&plus,&ngridaNL,
368  &skipR,&skipC,nchareRNL,&nsplitR,numRXNL);
369  initRCFFTholder(&bwdXPlanNL,&iopt,&nwork1,&nwork2,&scale,&mnus,&ngridaNL,
370  &skipR,&skipC,nchareRNL,&nsplitR,numRXNL);
371  }else{
372  initFFTholder (&fwdXPlanNLK, &iopt,&nwork1,&nwork2,&scale,&plus,&ngridaNL,
373  &unit,&skipC,nchareRNL,&nsplitR,numRXNL);
374  initFFTholder (&bwdXPlanNLK, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngridaNL,
375  &unit,&skipC,nchareRNL,&nsplitR,numRXNL);
376  }//endif
377  initFFTholder (&fwdZPlanNL,&iopt,&nwork1,&nwork2,&scale,&plus,&ngridcNL,
378  &unit, &ngridcNL,nchareGNL,&nsplitG,numGNL);
379  initFFTholder (&bwdZPlanNL,&iopt,&nwork1,&nwork2,&scale,&mnus,&ngridcNL,
380  &unit, &ngridcNL,nchareGNL,&nsplitG,numGNL);
381  }//endif : ees is on
382 
383  if(iopt==0 && ees_NL_on){
384  size[0] = ngridaNL; size[1] = ngridbNL; size[2] = 1;
385  if(config.doublePack){
386  fwdXPlanNL.rfftwPlan = rfftwnd_create_plan(1, (const int*)size, FFTW_COMPLEX_TO_REAL,
387  FFTW_MEASURE | FFTW_IN_PLACE|FFTW_USE_WISDOM);
388  bwdXPlanNL.rfftwPlan = rfftwnd_create_plan(1, (const int*)size,FFTW_REAL_TO_COMPLEX,
389  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
390  }else{
391  fwdXPlanNLK.fftwPlan = fftw_create_plan(ngridaNL, FFTW_FORWARD,
392  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
393  bwdXPlanNLK.fftwPlan = fftw_create_plan(ngridaNL, FFTW_BACKWARD,
394  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
395  }//endif
396  fwdYPlanNL.fftwPlan = fftw_create_plan(ngridbNL, FFTW_FORWARD,
397  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
398  bwdYPlanNL.fftwPlan = fftw_create_plan(ngridbNL, FFTW_BACKWARD,
399  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
400  if(!config.doublePack){
401  fwdYPlanNLLower.fftwPlan = fftw_create_plan(ngridbNL, FFTW_FORWARD,
402  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
403  bwdYPlanNLLower.fftwPlan = fftw_create_plan(ngridbNL, FFTW_BACKWARD,
404  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
405  }//endif
406  fwdZPlanNL.fftwPlan = fftw_create_plan(ngridcNL,FFTW_FORWARD,
407  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
408  bwdZPlanNL.fftwPlan = fftw_create_plan(ngridcNL,FFTW_BACKWARD,
409  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
410  }//end if
411 
412 ///////////////////////////////////////////////////////////////////////////////=
413 /// Euler spline plans : better names : non-cubic broken
414 
415  //-----------------------------------
416  // Double Pack Plan Eext :
417 
418  nf_max = (ngridaEext> ngridbEext) ? ngridaEext : ngridbEext;
419  nf_max = (ngridcEext > nf_max ) ? ngridcEext : nf_max;
420  nwork1 = 0; nwork2 = 0; nwork2T = 0;
421  if(iopt==1){
422  int iadd = (int) (2.3*(double)nf_max);
423  if(nf_max<=2048){iadd=0;}
424  nwork1 = 20000+iadd;
425  nwork2 = nwork1;
426  nwork2T = nwork1 + (2*nf_max+256)*64;
427  }//endif
428  skipR = ngridaEext+2;
429  skipC = ngridaEext/2+1;
430 
431  if(ees_eext_on){
432  if(rhoRsubPlanes==1){
433  initFFTholder (&fwdYPlanEext ,&iopt,&nwork1,&nwork2T,&scale,&plus,&ngridbEext,
434  &skipC,&unit,nchareREextTot,&nsplitR,numRYEext);
435  initFFTholder (&bwdYPlanEext ,&iopt,&nwork1,&nwork2T,&scale,&mnus,&ngridbEext,
436  &skipC,&unit,nchareREextTot,&nsplitR,numRYEext);
437  }else{
438  initFFTholder (&fwdYPlanEextS,&iopt,&nwork1,&nwork2,&scale,&plus,&ngridbEext,
439  &unit, &ngridbEext,nchareREextTot,&nsplitR,numRYEext);
440  initFFTholder (&bwdYPlanEextS,&iopt,&nwork1,&nwork2,&scale,&mnus,&ngridbEext,
441  &unit, &ngridbEext,nchareREextTot,&nsplitR,numRYEext);
442  }//endif
443  initCRFFTholder(&fwdXPlanEext, &iopt,&nwork1,&nwork2,&scale,&plus,&ngridaEext,
444  &skipR,&skipC,nchareREextTot,&nsplitR,numRXEext);
445  initRCFFTholder(&bwdXPlanEext, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngridaEext,
446  &skipR,&skipC,nchareREextTot,&nsplitR,numRXEext);
447  initFFTholder (&fwdZPlanEext, &iopt,&nwork1,&nwork2,&scale,&plus,&ngridcEext,
448  &unit, &ngridcEext,nchareGEext,&nsplitG,numGEext);
449  initFFTholder (&bwdZPlanEext, &iopt,&nwork1,&nwork2,&scale,&mnus,&ngridcEext,
450  &unit, &ngridcEext,nchareGEext,&nsplitG,numGEext);
451  }//endif
452  if(iopt==0){
453  size[0] = ngridaEext; size[1] = ngridbEext; size[2] = 1;
454  fwdXPlanEext.rfftwPlan = rfftwnd_create_plan(1, (const int*)size,FFTW_COMPLEX_TO_REAL,
455  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
456  bwdXPlanEext.rfftwPlan = rfftwnd_create_plan(1, (const int*)size,FFTW_REAL_TO_COMPLEX,
457  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
458  fwdYPlanEext.fftwPlan = fftw_create_plan(ngridbEext, FFTW_FORWARD,
459  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
460  bwdYPlanEext.fftwPlan = fftw_create_plan(ngridbEext, FFTW_BACKWARD,
461  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
462  fwdYPlanEextS.fftwPlan = fftw_create_plan(ngridbEext, FFTW_FORWARD,
463  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
464  bwdYPlanEextS.fftwPlan = fftw_create_plan(ngridbEext, FFTW_BACKWARD,
465  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
466  fwdZPlanEext.fftwPlan = fftw_create_plan(ngridcEext,FFTW_FORWARD,
467  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
468  bwdZPlanEext.fftwPlan = fftw_create_plan(ngridcEext,FFTW_BACKWARD,
469  FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
470  }//end switch
471  if(iopt==0 && CkMyPe()==0){
472  FILE *wisdomFile=fopen(wisstring,"w");
473  if(wisdomFile ==NULL){
474  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
475  CkPrintf("Can't open wisdom file %s for writing\n",wisstring);
476  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
477  }else{
478  fftw_export_wisdom_to_file(wisdomFile);
479  CkPrintf("Wisdom %s written\n",wisstring);
480  fclose(wisdomFile);
481  }//endif
482  }// end if opt==0
483  CmiUnlock(fftw_plan_lock);
484 //------------------------------------------------------------------------------
485  }//end routine
486 ///////////////////////////////////////////////////////////////////////////////=
487 
488 
489 
490 ///////////////////////////////////////////////////////////////////////////////
491 /// packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz
492 ///////////////////////////////////////////////////////////////////////////////=
493 ///////////////////////////////////////////////////////////////////////////////c
494 ///////////////////////////////////////////////////////////////////////////////=
495 
496 void FFTcache::expandGSpace(complex* data, complex *packedData,
497  RunDescriptor *runs, int numRuns, int numFull,
498  int numPoints, int nfftz)
499 
500 ///////////////////////////////////////////////////////////////////////////////=
501  {//begin routine
502 ///////////////////////////////////////////////////////////////////////////////=
503 /// Expand out lines of z so that an FFT can be performed on them
504 /// The ``run'' stores each line in two parts
505 /// -3,-2,-1 have offset, joff = r*nffz + nfftz
506 /// 0,1,2,3 have offset, joff = r*nfftz
507 /// runs[r].z stores kz if kz>0 and nfftz-kz if kz<0
508 /// The negative guys go 1st
509 /// Total size is nlines*nfftz where nlines=numRuns/2
510 
511 #define _BZERO_METH_
512 #ifdef _BZERO_METH_
513  bzero(data,sizeof(complex)*numFull);
514 #endif
515 
516  int nsub = (nfftz-runs[0].nz);
517  int koff = 0;
518  for (int r = 0,l=0; r < numRuns; r+=2,l++) {
519 
520  int joff1 = l*nfftz + runs[r].z + nsub; // k < 0
521  for (int i=0,j=joff1,k=koff; i<runs[r].length; i++,j++,k++) {
522  if(j<0 || j>nfftz*numRuns){printf("ecch!!! %d\n",j);CkExit();}
523  data[j] = packedData[k];
524  }//endfor
525  koff += runs[r].length;
526 
527  int r1=r+1;
528  int joff2 = l*nfftz + runs[r1].z; // k >= 0
529  for (int i=0,j=joff2,k=koff; i<runs[r1].length; i++,j++,k++) {
530  if(j<0 || j>nfftz*numRuns){printf("ecch!!! %d\n",j);CkExit();}
531  data[j] = packedData[k];
532  }//endfor
533  koff += runs[r1].length;
534 
535 #ifndef _BZERO_METH_
536  int joff3 = joff2+runs[r1].length;
537  for(int j=joff3;j<joff1;j++){data[j]=0.0;}
538 #endif
539  }//endfor
540 
541  CkAssert(numPoints == koff);
542 //------------------------------------------------------------------------------
543  }//end routine
544 ///////////////////////////////////////////////////////////////////////////////=
545 
546 
547 ///////////////////////////////////////////////////////////////////////////////
548 /// packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz
549 ///////////////////////////////////////////////////////////////////////////////=
550 ///////////////////////////////////////////////////////////////////////////////c
551 ///////////////////////////////////////////////////////////////////////////////=
552 
553 void FFTcache::packGSpace(complex* data, complex *packedData,
554  RunDescriptor *runs, int numRuns, int numFull,
555  int numPoints, int nfftz)
556 
557 ///////////////////////////////////////////////////////////////////////////////=
558  {//begin routine
559 ///////////////////////////////////////////////////////////////////////////////=
560 /// Contract lines of z after FFT has been performed
561 /// The ``run'' stores each line in two parts
562 /// -3,-2,-1 have offset, joff = r*nffz + nfftz
563 /// 0,1,2,3 have offset, joff = r*nfftz
564 /// runs[r].z stores kz if kz>0 and nfftz-kz if kz<0
565 /// The negative guys go 1st
566 /// Total size is nlines*nfftz where nlines=numRuns/2
567 
568  int koff = 0;
569  int nsub = (nfftz-runs[0].nz);
570  for (int r = 0,l=0; r < numRuns; r+=2,l++) {
571 
572  int joff = l*nfftz + runs[r].z + nsub;
573  for (int i=0,j=joff,k=koff; i<runs[r].length; i++,j++,k++) {
574  packedData[k] = data[j];
575  }//endfor
576  koff += runs[r].length;
577 
578  int r1=r+1;
579  joff = l*nfftz + runs[r1].z;
580  for (int i=0,j=joff,k=koff; i<runs[r1].length; i++,j++,k++) {
581  packedData[k] = data[j];
582  }//endfor
583  koff += runs[r1].length;
584  }//endfor
585 
586  CkAssert(numPoints == koff);
587 
588 //------------------------------------------------------------------------------
589  }//end routine
590 ///////////////////////////////////////////////////////////////////////////////=
591 
592 
593 ///////////////////////////////////////////////////////////////////////////////
594 /// non-local : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward
595 ///////////////////////////////////////////////////////////////////////////////=
596 ///////////////////////////////////////////////////////////////////////////////c
597 ///////////////////////////////////////////////////////////////////////////////=
598 void FFTcache::doNlFFTRtoG_Gchare(complex *data_in,complex *data_out,
599  int numFull, int numPoints,
600  int numLines, int numRuns, RunDescriptor *runs,
601  int nfftz, int index){
602 ///////////////////////////////////////////////////////////////////////////////=
603 /// FFT in expanded form
604 
605  fft_split(
606  &bwdZPlanNL, // Z-direction backward plan
607  numLines, // # of ffts : one for every line of z in the chare
608  (fftw_complex *)data_in, //input data
609  1, //stride
610  nfftz, //distance between z-data sets
611  NULL, 0, 0, // input is ouput
612  config.fftprogresssplit,// split parameter
613  index
614  );
615 
616 ///////////////////////////////////////////////////////////////////////////////=
617 /// Pack for computing
618 
619  packGSpace(data_in,data_out,runs,numRuns,numFull,numPoints,nfftz); // data_out is upacked
620 
621 //------------------------------------------------------------------------------
622  }//end routine
623 ///////////////////////////////////////////////////////////////////////////////=
624 
625 
626 
627 ///////////////////////////////////////////////////////////////////////////////
628 /// non-local : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward
629 ///////////////////////////////////////////////////////////////////////////////=
630 ///////////////////////////////////////////////////////////////////////////////c
631 ///////////////////////////////////////////////////////////////////////////////=
632 void FFTcache::doNlFFTGtoR_Gchare(complex *data_in,complex *data_out,
633  int numFull, int numPoints,int numLines, int numRuns,
634  RunDescriptor *runs, int nfftz, int index){
635 ///////////////////////////////////////////////////////////////////////////////=
636 /// Expand for ffting : Trickery
637 
638  expandGSpace(data_out,data_in,runs,numRuns,numFull,numPoints,nfftz);//data_out is expanded
639 
640 ///////////////////////////////////////////////////////////////////////////////=
641 /// FFT in expanded form
642 
643  fft_split(
644  &fwdZPlanNL, // Z-direction forward plan
645  numLines, // # of ffts : one for every line of z in the chare
646  (fftw_complex *)data_out, // data
647  1, //stride
648  nfftz, //distance between z-data sets
649  NULL, 0, 0, // input is ouput
650  config.fftprogresssplit, // split parameter
651  index
652  );
653 
654 //------------------------------------------------------------------------------
655  }//end routine
656 ///////////////////////////////////////////////////////////////////////////////=
657 
658 
659 ///////////////////////////////////////////////////////////////////////////////
660 /// Hartree : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward
661 ///////////////////////////////////////////////////////////////////////////////=
662 ///////////////////////////////////////////////////////////////////////////////c
663 ///////////////////////////////////////////////////////////////////////////////=
664 void FFTcache::doHartFFTGtoR_Gchare(complex *data_in,complex *data_out,
665  int numFull, int numPoints,int numLines, int numRuns,
666  RunDescriptor *runs, int nfftz, int index){
667 ///////////////////////////////////////////////////////////////////////////////=
668 /// Expand for ffting : Trickery
669 
670  expandGSpace(data_out,data_in,runs,numRuns,numFull,numPoints,nfftz);//data_out is expanded
671 
672 ///////////////////////////////////////////////////////////////////////////////=
673 /// FFT in expanded form
674 
675  fft_split(
676  &fwdZPlanRhoHart, // Z-direction forward plan
677  numLines, // # of ffts : one for every line of z in the chare
678  (fftw_complex *)data_out, // data
679  1, //stride
680  nfftz, //distance between z-data sets
681  NULL, 0, 0, // input is ouput
682  config.fftprogresssplit, // split parameter
683  index
684  );
685 
686 //------------------------------------------------------------------------------
687  }//end routine
688 ///////////////////////////////////////////////////////////////////////////////=
689 
690 
691 ///////////////////////////////////////////////////////////////////////////////
692 /// non-local : Rchare : data(x,y,z) -> data(gx,gy,z) : backward
693 ///////////////////////////////////////////////////////////////////////////////=
694 ///////////////////////////////////////////////////////////////////////////////c
695 ///////////////////////////////////////////////////////////////////////////////=
696 void FFTcache::doNlFFTRtoG_Rchare(complex *dataC,double *dataR,int nplane_x,
697  int sizeX,int sizeY, int index){
698 ///////////////////////////////////////////////////////////////////////////////=
699 /// FFT along x first
700 
701  int stride;
702  if (config.doublePack){
703  stride = sizeX/2+1;
704  rfftwnd_real_to_complex_split(
705  &bwdXPlanNL, // backward plan
706  sizeY, // these many 1D ffts
707  (fftw_real *)dataR, // data set
708  1, // stride
709  (sizeX+2), // spacing between data sets
710  NULL,0,0, // input array is output array
711  config.fftprogresssplitReal,//
712  index
713  );
714 
715  if(bwdXPlanNL.option==0){
716  for (int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
717  }//endif
718  }else{
719  stride = sizeX;
720  fft_split(
721  &bwdXPlanNLK, // x-plan
722  sizeY, // how many
723  (fftw_complex *)dataC, // input data
724  1, // stride (x is inner)
725  sizeX, // array separation
726  NULL,0,0, // output = input = dataR real
727  config.fftprogresssplitReal, // splitting parameter
728  index
729  );
730  }//endif
731 
732 ///////////////////////////////////////////////////////////////////////////////=
733 /// FFT along y
734 
735  fft_split(
736  &bwdYPlanNL, // backward plan
737  nplane_x, // these many 1D ffts
738  (fftw_complex *)dataC, // data set
739  stride, // stride
740  1, // spacing between data sets
741  NULL,0,0, // input is output
742  config.fftprogresssplitReal,// progress splitting for BG/L
743  index
744  );
745 
746  if(!config.doublePack){
747  int nplane_xm1 = nplane_x - 1;
748  int ioff = sizeX-nplane_x+1;
749  fft_split(
750  &bwdYPlanNLLower, // y-plan for negative Gx's
751  nplane_xm1, // how many < sizeX/2 + 1
752  (fftw_complex *)(&dataC[ioff]),//input data
753  sizeX, // stride betwen elements (x is inner)
754  1, // array separation (nffty elements)
755  NULL,0,0, // output data is input data
756  config.fftprogresssplitReal, // splitting parameter
757  index
758  );
759  }//endif
760 
761 //------------------------------------------------------------------------------
762  }//end routine
763 ///////////////////////////////////////////////////////////////////////////////=
764 
765 
766 ///////////////////////////////////////////////////////////////////////////////
767 /// non-local : Rchare : data(gx,gy,z) -> data(x,y,z) : Forward
768 ///////////////////////////////////////////////////////////////////////////////=
769 ///////////////////////////////////////////////////////////////////////////////c
770 ///////////////////////////////////////////////////////////////////////////////=
771 void FFTcache::doNlFFTGtoR_Rchare(complex *dataC,double *dataR,int nplane_x,
772  int sizeX,int sizeY, int index){
773 ///////////////////////////////////////////////////////////////////////////////=
774 /// FFT along Y direction : Y moves with stride sizex/2+1 through memory
775 /// : nplane_x is spherical cutoff <= sizeX/2+1
776 
777  int stride;
778  if (config.doublePack){
779  stride = sizeX/2+1;
780  }else{
781  stride = sizeX;
782  }//endif
783 
784  fft_split(
785  &fwdYPlanNL, // y-plan for NL Ees method
786  nplane_x, // how many < sizeX/2 + 1
787  (fftw_complex *)(dataC), //input data
788  stride, // stride betwen elements (x is inner)
789  1, // array separation (nffty elements)
790  NULL,0,0, // output data is input data
791  config.fftprogresssplitReal, // splitting parameter
792  index
793  );
794 
795  if(!config.doublePack){
796  int nplane_xm1 = nplane_x - 1;
797  int ioff = sizeX-nplane_x+1;
798  fft_split(
799  &fwdYPlanNLLower, // y-plan for negative Gx's
800  nplane_xm1, // how many < sizeX/2 + 1
801  (fftw_complex *)(&dataC[ioff]),//input data
802  stride, // stride betwen elements (x is inner)
803  1, // array separation (nffty elements)
804  NULL,0,0, // output data is input data
805  config.fftprogresssplitReal, // splitting parameter
806  index
807  );
808  }//endif
809 
810  // fftw only gives you one sign for real to complex : so do it yourself
811  if(config.doublePack){
812  if(fwdXPlanNL.option==0){// this plan is not intialized if doublePack==0
813  for(int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
814  }//endif
815  }//endif
816 
817 ///////////////////////////////////////////////////////////////////////////////=
818 /// FFT along X direction : X moves with stride 1 through memory
819 
820  if(config.doublePack){
821  rfftwnd_complex_to_real_split(
822  &fwdXPlanNL, // x-plan for NL Ees method
823  sizeY, // how many
824  (fftw_complex *)dataC, // input data
825  1, // stride (x is inner)
826  stride, // array separation
827  NULL,0,0, // output = input = dataR real
828  config.fftprogresssplitReal, // splitting parameter
829  index
830  );
831  }else{
832  fft_split(
833  &fwdXPlanNLK, // x-plan
834  sizeY, // how many
835  (fftw_complex *)dataC, // input data
836  1, // stride (x is inner)
837  stride, // array separation=sizeX
838  NULL,0,0, // output = input = dataR real
839  config.fftprogresssplitReal, // splitting parameter
840  index
841  );
842  }//endif
843 
844 //------------------------------------------------------------------------------
845  }//end routine
846 ///////////////////////////////////////////////////////////////////////////////=
847 
848 
849 
850 ///////////////////////////////////////////////////////////////////////////////
851 /// Eext : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward
852 ///////////////////////////////////////////////////////////////////////////////=
853 ///////////////////////////////////////////////////////////////////////////////c
854 ///////////////////////////////////////////////////////////////////////////////=
855 void FFTcache::doEextFFTRtoG_Gchare(complex *data,int numFull, int numPoints,
856  int numLines, int numRuns, RunDescriptor *runs,
857  int nfftz, int index){
858 ///////////////////////////////////////////////////////////////////////////////=
859 /// FFT in expanded form
860 
861  fft_split(
862  &bwdZPlanEext, // Z-direction backward plan
863  numLines, // # of ffts : one for every line of z in the chare
864  (fftw_complex *)data, // input data
865  1, // stride
866  nfftz, // distance between z-data sets
867  NULL, 0, 0, // input is ouput
868  config.fftprogresssplit, // split parameter
869  index
870  );
871 
872 ///////////////////////////////////////////////////////////////////////////////=
873 /// pack for computing
874 
875  getCacheMem("doEextFFTRtoG_Gchare");
876  packGSpace(data,tmpData,runs,numRuns,numFull,numPoints,nfftz); // tmpdata returns packed
877  CmiMemcpy(data,tmpData,sizeof(complex)*numPoints); // data is expanded; cp in tmpdata;
878  freeCacheMem("doEextFFTRtoG_Gchare");
879 
880 //------------------------------------------------------------------------------
881  }//end routine
882 ///////////////////////////////////////////////////////////////////////////////=
883 
884 
885 ///////////////////////////////////////////////////////////////////////////////
886 /// Eext : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward
887 ///////////////////////////////////////////////////////////////////////////////=
888 ///////////////////////////////////////////////////////////////////////////////c
889 ///////////////////////////////////////////////////////////////////////////////=
890 void FFTcache::doEextFFTGtoR_Gchare(complex *data_in, complex *data_out,
891  int numFull, int numPoints,int numLines,
892  int numRuns, RunDescriptor *runs,
893  int nfftz, int index){
894 ///////////////////////////////////////////////////////////////////////////////=
895 /// expand for ffting
896 
897  expandGSpace(data_out,data_in,runs,numRuns,numFull,numPoints,nfftz);// data_out expanded
898 
899 ///////////////////////////////////////////////////////////////////////////////=
900 /// FFT
901 
902  fft_split(
903  &fwdZPlanEext, // Z-direction forward plan
904  numLines, // # of ffts : one for every line of z in the chare
905  (fftw_complex *)data_out,//output data
906  1, //stride
907  nfftz, //distance between z-data sets
908  NULL, 0, 0, // input is ouput
909  config.fftprogresssplit, // split parameter
910  index
911  );
912 
913 //------------------------------------------------------------------------------
914  }//end routine
915 ///////////////////////////////////////////////////////////////////////////////=
916 
917 
918 ///////////////////////////////////////////////////////////////////////////////
919 /// Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward
920 ///////////////////////////////////////////////////////////////////////////////=
921 ///////////////////////////////////////////////////////////////////////////////c
922 ///////////////////////////////////////////////////////////////////////////////=
923 void FFTcache::doEextFFTRtoG_Rchare(complex *dataC,double *dataR,int nplane_x,
924  int sizeX,int sizeY, int index){
925 ///////////////////////////////////////////////////////////////////////////////=
926 /// FFT along x first
927 
928  rfftwnd_real_to_complex_split(
929  &bwdXPlanEext, // backward plan
930  sizeY, // these many 1D ffts
931  (fftw_real *)dataR, // data set
932  1, // stride
933  (sizeX+2), // spacing between data sets
934  NULL,0,0, // input array is output array
935  config.fftprogresssplitReal,//
936  index
937  );
938 
939  int stride = sizeX/2+1;
940 
941  if(bwdXPlanEext.option==0){
942  for (int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
943  }//endif
944 
945 ///////////////////////////////////////////////////////////////////////////////=
946 /// FFT along y
947 
948  fft_split(
949  &bwdYPlanEext, // backward plan
950  nplane_x, // these many 1D ffts
951  (fftw_complex *)dataC, // data set
952  stride, // stride
953  1, // spacing between data sets
954  NULL,0,0, // input is output
955  config.fftprogresssplitReal,// progress splitting for BG/L
956  index
957  );
958 
959 //------------------------------------------------------------------------------
960  }//end routine
961 ///////////////////////////////////////////////////////////////////////////////=
962 
963 
964 ///////////////////////////////////////////////////////////////////////////////
965 /// Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward
966 ///////////////////////////////////////////////////////////////////////////////=
967 ///////////////////////////////////////////////////////////////////////////////c
968 ///////////////////////////////////////////////////////////////////////////////=
969 void FFTcache::doEextFFTRxToGx_Rchare(complex *dataC,double *dataR,int nplane_x,
970  int sizeX,int sizeY, int index){
971 ///////////////////////////////////////////////////////////////////////////////=
972 /// FFT along x first : stride 1
973 
974  rfftwnd_real_to_complex_split(
975  &bwdXPlanEext, // backward plan
976  sizeY, // these many 1D ffts
977  (fftw_real *)dataR, // data set
978  1, // stride
979  (sizeX+2), // spacing between data sets
980  NULL,0,0, // input array is output array
981  config.fftprogresssplitReal,//
982  index
983  );
984 
985  int stride = sizeX/2+1;
986  if(bwdXPlanEext.option==0){
987  for (int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
988  }//endif
989 
990 //------------------------------------------------------------------------------
991  }//end routine
992 ///////////////////////////////////////////////////////////////////////////////=
993 
994 
995 
996 ///////////////////////////////////////////////////////////////////////////////
997 /// Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward
998 ///////////////////////////////////////////////////////////////////////////////=
999 ///////////////////////////////////////////////////////////////////////////////c
1000 ///////////////////////////////////////////////////////////////////////////////=
1001 void FFTcache::doEextFFTRyToGy_Rchare(complex *dataC,double *dataR,int nplane_x,
1002  int sizeX,int sizeY, int index){
1003 ///////////////////////////////////////////////////////////////////////////////=
1004 /// FFT along y : stride 1
1005 
1006  fft_split(
1007  &bwdYPlanEextS, // backward plan
1008  nplane_x, // these many 1D ffts
1009  (fftw_complex *)dataC, // data set
1010  1, // stride
1011  sizeY, // spacing between data sets
1012  NULL,0,0, // input is output
1013  config.fftprogresssplitReal,// progress splitting for BG/L
1014  index
1015  );
1016 
1017 //------------------------------------------------------------------------------
1018  }//end routine
1019 ///////////////////////////////////////////////////////////////////////////////=
1020 
1021 
1022 
1023 ///////////////////////////////////////////////////////////////////////////////
1024 /// Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward
1025 ///////////////////////////////////////////////////////////////////////////////=
1026 ///////////////////////////////////////////////////////////////////////////////c
1027 ///////////////////////////////////////////////////////////////////////////////=
1028 void FFTcache::doEextFFTGtoR_Rchare(complex *dataC,double *dataR,int nplane_x,
1029  int sizeX,int sizeY, int index){
1030 ///////////////////////////////////////////////////////////////////////////////=
1031 /// FFT along Y direction : Y moves with stride sizex/2+1 through memory
1032 /// : nplane_x is spherical cutoff <= sizeX/2+1
1033 
1034  int stride = sizeX/2+1;
1035  fft_split(
1036  &fwdYPlanEext, // y-plan for NL Ees method
1037  nplane_x, // how many < sizeX/2 + 1
1038  (fftw_complex *)(dataC), //input data
1039  stride, // stride betwen elements (x is inner)
1040  1, // array separation (nffty elements)
1041  NULL,0,0, // output data is input data
1042  config.fftprogresssplitReal, // splitting parameter
1043  index
1044  );
1045 
1046  // fftw only gives you one sign for real to complex : so do it yourself
1047  if(fwdXPlanEext.option==0){
1048  for(int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1049  }//endif
1050 
1051 ///////////////////////////////////////////////////////////////////////////////=
1052 /// FFT along X direction : X moves with stride 1 through memory
1053 
1054  rfftwnd_complex_to_real_split(
1055  &fwdXPlanEext, // x-plan for NL Ees method
1056  sizeY, // how many
1057  (fftw_complex *)dataC, // input data
1058  1, // stride (x is inner)
1059  stride, // array separation
1060  NULL,0,0, // output = input = dataR real
1061  config.fftprogresssplitReal, // splitting parameter
1062  index
1063  );
1064 
1065 //------------------------------------------------------------------------------
1066  }//end routine
1067 ///////////////////////////////////////////////////////////////////////////////=
1068 
1069 
1070 ///////////////////////////////////////////////////////////////////////////////
1071 /// Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward
1072 ///////////////////////////////////////////////////////////////////////////////=
1073 ///////////////////////////////////////////////////////////////////////////////c
1074 ///////////////////////////////////////////////////////////////////////////////=
1075 void FFTcache::doEextFFTGxToRx_Rchare(complex *dataC,double *dataR,int nplane_x,
1076  int sizeX,int sizeY, int index){
1077 ///////////////////////////////////////////////////////////////////////////////=
1078 /// FFT along X direction : X moves with stride 1 through memory
1079 
1080  int stride = sizeX/2+1;
1081  if(fwdXPlanEext.option==0){
1082  for(int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1083  }//endif
1084 
1085  rfftwnd_complex_to_real_split(
1086  &fwdXPlanEext, // x-plan for NL Ees method
1087  sizeY, // how many
1088  (fftw_complex *)dataC, // input data
1089  1, // stride (x is inner)
1090  stride, // array separation
1091  NULL,0,0, // output = input = dataR real
1092  config.fftprogresssplitReal, // splitting parameter
1093  index
1094  );
1095 
1096 //------------------------------------------------------------------------------
1097  }//end routine
1098 ///////////////////////////////////////////////////////////////////////////////=
1099 
1100 
1101 ///////////////////////////////////////////////////////////////////////////////
1102 /// Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward
1103 ///////////////////////////////////////////////////////////////////////////////=
1104 ///////////////////////////////////////////////////////////////////////////////c
1105 ///////////////////////////////////////////////////////////////////////////////=
1106 void FFTcache::doEextFFTGyToRy_Rchare(complex *dataC,double *dataR,int nplane_x,
1107  int sizeX,int sizeY, int index){
1108 ///////////////////////////////////////////////////////////////////////////////=
1109 /// FFT along Y direction : Y moves with stride 1 through memory
1110 /// : nplane_x is spherical cutoff <= sizeX/2+1
1111 
1112  fft_split(
1113  &fwdYPlanEextS, // y-plan for NL Ees method
1114  nplane_x, // how many < sizeX/2 + 1
1115  (fftw_complex *)(dataC), //input data
1116  1, // stride betwen elements (x is inner)
1117  sizeY, // array separation (nffty elements)
1118  NULL,0,0, // output data is input data
1119  config.fftprogresssplitReal, // splitting parameter
1120  index
1121  );
1122 
1123 //------------------------------------------------------------------------------
1124  }//end routine
1125 ///////////////////////////////////////////////////////////////////////////////=
1126 
1127 
1128 
1129 
1130 ///////////////////////////////////////////////////////////////////////////////
1131 /// StatePlane : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward
1132 ///////////////////////////////////////////////////////////////////////////////=
1133 ///////////////////////////////////////////////////////////////////////////////c
1134 ///////////////////////////////////////////////////////////////////////////////=
1135 void FFTcache::doStpFFTRtoG_Gchare(complex *data_in,complex *data_out,
1136  int numFull, int numPoints,
1137  int numLines, int numRuns, RunDescriptor *runs,
1138  int nfftz, int index){
1139 ///////////////////////////////////////////////////////////////////////////////=
1140 /// FFT in expanded form
1141 
1142  fft_split(
1143  &bwdZPlanState, // Z-direction backward plan
1144  numLines, // # of ffts : one for every line of z in the chare
1145  (fftw_complex *)data_in,//input data
1146  1, //stride
1147  nfftz, //distance between z-data sets
1148  NULL, 0, 0, // input is ouput
1149  config.fftprogresssplit,// split parameter
1150  index
1151  );
1152 
1153 ///////////////////////////////////////////////////////////////////////////////=
1154 /// Pack for computing
1155 
1156  packGSpace(data_in,data_out,runs,numRuns,numFull,numPoints,nfftz);
1157 
1158 //------------------------------------------------------------------------------
1159  }//end routine
1160 ///////////////////////////////////////////////////////////////////////////////=
1161 
1162 
1163 
1164 ///////////////////////////////////////////////////////////////////////////////
1165 /// StatePlane : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward
1166 ///////////////////////////////////////////////////////////////////////////////=
1167 ///////////////////////////////////////////////////////////////////////////////c
1168 ///////////////////////////////////////////////////////////////////////////////=
1169 void FFTcache::doStpFFTGtoR_Gchare(complex *data_in,complex *data_out,
1170  int numFull, int numPoints,
1171  int numLines, int numRuns, RunDescriptor *runs,
1172  int nfftz, int index){
1173 ///////////////////////////////////////////////////////////////////////////////=
1174 /// Expand for ffting
1175 
1176  // data_out is expanded
1177  // data_in is contracted
1178  expandGSpace(data_out,data_in,runs,numRuns,numFull,numPoints,nfftz);
1179 
1180 ///////////////////////////////////////////////////////////////////////////////=
1181 /// FFT in expanded form
1182 
1183  fft_split(
1184  &fwdZPlanState, // Z-direction forward plan
1185  numLines, // # of ffts : one for every line of z in the chare
1186  (fftw_complex *)data_out,//input data
1187  1, //stride
1188  nfftz, //distance between z-data sets
1189  NULL, 0, 0, // input is ouput
1190  config.fftprogresssplit, // split parameter
1191  index
1192  );
1193 
1194 //------------------------------------------------------------------------------
1195  }//end routine
1196 ///////////////////////////////////////////////////////////////////////////////=
1197 
1198 
1199 
1200 ///////////////////////////////////////////////////////////////////////////////=
1201 ///////////////////////////////////////////////////////////////////////////////c
1202 ///////////////////////////////////////////////////////////////////////////////=
1203 void FFTcache::doStpFFTGtoR_Rchare(complex *dataC,double *dataR,int nplane_x,
1204  int sizeX,int sizeY, int index){
1205 ///////////////////////////////////////////////////////////////////////////////=
1206 /// FFT along Y direction : Y moves with stride sizex/2+1 through memory
1207 /// : nplane_x is spherical cutoff <= sizeX/2+1
1208 
1209  int stride;
1210 
1211  if(config.doublePack){
1212  stride = sizeX/2+1;
1213  }else{
1214  stride = sizeX;
1215  }//endif
1216 
1217  fft_split(
1218  &fwdYPlanState, // y-plan
1219  nplane_x, // how many < sizeX/2 + 1
1220  (fftw_complex *)(dataC), //input data
1221  stride, // stride betwen elements (x is inner)
1222  1, // array separation (nffty elements)
1223  NULL,0,0, // output data is input data
1224  config.fftprogresssplitReal, // splitting parameter
1225  index
1226  );
1227 
1228  if(!config.doublePack){
1229  int nplane_xm1 = nplane_x - 1;
1230  int ioff = sizeX-nplane_x+1;
1231  fft_split(
1232  &fwdYPlanStateLower, // y-plan for negative Gx's
1233  nplane_xm1, // how many < sizeX/2 + 1
1234  (fftw_complex *)(&dataC[ioff]),//input data
1235  stride, // stride betwen elements (x is inner)
1236  1, // array separation (nffty elements)
1237  NULL,0,0, // output data is input data
1238  config.fftprogresssplitReal, // splitting parameter
1239  index
1240  );
1241  }//endif
1242 
1243  // fftw only gives you one sign for real to complex : so do it yourself
1244 
1245  if(config.doublePack){
1246  if(fwdXPlanState.option==0){ // this plan is unintialized if doublePack==0
1247  for(int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1248  }//endif
1249  }//endif
1250 
1251 ///////////////////////////////////////////////////////////////////////////////=
1252 /// FFT along X direction : X moves with stride 1 through memory
1253 
1254  if(config.doublePack){
1255  rfftwnd_complex_to_real_split(
1256  &fwdXPlanState, // x-plan
1257  sizeY, // how many
1258  (fftw_complex *)dataC, // input data
1259  1, // stride (x is inner)
1260  stride, // array separation
1261  NULL,0,0, // output = input = dataR real
1262  config.fftprogresssplitReal, // splitting parameter
1263  index
1264  );
1265  }else{
1266  fft_split(
1267  &fwdXPlanStateK, // x-plan
1268  sizeY, // how many
1269  (fftw_complex *)dataC, // input data
1270  1, // stride (x is inner)
1271  stride, // array separation = sizeX
1272  NULL,0,0, // output = input
1273  config.fftprogresssplitReal, // splitting parameter
1274  index
1275  );
1276  }//endif
1277 
1278 //------------------------------------------------------------------------------
1279  }//end routine
1280 ///////////////////////////////////////////////////////////////////////////////=
1281 
1282 
1283 ///////////////////////////////////////////////////////////////////////////////
1284 ///
1285 ///////////////////////////////////////////////////////////////////////////////=
1286 ///////////////////////////////////////////////////////////////////////////////c
1287 ///////////////////////////////////////////////////////////////////////////////=
1288 void FFTcache::doStpFFTRtoG_Rchare(complex *dataC,double *dataR,int nplane_x,
1289  int sizeX,int sizeY, int index){
1290 ///////////////////////////////////////////////////////////////////////////////=
1291 /// FFT along x first
1292 
1293  int stride;
1294 
1295  if(config.doublePack){
1296  stride = sizeX/2+1;
1297  rfftwnd_real_to_complex_split(
1298  &bwdXPlanState, // backward X plan
1299  sizeY, // these many 1D ffts
1300  (fftw_real *)dataR, // data set
1301  1, // stride
1302  (sizeX+2), // spacing between data sets
1303  NULL,0,0, // input array is output array
1304  config.fftprogresssplitReal,//
1305  index
1306  );
1307 
1308  if(bwdXPlanState.option==0){
1309  for (int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1310  }//endif
1311  }else{
1312  stride = sizeX;
1313  fft_split(
1314  &bwdXPlanStateK, // x-plan
1315  sizeY, // how many
1316  (fftw_complex *)dataC, // input data
1317  1, // stride (x is inner)
1318  sizeX, // array separation
1319  NULL,0,0, // output = input = dataR real
1320  config.fftprogresssplitReal, // splitting parameter
1321  index
1322  );
1323  }//endif
1324 
1325 ///////////////////////////////////////////////////////////////////////////////=
1326 /// FFT along y
1327 
1328  fft_split(
1329  &bwdYPlanState, // backward Y plan
1330  nplane_x, // these many 1D ffts
1331  (fftw_complex *)dataC, // data set
1332  stride, // stride
1333  1, // spacing between data sets
1334  NULL,0,0, // input is output
1335  config.fftprogresssplitReal,// progress splitting for BG/L
1336  index
1337  );
1338 
1339  if(!config.doublePack){
1340  int nplane_xm1 = nplane_x - 1;
1341  int ioff = sizeX - nplane_x + 1;
1342  fft_split(
1343  &bwdYPlanStateLower, // y-plan for negative Gx's
1344  nplane_xm1, // how many < sizeX/2 + 1
1345  (fftw_complex *)(&dataC[ioff]),//input data
1346  sizeX, // stride betwen elements (x is inner)
1347  1, // array separation (nffty elements)
1348  NULL,0,0, // output data is input data
1349  config.fftprogresssplitReal, // splitting parameter
1350  index
1351  );
1352  }//endif
1353 
1354 //------------------------------------------------------------------------------
1355  }//end routine
1356 ///////////////////////////////////////////////////////////////////////////////=
1357 
1358 
1359 
1360 ///////////////////////////////////////////////////////////////////////////////
1361 /// rho Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward
1362 ///////////////////////////////////////////////////////////////////////////////=
1363 ///////////////////////////////////////////////////////////////////////////////c
1364 ///////////////////////////////////////////////////////////////////////////////=
1365 void FFTcache::doRhoFFTGtoR_Gchare(complex *data_in,complex *data_out,
1366  int numFull, int numPoints,
1367  int numLines, int numRuns, RunDescriptor *runs,
1368  int nfftz,int iexpand, int index){
1369 ///////////////////////////////////////////////////////////////////////////////=
1370 /// Expand for ffting : data_out is expanded
1371 
1372  if(iexpand==1){
1373  expandGSpace(data_out,data_in,runs,numRuns,numFull,numPoints,nfftz);
1374  }//endif
1375 
1376 ///////////////////////////////////////////////////////////////////////////////=
1377 /// FFT in expanded form
1378 
1379  fft_split(
1380  &fwdZPlanRho, // Z-direction forward plan
1381  numLines, // # of ffts : one for every line of z in the chare
1382  (fftw_complex *)data_out,//input data
1383  1, //stride
1384  nfftz, //distance between z-data sets
1385  NULL, 0, 0, // input is ouput
1386  config.fftprogresssplit, // split parameter
1387  index
1388  );
1389 
1390 //------------------------------------------------------------------------------
1391  }//end routine
1392 ///////////////////////////////////////////////////////////////////////////////=
1393 
1394 
1395 
1396 ///////////////////////////////////////////////////////////////////////////////
1397 /// Rho Plane : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward
1398 ///////////////////////////////////////////////////////////////////////////////=
1399 ///////////////////////////////////////////////////////////////////////////////c
1400 ///////////////////////////////////////////////////////////////////////////////=
1401 void FFTcache::doRhoFFTRtoG_Gchare(complex *data_in,complex *data_out,
1402  int numFull, int numPoints,
1403  int numLines, int numRuns, RunDescriptor *runs,
1404  int nfftz, int ipack, int index){
1405 ///////////////////////////////////////////////////////////////////////////////=
1406 /// FFT in expanded form
1407 
1408  fft_split(
1409  &bwdZPlanRho, // Z-direction backward plan
1410  numLines, // # of ffts : one for every line of z in the chare
1411  (fftw_complex *)data_in,//input data
1412  1, //stride
1413  nfftz, //distance between z-data sets
1414  NULL, 0, 0, // input is ouput
1415  config.fftprogresssplit,// split parameter
1416  index
1417  );
1418 
1419 ///////////////////////////////////////////////////////////////////////////////=
1420 /// Pack for computing if necessary : data_in is expanded : dataout is contracted
1421 
1422  if(ipack==1){
1423  packGSpace(data_in,data_out,runs,numRuns,numFull,numPoints,nfftz);
1424  }//endif
1425 
1426 //------------------------------------------------------------------------------
1427  }//end routine
1428 ///////////////////////////////////////////////////////////////////////////////=
1429 
1430 
1431 
1432 
1433 
1434 ///////////////////////////////////////////////////////////////////////////////=
1435 ///////////////////////////////////////////////////////////////////////////////c
1436 ///////////////////////////////////////////////////////////////////////////////=
1437 void FFTcache::doRhoFFTRtoG_Rchare(complex *dataC,double *dataR,int nplane_x,
1438  int sizeX,int sizeY, int index)
1439 ///////////////////////////////////////////////////////////////////////////////=
1440  {//begin routine
1441 ///////////////////////////////////////////////////////////////////////////////=
1442 /// FFT along x first
1443 
1444  rfftwnd_real_to_complex_split(
1445  &bwdXPlanRho, // backward X plan
1446  sizeY, // these many 1D ffts
1447  (fftw_real *)dataR, // data set
1448  1, // stride
1449  (sizeX+2), // spacing between data sets
1450  NULL,0,0, // input array is output array
1451  config.fftprogresssplitReal,//
1452  index
1453  );
1454 
1455  int stride = sizeX/2+1;
1456  if(bwdXPlanRho.option==0){
1457  for (int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1458  }//endif
1459 
1460 ///////////////////////////////////////////////////////////////////////////////=
1461 /// FFT along y
1462 
1463  fft_split(
1464  &bwdYPlanRho, // backward Y plan
1465  nplane_x, // these many 1D ffts
1466  (fftw_complex *)dataC, // data set
1467  stride, // stride
1468  1, // spacing between data sets
1469  NULL,0,0, // input is output
1470  config.fftprogresssplitReal,// progress splitting for BG/L
1471  index
1472  );
1473 
1474 //------------------------------------------------------------------------------
1475  }//end routine
1476 ///////////////////////////////////////////////////////////////////////////////=
1477 
1478 
1479 ///////////////////////////////////////////////////////////////////////////////=
1480 ///////////////////////////////////////////////////////////////////////////////c
1481 ///////////////////////////////////////////////////////////////////////////////=
1482 void FFTcache::doRhoFFTRxToGx_Rchare(complex *dataC,double *dataR,int nplane_x,
1483  int sizeX,int sizeY, int index)
1484 ///////////////////////////////////////////////////////////////////////////////=
1485  {//begin routine
1486 ///////////////////////////////////////////////////////////////////////////////=
1487 /// FFT along x first
1488 
1489  rfftwnd_real_to_complex_split(
1490  &bwdXPlanRho, // backward X plan
1491  sizeY, // these many 1D ffts
1492  (fftw_real *)dataR, // data set
1493  1, // stride
1494  (sizeX+2), // spacing between data sets in doubles
1495  NULL,0,0, // input array is output array
1496  config.fftprogresssplitReal,//
1497  index
1498  );
1499 
1500  int stride = sizeX/2+1;
1501  if(bwdXPlanRho.option==0){
1502  for (int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1503  }//endif
1504 
1505 //------------------------------------------------------------------------------
1506  }//end routine
1507 ///////////////////////////////////////////////////////////////////////////////=
1508 
1509 
1510 
1511 ///////////////////////////////////////////////////////////////////////////////=
1512 ///////////////////////////////////////////////////////////////////////////////c
1513 ///////////////////////////////////////////////////////////////////////////////=
1514 void FFTcache::doRhoFFTRyToGy_Rchare(complex *dataC,double *dataR,int nplane_x,
1515  int sizeX,int sizeY, int index)
1516 ///////////////////////////////////////////////////////////////////////////////=
1517  {//begin routine
1518 ///////////////////////////////////////////////////////////////////////////////=
1519 /// FFT along y
1520 
1521  fft_split(
1522  &bwdYPlanRhoS, // backward Y plan
1523  nplane_x, // these many 1D ffts
1524  (fftw_complex *)dataC, // data set
1525  1, // stride
1526  sizeY, // spacing between data sets
1527  NULL,0,0, // input is output
1528  config.fftprogresssplitReal,// progress splitting for BG/L
1529  index
1530  );
1531 
1532 //------------------------------------------------------------------------------
1533  }//end routine
1534 ///////////////////////////////////////////////////////////////////////////////=
1535 
1536 
1537 
1538 ///////////////////////////////////////////////////////////////////////////////=
1539 ///////////////////////////////////////////////////////////////////////////////c
1540 ///////////////////////////////////////////////////////////////////////////////=
1541 void FFTcache::doRhoFFTGtoR_Rchare(complex *dataC,double *dataR,int nplane_x,
1542  int sizeX,int sizeY, int index){
1543 ///////////////////////////////////////////////////////////////////////////////=
1544 /// FFT along Y direction : Y moves with stride sizex/2+1 through memory
1545 /// : nplane_x is spherical cutoff <= sizeX/2+1
1546 
1547  int stride = sizeX/2+1;
1548  fft_split(
1549  &fwdYPlanRho, // y-plan
1550  nplane_x, // how many < sizeX/2 + 1
1551  (fftw_complex *)(dataC), //input data
1552  stride, // stride betwen elements (x is inner)
1553  1, // array separation (nffty elements)
1554  NULL,0,0, // output data is input data
1555  config.fftprogresssplitReal, // splitting parameter
1556  index
1557  );
1558 
1559  // fftw only gives you one sign for real to complex : so do it yourself
1560  if(fwdXPlanRho.option==0){
1561  for(int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1562  }//endif
1563 
1564 ///////////////////////////////////////////////////////////////////////////////=
1565 /// FFT along X direction : X moves with stride 1 through memory
1566 
1567  rfftwnd_complex_to_real_split(
1568  &fwdXPlanRho, // x-plan
1569  sizeY, // how many
1570  (fftw_complex *)dataC, // input data
1571  1, // stride (x is inner)
1572  stride, // array separation
1573  NULL,0,0, // output = input = dataR real
1574  config.fftprogresssplitReal, // splitting parameter
1575  index
1576  );
1577 
1578 //------------------------------------------------------------------------------
1579  }//end routine
1580 ///////////////////////////////////////////////////////////////////////////////=
1581 
1582 
1583 ///////////////////////////////////////////////////////////////////////////////=
1584 ///////////////////////////////////////////////////////////////////////////////c
1585 ///////////////////////////////////////////////////////////////////////////////=
1586 void FFTcache::doRhoFFTGyToRy_Rchare(complex *dataC,double *dataR,int nplane_x,
1587  int sizeX,int sizeY, int index){
1588 ///////////////////////////////////////////////////////////////////////////////=
1589 /// FFT along Y direction : Y moves with stride sizex/2+1 through memory
1590 /// : nplane_x is spherical cutoff <= sizeX/2+1
1591 
1592  fft_split(
1593  &fwdYPlanRhoS, // y-plan
1594  nplane_x, // how many < sizeX/2 + 1
1595  (fftw_complex *)(dataC), //input data
1596  1, // stride betwen elements (x is inner)
1597  sizeY, // array separation (nffty elements)
1598  NULL,0,0, // output data is input data
1599  config.fftprogresssplitReal, // splitting parameter
1600  index
1601  );
1602 
1603 //------------------------------------------------------------------------------
1604  }//end routine
1605 ///////////////////////////////////////////////////////////////////////////////=
1606 
1607 
1608 ///////////////////////////////////////////////////////////////////////////////=
1609 ///////////////////////////////////////////////////////////////////////////////c
1610 ///////////////////////////////////////////////////////////////////////////////=
1611 void FFTcache::doRhoFFTGxToRx_Rchare(complex *dataC,double *dataR,int nplane_x,
1612  int sizeX,int sizeY, int index){
1613 ///////////////////////////////////////////////////////////////////////////////=
1614 /// FFT along X direction : X moves with stride 1 through memory
1615 
1616  int stride = sizeX/2+1;
1617  if(fwdXPlanRho.option==0){
1618  for(int i=0;i<stride*sizeY;i++){dataC[i].im = -dataC[i].im;}
1619  }//endif
1620 
1621  rfftwnd_complex_to_real_split(
1622  &fwdXPlanRho, // x-plan
1623  sizeY, // how many
1624  (fftw_complex *)dataC, // input data
1625  1, // stride (x is inner)
1626  stride, // array separation
1627  NULL,0,0, // output = input = dataR real
1628  config.fftprogresssplitReal, // splitting parameter
1629  index
1630  );
1631 
1632 //------------------------------------------------------------------------------
1633  }//end routine
1634 ///////////////////////////////////////////////////////////////////////////////=
1635 
1636 
1637 //////////////////////////////////////////////////////////////////////////////
1638 //////////////////////////////////////////////////////////////////////////////
1639 //////////////////////////////////////////////////////////////////////////////
1640 /**
1641  * split up an fft call into multiple invocations with cmiprogress
1642  * calls between them.
1643  */
1644 //////////////////////////////////////////////////////////////////////////////
1645 void fft_split(FFTplanHolder *fftplanholder, int howmany,
1646  fftw_complex *in, int istride_in, int idist_in,
1647  fftw_complex *out_in, int ostride_in, int odist_in, int split,
1648  int index)
1649 //////////////////////////////////////////////////////////////////////////////
1650  {// begin routine
1651 //////////////////////////////////////////////////////////////////////////////
1652 
1653  fftw_plan plan = fftplanholder->fftwPlan;
1654  int iopt = fftplanholder->option;
1655  int isign = fftplanholder->isign;
1656  int nfft = fftplanholder->nfft;
1657  int nwork1 = fftplanholder->nwork1;
1658  int nwork2 = fftplanholder->nwork2;
1659  double scale = fftplanholder->scale;
1660  double *work1r = NULL;
1661  double *work2r = NULL;
1662  double *work1s = fftplanholder->work1;
1663  double *work2s = fftplanholder->work2;
1664 
1665  int zero = 0;
1666  int istride = istride_in;
1667  int idist = idist_in;
1668  int ostride = ostride_in;
1669  int odist = odist_in;
1670  fftw_complex *out = out_in;
1671  if(iopt==1){
1672  ostride = fftplanholder->ostride;
1673  odist = fftplanholder->odist;
1674  out = (out==NULL) ? in : out ;
1675  int kkk = fftplanholder->mapp[index];
1676  work1r = fftplanholder->essl_work[kkk].work1;
1677  work2r = fftplanholder->essl_work[kkk].work2;
1678  }//endif
1679 
1680  int thismany = split;
1681  int inleft = howmany;
1682  int numsplits = howmany/split;
1683  if(numsplits<=0){numsplits=1;}
1684  if(howmany>split && howmany%split!=0){numsplits++;}
1685 
1686 //////////////////////////////////////////////////////////////////////////////
1687 
1688  int inoff=0;
1689  int outoff=0;
1690 
1691  for(int i=0;i<numsplits;i++){
1692  thismany = split;
1693  double *work1 = work1s; double *work2 = work2s;
1694  if(inleft<split){thismany=inleft;work1=work1r,work2=work2r;}
1695 #ifdef TEST_ALIGN
1696  CkAssert((unsigned int) &(in[inoff]) % 16 ==0);
1697  CkAssert((unsigned int) &(out[outoff]) % 16 ==0);
1698 #endif
1699  fftw_complex *dummy = (out==NULL) ? NULL : &(out[outoff]) ;
1700  switch(iopt){
1701  case 0:
1702  fftw(
1703  plan, // da plan
1704  thismany, // how many
1705  &(in[inoff]), // input data
1706  istride, // stride betwen elements
1707  idist, // array separation
1708  dummy, // output data (null if inplace)
1709  ostride, // stride betwen elements (0 inplace)
1710  odist); // array separation (0 inplace)
1711  break;
1712  case 1: dcftWrap(&zero,(complex *)&(in[inoff]),&istride,&idist,
1713  (complex *)dummy,&ostride,&odist,
1714  &nfft,&thismany,&isign,&scale,work1,&nwork1,work2, &nwork2);
1715  break;
1716  default : CkAbort("impossible fft iopt"); break;
1717  }//end switch
1718  inoff += idist*(thismany);
1719  outoff += odist*(thismany);
1720  inleft -= thismany;
1721  }//endfor
1722 
1723 //------------------------------------------------------------------------------
1724  }//end routine
1725 //////////////////////////////////////////////////////////////////////////////
1726 
1727 
1728 //////////////////////////////////////////////////////////////////////////////
1729 //////////////////////////////////////////////////////////////////////////////
1730 //////////////////////////////////////////////////////////////////////////////
1731 /**
1732  * split up an fft call into multiple invocations with cmiprogress
1733  * calls between them.
1734  */
1735 //////////////////////////////////////////////////////////////////////////////
1736 void rfftwnd_complex_to_real_split(RFFTplanHolder *rfftplanholder, int howmany,
1737  fftw_complex *in, int istride_in, int idist_in,
1738  fftw_real *out_in, int ostride_in, int odist_in, int split, int index)
1739 //////////////////////////////////////////////////////////////////////////////
1740  {//begin routine
1741 //////////////////////////////////////////////////////////////////////////////
1742 
1743  rfftwnd_plan plan = rfftplanholder->rfftwPlan;
1744  int iopt = rfftplanholder->option;
1745  int isign = rfftplanholder->isign;
1746  int nfft = rfftplanholder->nfft;
1747  int nwork1 = rfftplanholder->nwork1;
1748  int nwork2 = rfftplanholder->nwork2;
1749  double scale = rfftplanholder->scale;
1750  double *work1r = NULL;
1751  double *work2r = NULL;
1752  double *work1s = rfftplanholder->work1;
1753  double *work2s = rfftplanholder->work2;
1754 
1755  int zero = 0;
1756  int istride = istride_in;
1757  int idist = idist_in;
1758  int ostride = ostride_in;
1759  int odist = odist_in;
1760  fftw_real *out = out_in;
1761  if(iopt==1){
1762  ostride = rfftplanholder->ostride;
1763  odist = rfftplanholder->odist;
1764  if(out==NULL){
1765  out = reinterpret_cast<double*> (in);
1766  }//endif
1767  int kkk = rfftplanholder->mapp[index];
1768  work1r = rfftplanholder->essl_work[kkk].work1;
1769  work2r = rfftplanholder->essl_work[kkk].work2;
1770  }//endif
1771 
1772  int thismany = split;
1773  int inleft = howmany;
1774  int numsplits = howmany/split;
1775  if(numsplits<=0){numsplits=1;}
1776  if(howmany>split && howmany%split!=0){numsplits++;}
1777 
1778 //////////////////////////////////////////////////////////////////////////////
1779 
1780  int inoff=0;
1781  int outoff=0;
1782 
1783  for(int i=0;i<numsplits;i++){
1784  thismany=split;
1785  double *work1 = work1s; double *work2 = work2s;
1786  if(inleft<split){thismany=inleft;work1=work1r,work2=work2r;}
1787 #ifdef TEST_ALIGN
1788  CkAssert((unsigned int) &(in[inoff]) % 16 ==0);
1789  CkAssert((unsigned int) &(out[outoff]) % 16 ==0);
1790 #endif
1791  fftw_real *dummy = (out==NULL) ? NULL : &(out[outoff]);
1792  switch(iopt){
1793  case 0:
1794  rfftwnd_complex_to_real(
1795  plan, // da plan
1796  thismany, // how many
1797  &(in[inoff]), // input data
1798  istride, // stride betwen elements
1799  idist, // array separation
1800  dummy, // output data (null if inplace)
1801  ostride, // stride betwen elements (0 inplace)
1802  odist); // array separation (0 inplace)
1803  break;
1804  case 1: dcrftWrap(&zero,(complex *)&(in[inoff]),&idist,(double *)dummy,&odist,
1805  &nfft,&thismany,&isign,&scale,work1,&nwork1,work2,&nwork2); break;
1806  default : CkAbort("impossible fft iopt"); break;
1807  }//end switch
1808  inoff += idist*(thismany);
1809  outoff += odist*(thismany);
1810  inleft -= thismany;
1811 
1812  }//endfor
1813 
1814 //------------------------------------------------------------------------------
1815  }// end routine
1816 //////////////////////////////////////////////////////////////////////////////
1817 
1818 
1819 //////////////////////////////////////////////////////////////////////////////
1820 //////////////////////////////////////////////////////////////////////////////
1821 //////////////////////////////////////////////////////////////////////////////
1822 /**
1823  * split up an fft call into multiple invocations with cmiprogress
1824  * calls between them.
1825  */
1826 //////////////////////////////////////////////////////////////////////////////
1827 void rfftwnd_real_to_complex_split(RFFTplanHolder *rfftplanholder, int howmany,
1828  fftw_real *in, int istride_in, int idist_in,
1829  fftw_complex *out_in, int ostride_in, int odist_in, int split, int index)
1830 //////////////////////////////////////////////////////////////////////////////
1831  {// begin routine
1832 //////////////////////////////////////////////////////////////////////////////
1833 
1834  rfftwnd_plan plan = rfftplanholder->rfftwPlan;
1835  int iopt = rfftplanholder->option;
1836  int isign = rfftplanholder->isign;
1837  int nfft = rfftplanholder->nfft;
1838  int nwork1 = rfftplanholder->nwork1;
1839  int nwork2 = rfftplanholder->nwork2;
1840  double scale = rfftplanholder->scale;
1841  double *work1r = NULL;
1842  double *work2r = NULL;
1843  double *work1s = rfftplanholder->work1;
1844  double *work2s = rfftplanholder->work2;
1845 
1846  int zero = 0;
1847  int istride = istride_in;
1848  int idist = idist_in;
1849  int ostride = ostride_in;
1850  int odist = odist_in;
1851  fftw_complex *out = out_in;
1852  if(iopt==1){
1853  ostride = rfftplanholder->ostride;
1854  odist = rfftplanholder->odist;
1855  if(out==NULL){
1856  out = reinterpret_cast<fftw_complex*> (in);
1857  }//endif
1858  int kkk = rfftplanholder->mapp[index];
1859  work1r = rfftplanholder->essl_work[kkk].work1;
1860  work2r = rfftplanholder->essl_work[kkk].work2;
1861  }//endif
1862 
1863  int thismany = split;
1864  int inleft = howmany;
1865  int numsplits = howmany/split;
1866  if(numsplits<=0){numsplits=1;}
1867  if(howmany>split && howmany%split!=0){numsplits++;}
1868 
1869 //////////////////////////////////////////////////////////////////////////////
1870 
1871  int inoff = 0;
1872  int outoff = 0;
1873  for(int i=0;i<numsplits;i++){
1874  thismany=split;
1875  double *work1 = work1s; double *work2 = work2s;
1876  if(inleft<split){thismany=inleft;work1=work1r,work2=work2r;}
1877 #ifdef TEST_ALIGN
1878  CkAssert((unsigned int) &(in[inoff]) % 16 ==0);
1879  CkAssert((unsigned int) &(out[outoff]) % 16 ==0);
1880 #endif
1881  fftw_complex *dummy = (out==NULL) ? NULL : &(out[outoff]);
1882  switch(iopt){
1883  case 0:
1884  rfftwnd_real_to_complex(
1885  plan, // da plan
1886  thismany, // how many
1887  &(in[inoff]), // input data
1888  istride, // stride betwen elements
1889  idist, // array separation
1890  dummy, // output data (null if inplace)
1891  ostride, // stride betwen elements (0 inplace)
1892  odist); // array separation (0 inplace)
1893  break;
1894  case 1: drcftWrap(&zero,(double *)&(in[inoff]),&idist,(complex *)dummy,&odist,
1895  &nfft,&thismany,&isign,&scale,work1,&nwork1,work2,&nwork2); break;
1896  default : CkAbort("impossible fft iopt"); break;
1897  }//end switch
1898  inoff += idist*(thismany);
1899  outoff += odist*(thismany);
1900  inleft -= thismany;
1901 
1902  }// endfor
1903 
1904 //----------------------------------------------------------------------------
1905  }//end routine
1906 //////////////////////////////////////////////////////////////////////////////
1907 
1908 
1909 
1910 //////////////////////////////////////////////////////////////////////////////
1911 //////////////////////////////////////////////////////////////////////////////
1912 //////////////////////////////////////////////////////////////////////////////
1913 void initFFTholder(FFTplanHolder *plan, int *iopt,int *nwork1,int *nwork2, double *scale,
1914  int *isign, int *nfft, int *stride, int *skip, int nchare,int *nsplit,
1915  int *num){
1916 //////////////////////////////////////////////////////////////////////////////
1917 
1918  plan->option = iopt[0];
1919  plan->nwork1 = nwork1[0];
1920  plan->nwork2 = nwork2[0];
1921  plan->scale = scale[0];
1922  plan->isign = isign[0];
1923  plan->nfft = nfft[0];
1924  plan->ostride = stride[0];
1925  plan->odist = skip[0];
1926  plan->work1 = NULL;
1927  plan->work2 = NULL;
1928  plan->nsplit = nsplit[0];
1929  plan->nchare = nchare;
1930 
1931  if(iopt[0] == 1){
1932  complex x[2];
1933  int nval,nmax;
1934  int unit = 1;
1935  int *index = new int[nchare];
1936  int *mappI = new int[nchare];
1937  int *mapp = new int[nchare];
1938  make_essl_work_map(nchare,num,index,mappI,mapp,&nval,&nmax,nsplit[0]);
1939  plan->essl_work = new ESSL_WORK[nval];
1940  plan->nval = nval;
1941  plan->mapp = mapp;
1942  for(int i = 0; i<nval;i++){
1943  plan->essl_work[i].num = index[i];
1944  plan->essl_work[i].work1 = NULL;
1945  plan->essl_work[i].work2 = NULL;
1946  if(index[i]>0){
1947  double *work1 = (double*) fftw_malloc(nwork1[0]*sizeof(double));
1948  double *work2 = (double*) fftw_malloc(nwork2[0]*sizeof(double));
1949  dcftWrap(&unit,x,stride,skip,x,stride,skip,nfft,&index[i],isign,scale,
1950  work1,nwork1,work2,nwork2);
1951  plan->essl_work[i].work1 = work1;
1952  plan->essl_work[i].work2 = work2;
1953  }//endif
1954  }//endfor
1955  if(nsplit[0]<=nmax){
1956  double *work1 = (double*) fftw_malloc(nwork1[0]*sizeof(double));
1957  double *work2 = (double*) fftw_malloc(nwork2[0]*sizeof(double));
1958  dcftWrap(&unit,x,stride,skip,x,stride,skip,nfft,nsplit,isign,scale,
1959  work1,nwork1,work2,nwork2);
1960  plan->work1 = work1;
1961  plan->work2 = work2;
1962  }//endif
1963  delete [] index;
1964  delete [] mappI;
1965  }//endif
1966 
1967 //----------------------------------------------------------------------------
1968  }//end routine
1969 //////////////////////////////////////////////////////////////////////////////
1970 
1971 
1972 //////////////////////////////////////////////////////////////////////////////
1973 //////////////////////////////////////////////////////////////////////////////
1974 //////////////////////////////////////////////////////////////////////////////
1975 void initRCFFTholder(RFFTplanHolder *plan, int *iopt,int *nwork1,int *nwork2, double *scale,
1976  int *isign, int *nfft, int *skipR, int *skipC, int nchare,int *nsplit,
1977  int *num){
1978 //////////////////////////////////////////////////////////////////////////////
1979 
1980  plan->option = iopt[0];
1981  plan->nwork1 = nwork1[0];
1982  plan->nwork2 = nwork2[0];
1983  plan->scale = scale[0];
1984  plan->isign = isign[0];
1985  plan->nfft = nfft[0];
1986  plan->ostride = 1;
1987  plan->odist = skipC[0];
1988  plan->work1 = NULL;
1989  plan->work2 = NULL;
1990  plan->nsplit = nsplit[0];
1991  plan->nchare = nchare;
1992 
1993  if(iopt[0] == 1){
1994  complex xc[2]; double xr[2];
1995  int nval,nmax;
1996  int unit = 1;
1997  int *index = new int[nchare];
1998  int *mappI = new int[nchare];
1999  int *mapp = new int[nchare];
2000  make_essl_work_map(nchare,num,index,mappI,mapp,&nval,&nmax,nsplit[0]);
2001  plan->essl_work = new ESSL_WORK[nval];
2002  plan->nval = nval;
2003  plan->mapp = mapp;
2004  for(int i = 0; i<nval;i++){
2005  plan->essl_work[i].num = index[i];
2006  plan->essl_work[i].work1 = NULL;
2007  plan->essl_work[i].work2 = NULL;
2008  if(index[i]>0){
2009  double *work1 = (double*) fftw_malloc(nwork1[0]*sizeof(double));
2010  double *work2 = (double*) fftw_malloc(nwork2[0]*sizeof(double));
2011  drcftWrap(&unit,xr,skipR,xc,skipC,nfft,&index[i],isign,scale,work1,nwork1,work2,nwork2);
2012  plan->essl_work[i].work1 = work1;
2013  plan->essl_work[i].work2 = work2;
2014  }//endif
2015  }//endfor
2016  if(nsplit[0]<=nmax){
2017  double *work1 = (double*) fftw_malloc(nwork1[0]*sizeof(double));
2018  double *work2 = (double*) fftw_malloc(nwork2[0]*sizeof(double));
2019  drcftWrap(&unit,xr,skipR,xc,skipC,nfft,nsplit,isign,scale,work1,nwork1,work2,nwork2);
2020  plan->work1 = work1;
2021  plan->work2 = work2;
2022  }//endif
2023  delete [] index;
2024  delete [] mappI;
2025  }//endif
2026 
2027 //----------------------------------------------------------------------------
2028  }//end routine
2029 //////////////////////////////////////////////////////////////////////////////
2030 
2031 
2032 //////////////////////////////////////////////////////////////////////////////
2033 //////////////////////////////////////////////////////////////////////////////
2034 //////////////////////////////////////////////////////////////////////////////
2035 void initCRFFTholder(RFFTplanHolder *plan, int *iopt,int *nwork1,int *nwork2, double *scale,
2036  int *isign, int *nfft, int *skipR, int *skipC,int nchare, int *nsplit,
2037  int *num){
2038 //////////////////////////////////////////////////////////////////////////////
2039 
2040  plan->option = iopt[0];
2041  plan->nwork1 = nwork1[0];
2042  plan->nwork2 = nwork2[0];
2043  plan->scale = scale[0];
2044  plan->isign = isign[0];
2045  plan->nfft = nfft[0];
2046  plan->ostride = 1;
2047  plan->odist = skipR[0];
2048  plan->work1 = NULL;
2049  plan->work2 = NULL;
2050  plan->nsplit = nsplit[0];
2051  plan->nchare = nchare;
2052 
2053  if(iopt[0] == 1){
2054  complex xc[2]; double xr[2];
2055  int nval,nmax;
2056  int unit = 1;
2057  int *index = new int[nchare];
2058  int *mappI = new int[nchare];
2059  int *mapp = new int[nchare];
2060  make_essl_work_map(nchare,num,index,mappI,mapp,&nval,&nmax,nsplit[0]);
2061  plan->essl_work = new ESSL_WORK[nval];
2062  plan->nval = nval;
2063  plan->mapp = mapp;
2064  for(int i = 0; i<nval;i++){
2065  plan->essl_work[i].num = index[i];
2066  plan->essl_work[i].work1 = NULL;
2067  plan->essl_work[i].work2 = NULL;
2068  if(index[i]>0){
2069  double *work1 = (double*) fftw_malloc(nwork1[0]*sizeof(double));
2070  double *work2 = (double*) fftw_malloc(nwork2[0]*sizeof(double));
2071  dcrftWrap(&unit,xc,skipC,xr,skipR,nfft,&index[i],isign,scale,work1,nwork1,work2,nwork2);
2072  plan->essl_work[i].work1 = work1;
2073  plan->essl_work[i].work2 = work2;
2074  }//endif
2075  }//endfor
2076  if(nsplit[0]<=nmax){
2077  double *work1 = (double*) fftw_malloc(nwork1[0]*sizeof(double));
2078  double *work2 = (double*) fftw_malloc(nwork2[0]*sizeof(double));
2079  dcrftWrap(&unit,xc,skipC,xr,skipR,nfft,nsplit,isign,scale,work1,nwork1,work2,nwork2);
2080  plan->work1 = work1;
2081  plan->work2 = work2;
2082  }//endif
2083  delete [] index;
2084  delete [] mappI;
2085  }//endif
2086 
2087 //----------------------------------------------------------------------------
2088  }//end routine
2089 //////////////////////////////////////////////////////////////////////////////
2090 
2091 ///////////////////////////////////////////////////////////////////////////=
2092 ///////////////////////////////////////////////////////////////////////////c
2093 ///////////////////////////////////////////////////////////////////////////=
2094 void make_essl_work_map(int nchare,int *num,int *index,int *mappI,int *mapp,
2095  int *nval_out, int *nmax_out,int nsplit){
2096 ///////////////////////////////////////////////////////////////////////////=
2097 
2098  int nmax = num[0];
2099  for(int i=0;i<nchare;i++){
2100  nmax = (nmax > num[i] ? nmax : num[i]);
2101  mappI[i] = i;
2102  index[i] = (num[i] % nsplit);
2103  }//endfor
2104 
2105  sort_commence(nchare,index,mappI);
2106 
2107  int j = 0;
2108  mapp[mappI[0]] = 0;
2109  for(int i=1;i<nchare;i++){
2110  if(index[i]!=index[j]){j++; index[j]=index[i];}
2111  mapp[mappI[i]] = j;
2112  }//endfor
2113 
2114  nval_out[0] = j+1;
2115  nmax_out[0] = nmax;
2116 }
2117 ///////////////////////////////////////////////////////////////////////////=
2118 
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
rfftwnd_complex_to_real_split
void rfftwnd_complex_to_real_split(RFFTplanHolder *rfftplanholder, int howmany, fftw_complex *in, int istride_in, int idist_in, fftw_real *out_in, int ostride_in, int odist_in, int split, int index)
split up an fft call into multiple invocations with cmiprogress calls between them.
Definition: fftCache.C:1736
UberCollection
holds the UberIndex and the offset for proxies
Definition: Uber.h:68
FFTcache::doStpFFTGtoR_Gchare
void doStpFFTGtoR_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int)
StatePlane : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward /////////////////////////////////////...
Definition: fftCache.C:1169
FFTcache::doEextFFTGxToRx_Rchare
void doEextFFTGxToRx_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward //////////////////////////////////////////////...
Definition: fftCache.C:1075
rfftwnd_real_to_complex_split
void rfftwnd_real_to_complex_split(RFFTplanHolder *rfftplanholder, int howmany, fftw_real *in, int istride_in, int idist_in, fftw_complex *out_in, int ostride_in, int odist_in, int split, int index)
split up an fft call into multiple invocations with cmiprogress calls between them.
Definition: fftCache.C:1827
FFTcache::doStpFFTRtoG_Rchare
void doStpFFTRtoG_Rchare(complex *, double *, int, int, int, int)
////////////////////////////////////////////////////////////////////////////= ///////////////////////...
Definition: fftCache.C:1288
FFTcache::doStpFFTRtoG_Gchare
void doStpFFTRtoG_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int)
StatePlane : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward ////////////////////////////////////...
Definition: fftCache.C:1135
FFTcache::doStpFFTGtoR_Rchare
void doStpFFTGtoR_Rchare(complex *, double *, int, int, int, int)
Definition: fftCache.C:1203
FFTcache::doEextFFTGtoR_Gchare
void doEextFFTGtoR_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int)
Eext : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward ///////////////////////////////////////////...
Definition: fftCache.C:890
FFTcache::doRhoFFTGtoR_Rchare
void doRhoFFTGtoR_Rchare(complex *, double *, int, int, int, int)
Definition: fftCache.C:1541
fft_split
void fft_split(FFTplanHolder *fftplanholder, int howmany, fftw_complex *in, int istride_in, int idist_in, fftw_complex *out_in, int ostride_in, int odist_in, int split, int index)
split up an fft call into multiple invocations with cmiprogress calls between them.
Definition: fftCache.C:1645
essl_work
= Holder classes for the plans : Allows many fft libaries to be used
Definition: fftCacheSlab.h:325
FFTcache::doEextFFTGtoR_Rchare
void doEextFFTGtoR_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward //////////////////////////////////////////////...
Definition: fftCache.C:1028
FFTcache::doRhoFFTRxToGx_Rchare
void doRhoFFTRxToGx_Rchare(complex *, double *, int, int, int, int)
Definition: fftCache.C:1482
FFTcache::doEextFFTRtoG_Rchare
void doEextFFTRtoG_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward /////////////////////////////////////////////...
Definition: fftCache.C:923
FFTcache::expandGSpace
void expandGSpace(complex *data, complex *packedData, RunDescriptor *runs, int numRuns, int numFull, int numPoints, int nfftz)
packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz ////////////////////////////...
Definition: fftCache.C:496
fftCacheSlab.h
Add type declarations for simulationConstants class (readonly vars) and once class for each type of o...
FFTcache::doNlFFTRtoG_Gchare
void doNlFFTRtoG_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int)
non-local : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward /////////////////////////////////////...
Definition: fftCache.C:598
FFTcache::doRhoFFTGyToRy_Rchare
void doRhoFFTGyToRy_Rchare(complex *, double *, int, int, int, int)
Definition: fftCache.C:1586
FFTcache::doEextFFTRtoG_Gchare
void doEextFFTRtoG_Gchare(complex *, int, int, int, int, RunDescriptor *, int, int)
Eext : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward //////////////////////////////////////////...
Definition: fftCache.C:855
FFTcache::doNlFFTGtoR_Rchare
void doNlFFTGtoR_Rchare(complex *, double *, int, int, int, int)
non-local : Rchare : data(gx,gy,z) -> data(x,y,z) : Forward /////////////////////////////////////////...
Definition: fftCache.C:771
FFTcache::FFTcache
FFTcache(int _ngrida, int _ngridb, int _ngridc, int _ngridaEext, int _ngridbEext, int _ngridcEext, int _ees_eext_on, int _ngridaNL, int _ngridbNL, int _ngridcNL, int _ees_NL_on, int _nlines_max, int _nlines_max_rho, int _nchareGState, int _nchareRState, int _nchareGNL, int _nchareRNL, int _nchareGRho, int _nchareRRho, int _nchareRRhoTot, int _nchareGEext, int _nchareREext, int _nchareREextTot, int *numGState, int *numRXState, int *numRYState, int *numRYStateLower, int *numGNL, int *numRXNL, int *numRYNL, int *numRYNLLower, int *numGRho, int *numRXRho, int *numRYRho, int *numGEext, int *numRXEext, int *numRYEext, int _fftopt, int _nsplitR, int _nsplitG, int _rhoRsubPlanes, UberCollection _thisInstance)
= FFTcache - sets up a fftcache on each processor
Definition: fftCache.C:39
FFTcache::packGSpace
void packGSpace(complex *data, complex *packedPlaneData, RunDescriptor *runs, int numRuns, int numFull, int numPoints, int nfftz)
packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz ////////////////////////////...
Definition: fftCache.C:553
RunDescriptor
== Index logic for lines of constant x,y in gspace.
Definition: RunDescriptor.h:22
cpaimd.h
Some basic data structures and the array map classes are defined here.
FFTcache::doRhoFFTRtoG_Gchare
void doRhoFFTRtoG_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int, int)
Rho Plane : Gchare : data(gx,gy,z) -> data(gx,gy,gz) : backward /////////////////////////////////////...
Definition: fftCache.C:1401
FFTcache::doRhoFFTGtoR_Gchare
void doRhoFFTGtoR_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int, int)
rho Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward //////////////////////////////////////////////...
Definition: fftCache.C:1365
FFTcache::doEextFFTRxToGx_Rchare
void doEextFFTRxToGx_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward /////////////////////////////////////////////...
Definition: fftCache.C:969
FFTcache::doEextFFTGyToRy_Rchare
void doEextFFTGyToRy_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward //////////////////////////////////////////////...
Definition: fftCache.C:1106
CP_State_Plane.h
debug_flags.h
Useful debugging flags.
FFTcache::doRhoFFTGxToRx_Rchare
void doRhoFFTGxToRx_Rchare(complex *, double *, int, int, int, int)
Definition: fftCache.C:1611
FFTcache::doRhoFFTRtoG_Rchare
void doRhoFFTRtoG_Rchare(complex *, double *, int, int, int, int)
Definition: fftCache.C:1437
FFTcache::doNlFFTGtoR_Gchare
void doNlFFTGtoR_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int)
non-local : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward //////////////////////////////////////...
Definition: fftCache.C:632
FFTcache::doEextFFTRyToGy_Rchare
void doEextFFTRyToGy_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward /////////////////////////////////////////////...
Definition: fftCache.C:1001
FFTcache::doHartFFTGtoR_Gchare
void doHartFFTGtoR_Gchare(complex *, complex *, int, int, int, int, RunDescriptor *, int, int)
Hartree : Gchare : data(gx,gy,gz) -> data(gx,gy,z) : forward ////////////////////////////////////////...
Definition: fftCache.C:664
FFTcache::doRhoFFTRyToGy_Rchare
void doRhoFFTRyToGy_Rchare(complex *, double *, int, int, int, int)
Definition: fftCache.C:1514