OpenAtom  Version1.5a
util.C
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////////==
2 ///////////////////////////////////////////////////////////////////////////////////cc
3 ///////////////////////////////////////////////////////////////////////////////////==
4 /** \file util.C
5  *
6  */
7 ///////////////////////////////////////////////////////////////////////////////////==
8 
9 #define CHARM_ON
10 #include "src_piny_physics_v1.0/include/class_defs/piny_constants.h"
11 #include "util.h"
12 #include "para_grp_parse.h"
13 #include "main/CPcharmParaInfoGrp.h"
14 #include "src_piny_physics_v1.0/friend_lib/proto_friend_lib_entry.h"
15 #include "src_mathlib/mathlib.h"
16 #include "src_piny_physics_v1.0/include/class_defs/allclass_gen.h"
17 #include "src_piny_physics_v1.0/include/class_defs/allclass_cp.h"
18 #include "src_piny_physics_v1.0/include/class_defs/PINY_INIT/PhysicsParamTrans.h"
19 #include <assert.h>
20 extern Config config;
21 extern int sizeX;
22 #if CMK_PROJECTIONS_USE_ZLIB
23 #include "zlib.h"
24 #endif
25 
26 #ifndef M_PI
27 #define M_PI 3.14159265358979323846
28 #endif
29 
30 ///////////////////////////////////////////////////////////////////////////////////==
31 ///////////////////////////////////////////////////////////////////////////////////c
32 ///////////////////////////////////////////////////////////////////////////////////==
33 
34 void make_rho_runs(CPcharmParaInfo *sim){
35 
36 ///////////////////////////////////////////////////////////////////////////////////==
37 
38  GENERAL_DATA *general_data = GENERAL_DATA::get();
39  CP *cp = CP::get();
40 
41 #include "src_piny_physics_v1.0/include/class_defs/allclass_strip_gen.h"
42 #include "src_piny_physics_v1.0/include/class_defs/allclass_strip_cp.h"
43 
44 ///////////////////////////////////////////////////////////////////////////////////==
45 /// set up the rho kvectors
46 
47  int *kx;
48  int *ky;
49  int *kz;
50  int nline_tot;
51  int nPacked;
52  int rhoRsubplanes = config.rhoRsubplanes;
53  int rhoGHelpers = config.rhoGHelpers;
54 
55  int sizeXEext = sim->ngrid_eext_a;
56  int sizeYEext = sim->ngrid_eext_b;
57  int sizeZEext = sim->ngrid_eext_c;
58  int sizeX = sim->sizeX;
59  int sizeY = sim->sizeY;
60  int sizeZ = sim->sizeZ;
61 
62  double *hmati = gencell->hmati;
63  double ecut4 = 8.0*cpcoeffs_info->ecut; // convert to Ryd extra factor of 2.0
64 
65  int ka_max = 2*(sim->kx_max);
66  int kb_max = 2*(sim->ky_max);
67  int kc_max = 2*(sim->kz_max);
68 
69  get_rho_kvectors(ecut4,hmati,&kx,&ky,&kz,&nline_tot,&nPacked,ka_max,kb_max,kc_max);
70 
71 ///////////////////////////////////////////////////////////////////////////////////==
72 /// Reorder the kvectors to produce better balance for the lines :
73 
74  int *kx_ind = new int[nline_tot];
75  int *kx_line = new int[nline_tot];
76  int *ky_line = new int[nline_tot];
77  int *kx_line_ext = new int[nline_tot];
78  int *ky_line_ext = new int[nline_tot];
79  int *istrt_line = new int [nline_tot];
80  int *iend_line = new int [nline_tot];
81  int *npts_line = new int [nline_tot];
82 
83  int nplane_x=0;
84  int ic = 0;
85  istrt_line[0] = 0;
86  for(int i = 1;i<nPacked;i++){
87  int iii = abs(kx[i]);
88  nplane_x = (iii > nplane_x ? iii : nplane_x);
89  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[(i-1)]){
90  iend_line[ic] = i;
91  npts_line[ic] = iend_line[ic]-istrt_line[ic];
92  ic++;
93  istrt_line[ic] = i;
94  }//endfor
95  }//endfor
96  iend_line[ic] = nPacked;
97  npts_line[ic] = iend_line[ic]-istrt_line[ic];
98  ic++;
99  if(ic!=nline_tot){
100  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
101  CkPrintf("Toasty Rho Line Flip-lines!\n");
102  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
103  CkExit();
104  }//endif
105  nplane_x += 1;
106 
107  double temp = ((double)nplane_x)*config.gExpandFactRho;
108  int nchareRhoG = (int)temp;
109  int nx = sim->sizeX;
110  int ny = sim->sizeY;
111  int nz = sim->sizeZ;
112  int nx_ext = sim->ngrid_eext_a;
113  int ny_ext = sim->ngrid_eext_b;
114  int nz_ext = sim->ngrid_eext_c;
115 
116  int *kxt = new int[nPacked];
117  int *kyt = new int[nPacked];
118  int *kzt = new int[nPacked];
119  CmiMemcpy(kxt,kx,(nPacked*sizeof(int)));
120  CmiMemcpy(kyt,ky,(nPacked*sizeof(int)));
121  CmiMemcpy(kzt,kz,(nPacked*sizeof(int)));
122 
123  int *mapl = new int[nline_tot];
124  for(int i=0;i<nline_tot;i++){mapl[i]=i;}
125 
126  if(config.rhoLineOrder==0){
127  int nsplit = (3*nplane_x)/2;
128  int jj=0;
129  for(int i=0;i<nsplit; i++){
130  for(int j=i;j<nline_tot;j+=nsplit){
131  mapl[jj]=j; jj++;
132  }}//endfor
133  }//end switch
134  if(config.rhoLineOrder==1){
135  long seed = 174571;
136  for(int j=0;j<4;j++){
137  for(int i=0;i<nline_tot;i++){
138  double stuff = altRandom(&seed);
139  int index = (int)( ((double) nline_tot)*stuff );
140  index = MIN(index,nline_tot-1);
141  int itemp = mapl[i];
142  mapl[i] = mapl[index];
143  mapl[index] = itemp;
144  }}//endfor
145  }//endif
146 
147  int jc = 0;
148  for(int i=0;i<nline_tot;i++){
149  int j = mapl[i];
150  for(int lt=istrt_line[j],l=jc;lt<iend_line[j];lt++,l++){
151  kx[l] = kxt[lt];
152  ky[l] = kyt[lt];
153  kz[l] = kzt[lt];
154  }//endfor
155  jc+=npts_line[j];
156  }//endfor
157  if(jc!=nPacked) {
158  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
159  CkPrintf("Toasty Rho Line Flip-pts!\n");
160  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
161  CkExit();
162  }//endif
163 
164  delete [] mapl;
165  delete [] kxt;
166  delete [] kyt;
167  delete [] kzt;
168 
169  ic = 0;
170  istrt_line[0] = 0;
171  kx_line[ic] = kx[ic];
172  ky_line[ic] = ky[ic];
173  kx_line_ext[ic] = kx[ic];
174  ky_line_ext[ic] = ky[ic];
175  if(kx_line[ic] <0){kx_line[ic] +=nx;}
176  if(ky_line[ic] <0){ky_line[ic] +=ny;}
177  if(kx_line_ext[ic]<0){kx_line_ext[ic]+=nx_ext;}
178  if(ky_line_ext[ic]<0){ky_line_ext[ic]+=ny_ext;}
179  for(int i = 1;i<nPacked;i++){
180  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[(i-1)]){
181  iend_line[ic] = i;
182  npts_line[ic] = iend_line[ic]-istrt_line[ic];
183  ic++;
184  istrt_line[ic] = i;
185  kx_line[ic] = kx[i];
186  ky_line[ic] = ky[i];
187  kx_line_ext[ic] = kx[i];
188  ky_line_ext[ic] = ky[i];
189  if(kx_line[ic] <0){kx_line[ic] +=nx;}
190  if(ky_line[ic] <0){ky_line[ic] +=ny;}
191  if(kx_line_ext[ic]<0){kx_line_ext[ic]+=nx_ext;}
192  if(ky_line_ext[ic]<0){ky_line_ext[ic]+=ny_ext;}
193  }//endfor
194  }//endfor
195  iend_line[ic] = nPacked;
196  npts_line[ic] = iend_line[ic]-istrt_line[ic];
197  ic++;
198 
199  if(ic!=nline_tot){
200  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
201  CkPrintf("Toasty Line Rho Flip-lines.b!\n");
202  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
203  CkExit();
204  }//endif
205 
206 ///////////////////////////////////////////////////////////////////////////////////==
207 /// Create runs using reorder k-vectors
208 
209  int nrun_tot = 1;
210  int run_length_sum = 0;
211  int run_length = 1;
212  int curr_x = kx[0];
213  int curr_y = ky[0];
214  int curr_z = kz[0];
215  if(curr_x<0){curr_x+=nx;}
216  if(curr_y<0){curr_y+=ny;}
217  if(curr_z<0){curr_z+=nz;}
218  int tmpz = curr_z;
219  int nline_tot_now = 1;
220 
221  CkVec<RunDescriptor> runs;
222 
223  if(kz[0]>=0){ // add neg half-line of size 0
224  runs.push_back(RunDescriptor(curr_x,curr_y,0,run_length_sum,0,1,nz));
225  nrun_tot +=1;
226  }//endif
227  for(int pNo=1;pNo<nPacked;pNo++) {
228  int x = kx[pNo];
229  int y = ky[pNo];
230  int z = kz[pNo];
231  if (x<0){x+=nx;}
232  if (y<0){y+=ny;}
233  if (z<0){z+=nz;}
234  // Count the lines of z by twos
235  if (z == tmpz + 1 && x == curr_x && y == curr_y) {
236  // same half line : keep counting
237  run_length++;
238  tmpz += 1;
239  }else{
240  // We have changed half lines so increment run index
241  // Each line of z, constant x,y is stored in 2 run descriptors
242  // Example : -3 -2 -1 is a separate ``run of z''
243  // 0 1 2 3 is a separte ``run of z''
244  // for a line with only a 0 add a zero length descriptor
245  // to represent the missing negative part of the line.
246  runs.push_back(RunDescriptor(curr_x,curr_y,curr_z,run_length_sum,run_length,1,nz));
247  nrun_tot +=1;
248  run_length_sum += run_length;
249  if((curr_x!=x || curr_y!=y) && kz[pNo-1]<0){ // add pos half-line of size 0
250  runs.push_back(RunDescriptor(curr_x,curr_y,0,run_length_sum,0,1,nz));
251  nrun_tot +=1;
252  }//endif
253  if((curr_x!=x || curr_y!=y) && kz[pNo]>=0){ // add neg half-line of size 0
254  runs.push_back(RunDescriptor(x,y,0,run_length_sum,0,1,nz));
255  nrun_tot +=1;
256  }//endif
257  curr_x = x;
258  curr_y = y;
259  curr_z = z;
260  tmpz = z;
261  run_length = 1;
262  }//endif
263  if(kx[pNo]!=kx[(pNo-1)] || ky[pNo]!=ky[(pNo-1)] ){
264  nline_tot_now++;
265  if( (nrun_tot-1)/2 != nline_tot_now-1){
266  CkPrintf("\n");
267  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
268  CkPrintf("Broken Rho Run Desscriptor : %d %d %d : %d %d %d: %d %d\n",
269  kx[pNo],ky[pNo],kz[pNo],kx[(pNo-1)],ky[(pNo-1)],kz[(pNo-1)],
270  nrun_tot-1,nline_tot_now-1);
271  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
272  CkExit();
273  }//endif
274  }//endif
275  }//endfor
276  // Record the last run of z.
277  runs.push_back(RunDescriptor(curr_x,curr_y,curr_z,run_length_sum,run_length,1,nz));
278  run_length_sum += run_length;
279  if(kz[nPacked-1]<0){ // add pos half-line of size 0
280  runs.push_back(RunDescriptor(curr_x,curr_y,0,run_length_sum,0,1,nz));
281  nrun_tot +=1;
282  }//endif
283 
284 
285  if(run_length_sum!=nPacked){
286  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
287  CkPrintf("The rho rundescriptor didn't assign all pts to rho runs \n");
288  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
289  CkExit();
290  }//endif
291 
292  if((nrun_tot %2)!=0 || nrun_tot != 2*nline_tot){
293  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
294  CkPrintf("The rho rundescriptor MUST find an even number of half-lines\n");
295  CkPrintf("The rho rundescriptor MUST find the correct number of lines\n");
296  CkPrintf("%d %d %d %d\n",nrun_tot,nrun_tot/2,nline_tot,
297  nline_tot_now);
298  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
299  CkExit();
300  }//endif
301 
302  for(int i=0;i<nline_tot;i+=2){
303  RunDescriptor *Desi = &runs[i];
304  RunDescriptor *Desi1 = &runs[(i+1)];
305  //CkPrintf("i %d kx %d kx1 %d ky %d ky1 %d\n",i,Desi->x,Desi1->x,Desi->y,Desi1->y);
306  if( (Desi->x != Desi1->x) || (Desi->y != Desi1->y) ){
307  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
308  CkPrintf("The rho rundescriptor MUST pair up the half-lines\n");
309  CkPrintf("i %d kx %d kx1 %d ky %d ky1 %d\n",i,Desi->x,Desi1->x,Desi->y,Desi1->y);
310  CkPrintf("or you will not be a happy camper :\n");
311  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
312  CkExit();
313  }//endif
314  if(Desi->z==0 && Desi->length != 0){
315  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
316  CkPrintf("The rho rundescriptor with 1st z == 0 must have 0 lngth\n");
317  CkPrintf("%d %d %d : %d %d %d\n",Desi->x,Desi->y,Desi->z,
318  Desi1->x,Desi1->y,Desi1->z);
319  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
320  CkExit();
321  }//endif
322  if(Desi1->z<0 || (Desi->z<=nz/2 && Desi->length != 0)){
323  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
324  CkPrintf("The rho rundescriptor MUST have 1st z <=0 and 2nd >=0\n");
325  CkPrintf("%d %d %d : %d %d %d\n",Desi->x,Desi->y,Desi->z,
326  Desi1->x,Desi1->y,Desi1->z);
327  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
328  CkExit();
329  }//endif
330  }//endfor
331 
332 ///////////////////////////////////////////////////////////////////////////////////==
333 /// Decompose lines to balance points
334 
335  int nchareRhoGEext = nchareRhoG*rhoGHelpers;
336  int *istrt_lgrp = new int [nchareRhoG];
337  int *iend_lgrp = new int [nchareRhoG];
338  int *npts_lgrp = new int [nchareRhoG];
339  int *nline_lgrp = new int [nchareRhoG];
340  int *nline_lgrp_eext = new int [nchareRhoGEext];
341 
342  ParaGrpParse::get_chareG_line_prms(nPacked,nchareRhoG,nline_tot,npts_line,
343  istrt_lgrp,iend_lgrp,npts_lgrp,nline_lgrp,false);
344 
345 //////////////////////////////////////////////////////////////////////////////
346 /// Create the line decomposition and a sorted run descriptor
347 /// There are two rundescriptors per line : Noah's arc sort
348 
349  int nlines_min=nline_lgrp[0];
350  int nlines_max=0;
351  for(int i=0;i<nchareRhoG;i++){
352  nlines_max=MAX(nlines_max,nline_lgrp[i]);
353  nlines_min=MIN(nlines_min,nline_lgrp[i]);
354  }//endfor
355  int **index_tran_upack_rho = cmall_int_mat(0,nchareRhoG,0,nlines_max,"util.C");
356  int **index_tran_upack_eext = cmall_int_mat(0,nchareRhoGEext,0,nlines_max,"util.C");
357 
358  int yspace=sizeX/2+1;
359  for(int igrp=0;igrp<nchareRhoG;igrp++){
360  for(int i=istrt_lgrp[igrp],j=0;i<iend_lgrp[igrp];i++,j++){
361  index_tran_upack_rho[igrp][j] = kx_line[i] + ky_line[i]*yspace;
362  }//endfor
363  }//endfor
364 
365  int nlines_max_eext = 0;
366  int yspaceEext=sizeXEext/2+1;
367  for(int igrp=0,jgrp=0;igrp<nchareRhoG;igrp++){
368  int nlTot = nline_lgrp[igrp];
369  int istrt = istrt_lgrp[igrp];
370  for(int k=0;k<rhoGHelpers;k++,jgrp++){
371  int kstrt,kend,nl;
372  getSplitDecomp(&kstrt,&kend,&nl,nlTot,rhoGHelpers,k);
373  kstrt += istrt;
374  nline_lgrp_eext[jgrp] = nl;
375  nlines_max_eext = MAX(nlines_max_eext,nl);
376  for(int j=0,i=kstrt;j<nl;j++,i++){
377  index_tran_upack_eext[jgrp][j] = kx_line_ext[i] + ky_line_ext[i]*yspaceEext;
378  }//endfor
379  }//endfor
380  }//endfor
381 
382  CkVec<RunDescriptor> *RhosortedRunDescriptors;
383  RhosortedRunDescriptors = new CkVec<RunDescriptor> [nchareRhoG];
384  for(int igrp = 0; igrp < nchareRhoG; igrp++){
385  for(int i=istrt_lgrp[igrp];i<iend_lgrp[igrp];i++){
386  int j = 2*i;
387  int j1 = 2*i+1;
388  RhosortedRunDescriptors[igrp].push_back(runs[j]);
389  RhosortedRunDescriptors[igrp].push_back(runs[j1]);
390  }//endfor
391  }//endfor
392 
393  for(int x = 0; x < nchareRhoG; x ++) {
394  int runsToBeSent = RhosortedRunDescriptors[x].size();
395  int numPoints = 0;
396  for (int j = 0; j < RhosortedRunDescriptors[x].size(); j++){
397  numPoints += RhosortedRunDescriptors[x][j].length;
398  }//endfor
399  }//endfor
400 
401 //////////////////////////////////////////////////////////////////////////////
402 /// variables that could be used for mapping but aren't yet.
403 
404  double *pts_per_chare = new double[nchareRhoG];
405  double *lines_per_chare = new double[nchareRhoG];
406  for(int i=0;i<nchareRhoG;i++){
407  pts_per_chare[i] =(double) npts_lgrp[i];
408  lines_per_chare[i]=(double) nline_lgrp[i];
409  }//endfor
410 
411 //////////////////////////////////////////////////////////////////////////////
412 /// Check for rhoghelper size
413 
414  if(nlines_min==0){
415  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
416  CkPrintf("No lines in a RhoG collection. Your RhoG decomp stinks.\n");
417  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
418  CkExit();
419  }//endif
420 
421  if(nlines_min<config.rhoGHelpers){
422  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
423  CkPrintf("RhoGHelper parameter, %d, must be <= minimum number\n",config.rhoGHelpers);
424  CkPrintf("of lines in any RhoG chare array element, %d.\n",nlines_min);
425  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
426  CkExit();
427  }//endif
428 
429  if(rhoRsubplanes > nplane_x){
430  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
431  CkPrintf("RhoGsubplanes parameter, %d, must be <= number\n",rhoRsubplanes);
432  CkPrintf("number of gx values which span 0<=gx <=%d.\n",nplane_x);
433  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
434  CkExit();
435  }//endif
436 
437 //////////////////////////////////////////////////////////////////////////////
438 /// Rejigger the stuff for subplane parallelization : RHOG
439 
440  int ***index_tran_upack_rho_y = NULL;
441  int ***index_tran_pack_rho_y = NULL;
442  int **nline_send_rho_y = NULL;
443  int ***index_tran_upack_eext_y = NULL;
444  int ***index_tran_upack_eext_ys= NULL;
445  int ***index_tran_pack_eext_y = NULL;
446  int ***index_tran_pack_eext_ys = NULL;
447  int **nline_send_eext_y = NULL;
448  int *numSubGx = NULL;
449  int **listSubGx = NULL;
450  int listSubFlag = 0;
451  int ngxSubMax = nplane_x/rhoRsubplanes+1;
452  if(rhoRsubplanes==1){ngxSubMax=0;}
453 
454  if(rhoRsubplanes>1){
455  CkPrintf("\n");
456  PRINT_LINE_STAR;
457  CkPrintf("Creating the subPlane maps\n");
458  PRINT_LINE_DASH;
459  }//endif
460 
461  if(rhoRsubplanes>1){//subplanes on
462  //----------------------------------------------------------------
463  // for normal rho work
464  index_tran_upack_rho_y = cmall_itens3(0,nchareRhoG,0,rhoRsubplanes,
465  0,nlines_max,"util.C");
466  index_tran_pack_rho_y = cmall_itens3(0,nchareRhoG,0,rhoRsubplanes,
467  0,nlines_max,"util.C");
468  nline_send_rho_y = cmall_int_mat(0,nchareRhoG,0,rhoRsubplanes,"util.C");
469  //----------------------------------------------------------------
470  // for eext-ees rho work
471  index_tran_upack_eext_y = cmall_itens3(0,nchareRhoGEext,0,rhoRsubplanes,
472  0,nlines_max_eext,"util.C");
473  index_tran_upack_eext_ys = cmall_itens3(0,nchareRhoGEext,0,rhoRsubplanes,
474  0,nlines_max_eext,"util.C");
475  index_tran_pack_eext_y = cmall_itens3(0,nchareRhoGEext,0,rhoRsubplanes,
476  0,nlines_max_eext,"util.C");
477  index_tran_pack_eext_ys = cmall_itens3(0,nchareRhoGEext,0,rhoRsubplanes,
478  0,nlines_max_eext,"util.C");
479  nline_send_eext_y = cmall_int_mat(0,nchareRhoGEext,0,rhoRsubplanes,"util.C");
480  //----------------------------------------------------------------
481  // Balanced SubPlane decomp
482  numSubGx = new int [rhoRsubplanes];
483  listSubGx = cmall_int_mat(0,rhoRsubplanes,0,ngxSubMax,"util.C");
484  listSubFlag = 0;
485 
486  //----------------------------------------------------------------
487  // Group the Gx in the subplanes so as to minimize # of message sent from RtoG
488  int *listGx = new int [nplane_x];
489  int *mapGrpGx = new int [nplane_x];
490  int *mapMemGx = new int [nplane_x];
491 
492  int div = (nplane_x/rhoRsubplanes);
493  int rem = (nplane_x % rhoRsubplanes);
494  for(int igrp=0;igrp<rhoRsubplanes;igrp++){
495  int max = (igrp < rem ? 1 : 0);
496  numSubGx[igrp] = max+div;
497  }//endfor
498  create_subPlane_decomp(nplane_x,listGx,mapGrpGx,nchareRhoGEext,numSubGx,
499  nline_lgrp_eext,kx_line,nline_send_eext_y,rhoRsubplanes);
500  if(config.rhoSubPlaneBalance==1){
501  for(int i=0;i<nplane_x;i++){mapGrpGx[listGx[i]] = i;}
502  create_gx_decomp(nline_tot,nplane_x,kx_line,mapGrpGx,rhoRsubplanes,numSubGx);
503  }//endfif
504 
505  listSubFlag=0;
506  for(int i=0;i<nplane_x;i++){
507  if(listGx[i]!=i){listSubFlag==1;}
508  }//endif
509  if(listSubFlag==0){
510  CkPrintf("This is a straight run through gx on the subplanes\n");
511  }//endif
512  int iii = 0;
513  for(int igrp=0;igrp<rhoRsubplanes;igrp++){
514  int num = numSubGx[igrp];
515  if(num>1){sort_me(num,&listGx[iii]);} //order the gx you have
516  CkPrintf("subplane[%d] : %d : Gx { ",igrp,num);
517  for(int jc=0,ic=iii;ic<iii+num;ic++,jc++){
518  listSubGx[igrp][jc] = listGx[ic];
519  mapGrpGx[listGx[ic]] = igrp;
520  mapMemGx[listGx[ic]] = jc;
521  CkPrintf("%d ",listGx[ic]);
522  }//endfor
523  CkPrintf("}\n");
524  iii += num;
525  }//endfor
526 
527  //----------------------------------------------------------------
528  // RhoR(gx,gy,z) parallelized by gx(subPlane) and z
529  // RhoG(gx,gy,z) parallelized by collections of {gx,gy}
530  // pack and upack indicies for rhoG(gx,gy,z) <-> rhoR(gx,gy,z)
531  for(int igrp=0,i=0;igrp<nchareRhoG;igrp++){
532  for(int ic=0;ic<rhoRsubplanes;ic++){
533  nline_send_rho_y[igrp][ic]=0;
534  }}//endfor
535 
536  for(int igrp=0;igrp<nchareRhoG;igrp++){
537  for(int i=istrt_lgrp[igrp],j=0;i<iend_lgrp[igrp];i++,j++){
538  int ic = mapGrpGx[kx_line[i]]; //subPlane index where kx is located
539  int ip = mapMemGx[kx_line[i]]; //which kx in the subplane this is
540  int jc = nline_send_rho_y[igrp][ic]; ///nt pts sent from subplane to g-chare
541  nline_send_rho_y[igrp][ic]++;
542  index_tran_pack_rho_y[igrp][ic][jc] = j*sizeZ;
543  index_tran_upack_rho_y[igrp][ic][jc] = ip*sizeY + ky_line[i];
544  }//endfor
545  }//endfor
546 
547  int nsend_min=10000000;
548  int nsend_max=0;
549  for(int igrp=0;igrp<nchareRhoG;igrp++){
550  for(int ic=0;ic<rhoRsubplanes;ic++){
551  if(nline_send_rho_y[igrp][ic]>0){
552  nsend_min= MIN(nsend_min, nline_send_rho_y[igrp][ic]);
553  nsend_max= MAX(nsend_max, nline_send_rho_y[igrp][ic]);
554  }//endif
555  }//endif
556  }//endfor
557  CkPrintf("Msg size Imbalance : rhoG <-> rhoR min %d max %d\n",
558  nsend_min, nsend_max);
559 
560  //----------------------------------------------------------------
561  // EextR(gx,gy,z) parallelized by gx(subPlane) and z : bigger Z than rho
562  // EextG(gx,gy,z) parallelized by collections of {gx,gy} (subdivided)
563  // pack and upack indicies for EextG(gx,gy,z) <-> EextR(gx,gy,z)
564  for(int igrp=0,i=0;igrp<nchareRhoGEext;igrp++){
565  for(int ic=0;ic<rhoRsubplanes;ic++){
566  nline_send_eext_y[igrp][ic]=0;
567  }}//endfor
568 
569  for(int igrp=0,i=0;igrp<nchareRhoGEext;igrp++){
570  for(int j=0;j<nline_lgrp_eext[igrp];j++,i++){
571  int ic = mapGrpGx[kx_line[i]]; //subPlane index
572  int ip = mapMemGx[kx_line[i]]; //which kx in the subplane
573  int jc = nline_send_eext_y[igrp][ic]; // cnt num pts sent
574  nline_send_eext_y[igrp][ic]++;
575  index_tran_pack_eext_y[igrp][ic][jc] = j*sizeZEext;
576  index_tran_pack_eext_ys[igrp][ic][jc] = j*sizeZ;
577  index_tran_upack_eext_y[igrp][ic][jc] = ip*sizeYEext + ky_line_ext[i];
578  index_tran_upack_eext_ys[igrp][ic][jc]= ip*sizeY + ky_line[i];
579  }//endfor
580  }//endfor
581 
582  nsend_min=10000000;
583  nsend_max=0;
584  for(int igrp=0;igrp<nchareRhoGEext;igrp++){
585  for(int ic=0;ic<rhoRsubplanes;ic++){
586  if(nline_send_eext_y[igrp][ic]>0){
587  nsend_min= MIN(nsend_min, nline_send_eext_y[igrp][ic]);
588  nsend_max= MAX(nsend_max, nline_send_eext_y[igrp][ic]);
589  }//endif
590  }//endfor
591  }//endfor
592  CkPrintf("Msg size Imbalance : rhoGhart <-> rhoRhart min %d max %d\n",
593  nsend_min, nsend_max);
594 
595  delete [] listGx;
596  delete [] mapGrpGx;
597  delete [] mapMemGx;
598 
599  }//endif : subplanes are on.
600 
601 //////////////////////////////////////////////////////////////////////////////
602 /// Pack up the stuff
603 
604  config.nchareRhoG = nchareRhoG;
605 
606  sim->nplane_rho_x = nplane_x;
607  sim->nchareRhoG = nchareRhoG;
608  sim->nchareRhoGEext = nchareRhoGEext;
609  sim->npts_per_chareRhoG = npts_lgrp;
610  sim->index_tran_upack_rho = index_tran_upack_rho;
611  sim->index_tran_upack_eext = index_tran_upack_eext;
612  sim->nlines_max_rho = nlines_max;
613  sim->nlines_per_chareRhoG = nline_lgrp;
614  sim->nlines_per_chareRhoGEext= nline_lgrp_eext;
615  sim->RhosortedRunDescriptors = RhosortedRunDescriptors;
616  sim->npts_tot_rho = nPacked;
617  sim->nlines_tot_rho = nline_tot;
618  sim->lines_per_chareRhoG = lines_per_chare;
619  sim->pts_per_chareRhoG = pts_per_chare;
620 
621  sim->rhoRsubplanes = config.rhoRsubplanes;
622  sim->nlines_max_eext = nlines_max_eext;
623  sim->nline_send_eext_y = nline_send_eext_y;
624  sim->nline_send_rho_y = nline_send_rho_y;
625  sim->index_tran_upack_rho_y = index_tran_upack_rho_y;
626  sim->index_tran_upack_eext_y = index_tran_upack_eext_y;
627  sim->index_tran_upack_eext_ys= index_tran_upack_eext_ys;
628  sim->index_tran_pack_rho_y = index_tran_pack_rho_y;
629  sim->index_tran_pack_eext_y = index_tran_pack_eext_y;
630  sim->index_tran_pack_eext_ys = index_tran_pack_eext_ys;
631 
632  sim->ngxSubMax = ngxSubMax; // max number of gx values in any grp
633  sim->numSubGx = numSubGx; // number of gx values in each grp
634  sim->listSubGx = listSubGx; // gx values in grp
635  sim->listSubFlag = listSubFlag;
636 
637 ///////////////////////////////////////////////////////////////////////////////////
638 /// Analyze the Send and Receives
639 
640  if(rhoRsubplanes>1){
641  int Rhart_max=0;
642  int Rhart_min=10000000;
643  double total=0;
644  for( int j=0; j< nchareRhoGEext;j++){
645  int recvCountFromRHartExt = 0;
646  for(int i=0;i<rhoRsubplanes;i++){
647  if(sim->nline_send_eext_y[j][i]>0)
648  recvCountFromRHartExt++;
649  }//endfor
650  recvCountFromRHartExt*=sizeZEext;
651  total += recvCountFromRHartExt;
652  Rhart_max=MAX(Rhart_max,recvCountFromRHartExt);
653  Rhart_min=MIN(Rhart_min,recvCountFromRHartExt);
654  }//endfor
655  total /= (double)nchareRhoGEext;
656  CkPrintf("GHart recv %d min msg %d max msg avg %g from RHart\n",
657  Rhart_min, Rhart_max,total);
658  }//endif
659 
660  if(rhoRsubplanes>1){
661  int Rho_max=0;
662  int Rho_min=10000000;
663  double total=0;
664  for( int j=0; j< nchareRhoGEext;j++){
665  int recvCountFromRho = 0;
666  for(int i=0;i<rhoRsubplanes;i++){
667  if(sim->nline_send_eext_y[j][i]>0)
668  recvCountFromRho++;
669  }//endfor
670  recvCountFromRho*=sizeZ;
671  total += recvCountFromRho;
672  Rho_max=MAX(Rho_max,recvCountFromRho);
673  Rho_min=MIN(Rho_min,recvCountFromRho);
674  }//endfor
675  total /= (double)nchareRhoGEext;
676  CkPrintf("GHart recv %d min msg %d max msg avg %g from RRho\n",
677  Rho_min, Rho_max,total);
678  }//endif
679 
680  if(rhoRsubplanes>1){
681  int RRho_max=0;
682  int RRho_min=10000000;
683  double total=0;
684  for( int j=0; j< nchareRhoG;j++){
685  int recvCountFromRRho = 0;
686  for(int i=0;i<rhoRsubplanes;i++){
687  if(sim->nline_send_rho_y[j][i]>0)
688  recvCountFromRRho++;
689  }//endfor
690  recvCountFromRRho*=sizeZ;
691  total += recvCountFromRRho;
692  RRho_max=MAX(RRho_max,recvCountFromRRho);
693  RRho_min=MIN(RRho_min,recvCountFromRRho);
694  }//endfor
695  total /= (double)nchareRhoG;
696  CkPrintf("GRho recv %d min msg %d max msg avg %g from RRho\n",
697  RRho_min,RRho_max,total);
698  }//endif
699 
700  if(rhoRsubplanes>1){
701  int GRho_max=0;
702  int GRho_min=10000000;
703  for( int j=0; j< rhoRsubplanes;j++){
704  int recvCountFromGRho = 0;
705  for(int i=0;i<nchareRhoG;i++){
706  if(sim->nline_send_rho_y[i][j]>0)
707  recvCountFromGRho++;
708  }//endfor
709  GRho_max=MAX(GRho_max,recvCountFromGRho);
710  GRho_min=MIN(GRho_min,recvCountFromGRho);
711  }//endfor
712  CkPrintf("RRho recv %d min msg %d max msg from RhoG\n",GRho_min, GRho_max);
713  }//endif
714 
715  if(rhoRsubplanes>1){
716  int GHart_max=0;
717  int GHart_min=10000000;
718  for( int j=0; j< rhoRsubplanes;j++){
719  int recvCountFromGHartExt = 0;
720  for(int i=0;i<nchareRhoGEext;i++){
721  if(sim->nline_send_eext_y[i][j]>0)
722  recvCountFromGHartExt++;
723  }//endfor
724  GHart_max=MAX(GHart_max,recvCountFromGHartExt);
725  GHart_min=MIN(GHart_min,recvCountFromGHartExt);
726  }//endfor
727  CkPrintf("RRho recv %d min msg %d max msg from GHart\n",GHart_min, GHart_max);
728  }//endif
729 
730  if(rhoRsubplanes>1){
731  int GHart_max=0;
732  int GHart_min=10000000;
733  for( int j=0; j< rhoRsubplanes;j++){
734  int recvCountFromGHartExt = 0;
735  for(int i=0;i<nchareRhoGEext;i++){
736  if(sim->nline_send_eext_y[i][j]>0)
737  recvCountFromGHartExt++;
738  }//endfor
739  GHart_max=MAX(GHart_max,recvCountFromGHartExt);
740  GHart_min=MIN(GHart_min,recvCountFromGHartExt);
741  }//endfor
742  CkPrintf("RHart recv %d min msg %d max msg from GHart\n",GHart_min, GHart_max);
743  }//endif
744 
745  if(rhoRsubplanes>1){
746  PRINT_LINE_DASH;
747  CkPrintf("Completed subPlane map creation.\n");
748  PRINT_LINE_STAR; CkPrintf("\n\n");
749  }//endif
750 
751 //////////////////////////////////////////////////////////////////////////////
752 /// Clean up the memory
753 
754  delete [] kx_ind;
755  delete [] kx_line;
756  delete [] ky_line;
757  delete [] kx_line_ext;
758  delete [] ky_line_ext;
759  delete [] istrt_line;
760  delete [] iend_line;
761  delete [] npts_line;
762  delete [] kx;
763  delete [] ky;
764  delete [] kz;
765  delete [] istrt_lgrp;
766  delete [] iend_lgrp;
767 
768 //////////////////////////////////////////////////////////////////////////////
769  }//end routine
770 //////////////////////////////////////////////////////////////////////////////
771 
772 
773 //////////////////////////////////////////////////////////////////////////////
774 //////////////////////////////////////////////////////////////////////////////
775 //////////////////////////////////////////////////////////////////////////////
776 
777 void get_rho_kvectors(double ecut4, double *hmati, int **kx_ret, int **ky_ret,
778  int **kz_ret, int *nline_tot_ret, int *nPacked_ret,
779  int ka_max, int kb_max, int kc_max)
780 
781 //////////////////////////////////////////////////////////////////////////////
782  {//begin routine
783 //////////////////////////////////////////////////////////////////////////////
784 /// count the k-vectors : preserve nice symmetry for non-cubic boxes even
785 /// though you need more kvectors
786 
787  int iii;
788  double tpi = 2.0*M_PI;
789 
790  int nPacked = 0;
791  int nPacked_np = 0;
792  for(int ka=0;ka<=ka_max;ka++){
793  for(int kb=-kb_max;kb<=kb_max;kb++){
794  for(int kc=-kc_max;kc<=kc_max;kc++){
795  double gx = tpi*(ka*hmati[1] + kb*hmati[2] + kc*hmati[3]);
796  double gy = tpi*(ka*hmati[4] + kb*hmati[5] + kc*hmati[6]);
797  double gz = tpi*(ka*hmati[7] + kb*hmati[8] + kc*hmati[9]);
798  double g2 = gx*gx+gy*gy+gz*gz;
799  if(g2<=ecut4){nPacked++;}
800  if(g2<=ecut4){
801  nPacked_np++;
802  if(ka==0 && kb<0){nPacked_np--;}
803  if(ka==0 && kb==0 && kc<=0){nPacked_np--;}
804  }
805  }//endfor:kc
806  }//endfor:kb
807  }//endfor:ka
808 
809  CkPrintf("Rho kvectors perfect sphere %d : %d %d %d\n",nPacked_np,ka_max,kb_max,kc_max);
810 
811 //////////////////////////////////////////////////////////////////////////////
812 /// fill the k-vectors
813 
814  int *kx = new int[nPacked];
815  int *ky = new int[nPacked];
816  int *kz = new int[nPacked];
817 
818  int ic = 0;
819  for(int ka=0;ka<=ka_max;ka++){
820  for(int kb=-kb_max;kb<=kb_max;kb++){
821  for(int kc=-kc_max;kc<=kc_max;kc++){
822  double gx = tpi*(ka*hmati[1] + kb*hmati[2] + kc*hmati[3]);
823  double gy = tpi*(ka*hmati[4] + kb*hmati[5] + kc*hmati[6]);
824  double gz = tpi*(ka*hmati[7] + kb*hmati[8] + kc*hmati[9]);
825  double g2 = gx*gx+gy*gy+gz*gz;
826  if(g2<=ecut4){
827  kx[ic]=ka;
828  ky[ic]=kb;
829  kz[ic]=kc;
830  ic++;
831  }/*endif*/
832  }//endfor
833  }//endfor
834  }//endfor
835 
836 //////////////////////////////////////////////////////////////////////////////
837 /// Count the lines
838 
839  int nline_tot = 1;
840  for(int i=1; i<nPacked;i++){
841  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[i-1]){nline_tot++;}
842  }//endfor
843 
844 //////////////////////////////////////////////////////////////////////////////
845 /// Set return values
846 
847  *kx_ret = kx;
848  *ky_ret = ky;
849  *kz_ret = kz;
850  *nline_tot_ret = nline_tot;
851  *nPacked_ret = nPacked;
852 
853 //////////////////////////////////////////////////////////////////////////////
854  }//end routine
855 //////////////////////////////////////////////////////////////////////////////
856 
857 
858 ///////////////////////////////////////////////////////////////////////////////////==
859 ///////////////////////////////////////////////////////////////////////////////////cc
860 ///////////////////////////////////////////////////////////////////////////////////==
861 
862 void readStateIntoRuns(int nPacked, int ncoef,complex *arrCP, CkVec<RunDescriptor> &runs,
863  const char *fromFile,int ibinary_opt,
864  int *nline_tot_ret,int *nplane_ret,
865  int *istrt_lgrp,int *iend_lgrp,
866  int *npts_lgrp,int *nline_lgrp,
867  int **kx_line_ret, int **ky_line_ret,
868  int **kx_ret,int **ky_ret, int **kz_ret, int iget_decomp,
869  int gen_wave,int nx_in, int ny_in , int nz_in)
870 
871 ///////////////////////////////////////////////////////////////////////////////////==
872  {//begin routine
873 ///////////////////////////////////////////////////////////////////////////////////==
874 /// A little screen output for the fans
875 
876 #ifdef _CP_DEBUG_UTIL_VERBOSE_
877  CkPrintf("Reading state from file: %s\n",fromFile);
878 #endif
879  if(ibinary_opt < 0 || ibinary_opt > 3){
880  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
881  CkPrintf("Bad binary option : %d %s\n",ibinary_opt,fromFile);
882  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
883  CkExit();
884  }//endif
885 
886 ///////////////////////////////////////////////////////////////////////////////////==
887 /// First read in the state and k-vectors : allows parsing of doublePack option
888 
889  int doublePack = config.doublePack;
890  int nx,ny,nz;
891  int *kx= new int [nPacked];
892  int *ky= new int [nPacked];
893  int *kz= new int [nPacked];
894 
895  if(gen_wave==0){
896  readState(nPacked,arrCP,fromFile,ibinary_opt,nline_tot_ret,
897  nplane_ret,kx,ky,kz,&nx,&ny,&nz,
898  istrt_lgrp,iend_lgrp,npts_lgrp,nline_lgrp,iget_decomp,0);
899  }else{
900  nx = nx_in;
901  ny = ny_in;
902  nz = nz_in;
903  kx -= 1; ky -= 1; kz -=1;
904  PhysicsParamTransfer::fetch_state_kvecs(kx,ky,kz,ncoef,doublePack);
905  kx += 1; ky += 1; kz +=1;
906  processState(nPacked,ncoef,arrCP,fromFile,ibinary_opt,nline_tot_ret,
907  nplane_ret,kx,ky,kz,
908  istrt_lgrp,iend_lgrp,npts_lgrp,nline_lgrp,iget_decomp,0,0,ny);
909  }//endif
910 
911  int nplane = (*nplane_ret);
912  int nline_tot = (*nline_tot_ret);
913  int nchareG = config.nchareG;
914 
915  if(nx!=nx_in || ny!=ny_in || nz!=nz_in){
916  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
917  CkPrintf("Bad State fft size Input\n");
918  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
919  CkExit();
920  }//endif
921 
922 ///////////////////////////////////////////////////////////////////////////////////==
923 /// Read the state into the rundescriptor puppy dog
924 
925 #ifdef _INPUT_FOR_KPTS_MESSES_UP_
926  if(!config.doublePack){
927  CkPrintf("@@@@@@@@@@@@@@@@@@@@_warning_@@@@@@@@@@@@@@@@@@@@\n");
928  CkPrintf("The rundescriptor needs some love for the non-double pack\n");
929  CkPrintf("It is not consistent with new FFT logic due to input data order\n");
930  CkPrintf("If the data is just reordered all should be well, %s\n",fromFile);
931  CkPrintf("Raz is on the job.\n");
932  CkPrintf("@@@@@@@@@@@@@@@@@@@@_warning_@@@@@@@@@@@@@@@@@@@@\n");
933  }//endif
934 #endif
935 
936  int nrun_tot = 1;
937  int run_length_sum = 0;
938  int run_length = 1;
939  int curr_x = kx[0];
940  int curr_y = ky[0];
941  int curr_z = kz[0];
942  if(curr_x<0){curr_x+=nx;}
943  if(curr_y<0){curr_y+=ny;}
944  if(curr_z<0){curr_z+=nz;}
945  int tmpz = curr_z;
946  int nline_tot_now = 1;
947 
948  if(kz[0]>=0){ // add neg half-line of size 0
949  runs.push_back(RunDescriptor(curr_x,curr_y,curr_z,run_length_sum,0,1,nz));
950  nrun_tot +=1;
951  }//endif
952 
953  for(int pNo=1;pNo<nPacked;pNo++) {
954  int x = kx[pNo];
955  int y = ky[pNo];
956  int z = kz[pNo];
957  if (x<0){x+=nx;}
958  if (y<0){y+=ny;}
959  if (z<0){z+=nz;}
960  // Count the lines of z by twos
961  if (z == tmpz + 1 && x == curr_x && y == curr_y) {
962  // same half line : keep counting
963  run_length++;
964  tmpz += 1;
965  }else{
966  // We have changed half lines so increment run index
967  // Each line of z, constant x,y is stored in 2 run descriptors
968  // Example : -3 -2 -1 is a separate ``run of z''
969  // 0 1 2 3 is a separte ``run of z''
970  // for a line with only a 0 add a zero length descriptor
971  // to represent the missing negative part of the line.
972  runs.push_back(RunDescriptor(curr_x,curr_y,curr_z,run_length_sum,run_length,1,nz));
973  nrun_tot +=1;
974  run_length_sum += run_length;
975  if((curr_x!=x || curr_y!=y) && kz[pNo-1]<0){ // add pos half-line of size 0
976  runs.push_back(RunDescriptor(curr_x,curr_y,0,run_length_sum,0,1,nz));
977  nrun_tot +=1;
978  }//endif
979  if((curr_x!=x || curr_y!=y) && kz[pNo]>=0){ // add neg half-line of size 0
980  runs.push_back(RunDescriptor(x,y,0,run_length_sum,0,1,nz));
981  nrun_tot +=1;
982 
983  }//endif
984  curr_x = x;
985  curr_y = y;
986  curr_z = z;
987  tmpz = z;
988  run_length = 1;
989  }//endif
990  if(kx[pNo]!=kx[(pNo-1)] || ky[pNo]!=ky[(pNo-1)] ){
991  nline_tot_now++;
992  if( (nrun_tot-1)/2 != nline_tot_now-1){
993  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
994  CkPrintf("Broken State Run Desscriptor : %d %d %d : %d %d %d: %d %d\n",
995  kx[pNo],ky[pNo],kz[pNo],kx[(pNo-1)],ky[(pNo-1)],kz[(pNo-1)],
996  nrun_tot-1,nline_tot_now-1);
997  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
998  CkExit();
999  }//endif
1000  }//endif
1001  }//endfor
1002  // Record the last run of z.
1003  runs.push_back(RunDescriptor(curr_x,curr_y,curr_z,run_length_sum,run_length,1,nz));
1004  run_length_sum += run_length;
1005  if(kz[nPacked-1]<0){ // add pos half-line of size 0
1006  runs.push_back(RunDescriptor(curr_x,curr_y,0,run_length_sum,0,1,nz));
1007  nrun_tot +=1;
1008  }//endif
1009 
1010  if(run_length_sum!=nPacked){
1011  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1012  CkPrintf("The state rundescriptor didn't assign all pts to runs %s\n",fromFile);
1013  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1014  CkExit();
1015  }//endif
1016 
1017  if((nrun_tot %2)!=0 || nrun_tot != 2*nline_tot){
1018  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1019  CkPrintf("The state rundescriptor MUST find an even number of half-lines\n");
1020  CkPrintf("The rundescriptor MUST find the correct number of lines\n");
1021  CkPrintf("%d %d %d %d %s\n",nrun_tot,nrun_tot/2,nline_tot,
1022  nline_tot_now,fromFile);
1023  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1024  CkExit();
1025  }//endif
1026 
1027  for(int i=0;i<nline_tot;i+=2){
1028  RunDescriptor *Desi = &runs[i];
1029  RunDescriptor *Desi1 = &runs[(i+1)];
1030  if( (Desi->x != Desi1->x) || (Desi->y != Desi1->y) ){
1031  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1032  CkPrintf("The state rundescriptor MUST pair up the half-lines\n");
1033  CkPrintf("or you will not be a happy camper : %s\n",fromFile);
1034  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1035  CkExit();
1036  }//endif
1037  if(Desi->z==0 && Desi->length != 0){
1038  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1039  CkPrintf("The state rundescriptor with 1st z == 0 must have 0 lngth\n");
1040  CkPrintf("%d %d %d : %d %d %d\n",Desi->x,Desi->y,Desi->z,
1041  Desi1->x,Desi1->y,Desi1->z);
1042  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1043  CkExit();
1044  }//endif
1045  if(Desi1->z<0 || (Desi->z<=nz/2 && Desi->length != 0)){
1046  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1047  CkPrintf("The state rundescriptor MUST have 1st z <=0 and 2nd >=0\n");
1048  CkPrintf("%d %d %d : %d %d %d\n",Desi->x,Desi->y,Desi->z,
1049  Desi1->x,Desi1->y,Desi1->z);
1050  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1051  CkExit();
1052  }//endif
1053  }//endfor
1054 
1055 ///////////////////////////////////////////////////////////////////////////////////==
1056 
1057  int *kx_line = new int [nline_tot];
1058  int *ky_line = new int [nline_tot];
1059  int ic = 0;
1060  kx_line[0] = kx[0];
1061  ky_line[0] = ky[0];
1062  if(kx_line[0]<0){kx_line[0]+=nx;}
1063  if(ky_line[0]<0){ky_line[0]+=ny;}
1064 
1065  for(int i = 1;i<nPacked;i++){
1066  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[(i-1)]){
1067  ic++;
1068  kx_line[ic] = kx[i];
1069  ky_line[ic] = ky[i];
1070  if(ky_line[ic]<0){ky_line[ic]+=ny;}
1071  if(kx_line[ic]<0){kx_line[ic]+=nx;}
1072  }//endfor
1073  }//endfor
1074  ic++;
1075 
1076  if(ic!=nline_tot){
1077  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1078  CkPrintf("Incorrect number of (state) lines : %s\n",fromFile);
1079  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1080  CkExit();
1081  }//endif
1082 
1083  *kx_line_ret = kx_line;
1084  *ky_line_ret = ky_line;
1085 
1086 ///////////////////////////////////////////////////////////////////////////////////==
1087 
1088  *kx_ret = kx;
1089  *ky_ret = ky;
1090  *kz_ret = kz;
1091 
1092 ///////////////////////////////////////////////////////////////////////////////////==
1093 /// A little output to the screen!
1094 
1095 #ifdef _CP_DEBUG_UTIL_VERBOSE_
1096  CkPrintf("Done reading state from file: %s\n",fromFile);
1097 #endif
1098 
1099 //----------------------------------------------------------------------------------
1100  }//end routine
1101 ///////////////////////////////////////////////////////////////////////////////////==
1102 
1103 
1104 ///////////////////////////////////////////////////////////////////////////////////==
1105 ///////////////////////////////////////////////////////////////////////////////////cc
1106 ///////////////////////////////////////////////////////////////////////////////////==
1107 
1108 void readState(int nPacked, complex *arrCP, const char *fromFile,int ibinary_opt,
1109  int *nline_tot_ret,int *nplane_ret, int *kx, int *ky, int *kz,
1110  int *nx_ret, int *ny_ret, int *nz_ret,
1111  int *istrt_lgrp,int *iend_lgrp,int *npts_lgrp,int *nline_lgrp,
1112  int iget_decomp,int iget_vstate)
1113 
1114 ///////////////////////////////////////////////////////////////////////////////////==
1115  {//begin routine
1116 ///////////////////////////////////////////////////////////////////////////////////==
1117 /// A little screen output for the fans
1118 
1119  char stuff[25];
1120  if(iget_vstate==0){
1121  strcpy(stuff,"coef");
1122  }else{
1123  strcpy(stuff,"velocity");
1124  }
1125 #ifdef _CP_DEBUG_UTIL_VERBOSE_
1126  CkPrintf("Reading %s state file: %s\n",stuff,fromFile);
1127 #endif
1128  if(ibinary_opt < 0 || ibinary_opt > 3){
1129  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1130  CkPrintf("Bad binary option : %d\n",ibinary_opt);
1131  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1132  CkExit();
1133  }//endif
1134 
1135 ///////////////////////////////////////////////////////////////////////////////////==
1136 /// First read in the state and k-vectors : allows parsing of doublePack option
1137 
1138 #if !CMK_PROJECTIONS_USE_ZLIB
1139  if(ibinary_opt>1)
1140  {
1141  CkPrintf("Attempt to use ZLIB Failed! Please review compilation\n");
1142  //CkPrintf("Macro cmk-projections-use-zlib is %d \n", CMK_PROJECTIONS_USE_ZLIB);
1143  CkExit();
1144  }
1145 #endif
1146  int nx,ny,nz;
1147  int nktot = 0;
1148  if(ibinary_opt==0){
1149  FILE *fp=fopen(fromFile,"r");
1150  if (fp==NULL){
1151  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1152  CkPrintf("Can't open %s state file %s\n",stuff,fromFile);
1153  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1154  CkExit();
1155  }//endif
1156  int nPackedLoc;
1157  if(4!=fscanf(fp,"%d%d%d%d",&nPackedLoc,&nx,&ny,&nz)){
1158  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1159  CkPrintf("Can't parse size line of file %s\n",fromFile);
1160  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1161  CkExit();
1162  }//endif
1163  for(int pNo=0;pNo<nPacked;pNo++) {
1164  int x,y,z;
1165  double re,im;
1166  if(5!=fscanf(fp,"%lg%lg%d%d%d",&re,&im,&x,&y,&z)){
1167  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1168  CkPrintf("Can't parse packed %s state location %s\n",stuff,fromFile);
1169  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1170  CkExit();
1171  }//endif
1172  arrCP[pNo] = complex(re, im);
1173  kx[pNo] = x;
1174  ky[pNo] = y;
1175  kz[pNo] = z;
1176  nktot++;
1177  if(config.doublePack && x==0 && y==0 && z==0){break;}
1178  }//endfor
1179  fclose(fp);
1180  }
1181 #if CMK_PROJECTIONS_USE_ZLIB
1182  else if(ibinary_opt==2){
1183  // CkPrintf("Using ZLIB to load ascii states\n");
1184  char bigenough[1000]; //we know our lines are shorter than this
1185  char localFile[1000]; // fromFile is const
1186  strcpy(localFile,fromFile);
1187  strcat(localFile,".gz");
1188  gzFile zfp=gzopen(localFile,"rb");
1189  if (zfp==NULL){
1190  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1191  CkPrintf("Can't open %s state file %s\n",stuff,localFile);
1192  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1193  CkExit();
1194  }//endif
1195  int nPackedLoc;
1196  if(gzgets(zfp,bigenough,1000)!=Z_NULL)
1197  {
1198  if(4!=sscanf(bigenough,"%d%d%d%d",&nPackedLoc,&nx,&ny,&nz)){
1199  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1200  CkPrintf("Can't parse size line of file %s\n",localFile);
1201  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1202  CkExit();
1203  }//endif
1204  }
1205  else
1206  {
1207  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1208  CkPrintf("Can't parse size line of file %s\n",localFile);
1209  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1210  CkExit();
1211  }
1212  for(int pNo=0;pNo<nPacked;pNo++) {
1213  int x,y,z;
1214  double re,im;
1215  if(gzgets(zfp,bigenough,1000)!=Z_NULL)
1216  {
1217  if(5!=sscanf(bigenough,"%lg%lg%d%d%d",&re,&im,&x,&y,&z)){
1218  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1219  CkPrintf("Can't parse packed %s state location %s\n",stuff,localFile);
1220  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1221  CkExit();
1222  }//endif
1223  }
1224  else
1225  {
1226  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1227  CkPrintf("Can't parse packed %s state location %s\n",stuff,localFile);
1228  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1229  CkExit();
1230  }
1231  arrCP[pNo] = complex(re, im);
1232  kx[pNo] = x;
1233  ky[pNo] = y;
1234  kz[pNo] = z;
1235  nktot++;
1236  if(config.doublePack && x==0 && y==0 && z==0){break;}
1237  }//endfor
1238  gzclose(zfp);
1239  }
1240  else if (ibinary_opt==3)
1241  {
1242  // CkPrintf("Using ZLIB to load binary states\n");
1243  char localFile[1000]; // fromFile is const
1244  strcpy(localFile,fromFile);
1245  strcat(localFile,".gz");
1246  gzFile zfp=gzopen(localFile,"rb");
1247  if (zfp==NULL){
1248  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1249  CkPrintf("Can't open %s state file %s\n",stuff,localFile);
1250  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1251  CkExit();
1252  }
1253  int nPackedLoc;
1254  int n=1;
1255  gzread(zfp,&(nPackedLoc),sizeof(int));
1256  gzread(zfp,&(nx),sizeof(int));
1257  gzread(zfp,&(ny),sizeof(int));
1258  gzread(zfp,&(nz),sizeof(int));
1259  for(int pNo=0;pNo<nPacked;pNo++) {
1260  int x,y,z;
1261  double re,im;
1262  gzread(zfp,&(re),sizeof(double));
1263  gzread(zfp,&(im),sizeof(double));
1264  gzread(zfp,&(x),sizeof(int));
1265  gzread(zfp,&(y),sizeof(int));
1266  gzread(zfp,&(z),sizeof(int));
1267  arrCP[pNo] = complex(re, im);
1268  kx[pNo] = x;
1269  ky[pNo] = y;
1270  kz[pNo] = z;
1271  nktot++;
1272  if(config.doublePack && x==0 && y==0 && z==0){break;}
1273  }
1274  gzclose(zfp);
1275 
1276  }
1277 #endif
1278  else{
1279  FILE *fp=fopen(fromFile,"rb");
1280  if (fp==NULL){
1281  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1282  CkPrintf("Can't open %s state file %s\n",stuff,fromFile);
1283  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1284  CkExit();
1285  }
1286  int nPackedLoc;
1287  int n=1;
1288  assert(fread(&(nPackedLoc),sizeof(int),n,fp)>0);
1289  assert(fread(&(nx),sizeof(int),n,fp));
1290  assert(fread(&(ny),sizeof(int),n,fp));
1291  assert(fread(&(nz),sizeof(int),n,fp));
1292  for(int pNo=0;pNo<nPacked;pNo++) {
1293  int x,y,z;
1294  double re,im;
1295  assert(fread(&(re),sizeof(double),n,fp));
1296  assert(fread(&(im),sizeof(double),n,fp));
1297  assert(fread(&(x),sizeof(int),n,fp));
1298  assert(fread(&(y),sizeof(int),n,fp));
1299  assert(fread(&(z),sizeof(int),n,fp));
1300  arrCP[pNo] = complex(re, im);
1301  kx[pNo] = x;
1302  ky[pNo] = y;
1303  kz[pNo] = z;
1304  nktot++;
1305  if(config.doublePack && x==0 && y==0 && z==0){break;}
1306  }//endfor
1307  fclose(fp);
1308  }//endif
1309 
1310 #ifdef _CP_DEBUG_UTIL_VERBOSE_
1311  CkPrintf("Done reading %s state from file: %s\n",stuff,fromFile);
1312 #endif
1313 
1314 ///////////////////////////////////////////////////////////////////////////////////==
1315 
1316  processState(nPacked,nktot,arrCP,fromFile,ibinary_opt,nline_tot_ret,nplane_ret,kx,ky,kz,
1317  istrt_lgrp,iend_lgrp,npts_lgrp,nline_lgrp,iget_decomp,iget_vstate,1,ny);
1318 
1319  *nx_ret = nx;
1320  *ny_ret = ny;
1321  *nz_ret = nz;
1322 
1323 //----------------------------------------------------------------------------------
1324  }//end routine
1325 ///////////////////////////////////////////////////////////////////////////////////==
1326 
1327 
1328 ///////////////////////////////////////////////////////////////////////////////////==
1329 ///////////////////////////////////////////////////////////////////////////////////cc
1330 ///////////////////////////////////////////////////////////////////////////////////==
1331 
1332 void processState(int nPacked, int nktot, complex *arrCP, const char *fromFile,int ibinary_opt,
1333  int *nline_tot_ret,int *nplane_ret, int *kx, int *ky, int *kz,
1334  int *istrt_lgrp,int *iend_lgrp,int *npts_lgrp,int *nline_lgrp,
1335  int iget_decomp,int iget_vstate, int iopt,int ny) {
1336 
1337 ///////////////////////////////////////////////////////////////////////////////////==
1338 /// Set up some variables
1339 
1340 #ifdef _CP_DEBUG_UTIL_VERBOSE_
1341  CkPrintf("Proceesing %s state %s\n",stuff,fromFile);
1342 #endif
1343 
1344  int nchareG = config.nchareG;
1345 
1346  char stuff[25];
1347  if(iget_vstate==0){
1348  strcpy(stuff,"coef");
1349  }else{
1350  strcpy(stuff,"velocity");
1351  }/*endif*/
1352 
1353 ///////////////////////////////////////////////////////////////////////////////////==
1354 /// Add the bottom half of plane zero because the code likes to have it.
1355 
1356  if(config.doublePack){
1357  int n_ret;
1358  ParaGrpParse::flip_data_set(nktot,&n_ret,kx,ky,kz,arrCP,iopt);
1359  if(n_ret!=nPacked){
1360  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1361  CkPrintf("Bad num pts in readState() %s: %d %d \n",nktot,n_ret,fromFile);
1362  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1363  CkExit();
1364  }//endif
1365  }//endif
1366 
1367  // We are hoping here that the sensible input for !doublePack works corrrectly
1368 
1369 ///////////////////////////////////////////////////////////////////////////////////==
1370 /// Process the state data
1371 
1372  int nline_tot = 1;
1373  int nplane = 1;
1374  for(int i = 1;i<nPacked;i++){
1375  if(kx[i]!=kx[(i-1)]){nplane++;}
1376  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[(i-1)]){nline_tot++;}
1377  if(kx[i]<kx[(i-1)]){
1378  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1379  CkPrintf("Bad x-flip in processState() %s\n",fromFile);
1380  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1381  CkExit();
1382  }//endif
1383  if(kx[i]==kx[(i-1)]){
1384  if(ky[i]<ky[(i-1)]){
1385  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1386  CkPrintf("Bad y-flip in readState() %s\n",fromFile);
1387  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1388  CkExit();
1389  }//endif
1390  }//endif
1391  if(kx[i]==kx[(i-1)] && ky[i]==ky[(i-1)]){
1392  if(kz[i]!=(kz[(i-1)]+1)){
1393  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1394  CkPrintf("Bad y-flip in readState() %s\n",fromFile);
1395  CkPrintf(" %d %d %d\n",kx[i],ky[i],kz[i]);
1396  CkPrintf(" %d %d %d\n",kx[(i-1)],ky[(i-1)],kz[(i-1)]);
1397  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1398  CkExit();
1399  }//endif
1400  }//endif
1401  }//endfor : pts in g-space
1402 
1403  if(config.nchareG>nline_tot){
1404  CkPrintf("@@@@@@@@@@@@@@@@@@@@_Error_@@@@@@@@@@@@@@@@@@@@\n");
1405  CkPrintf("Dude, too much juice on the chunks. Chill on gExpandFact\n");
1406  CkPrintf("@@@@@@@@@@@@@@@@@@@@_Error_@@@@@@@@@@@@@@@@@@@@\n");
1407  CkExit();
1408  }//endif
1409 
1410 ///////////////////////////////////////////////////////////////////////////////////==
1411 /// Reorder the data to produce better balance for the lines :
1412 
1413  int *kx_ind = new int[nline_tot];
1414  int *kx_line = new int[nline_tot];
1415  int *ky_line = new int[nline_tot];
1416  int *istrt_line = new int [nline_tot];
1417  int *iend_line = new int [nline_tot];
1418  int *npts_line = new int [nline_tot];
1419  int ic = 0;
1420  istrt_line[0] = 0;
1421  for(int i = 1;i<nPacked;i++){
1422  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[(i-1)]){
1423  iend_line[ic] = i;
1424  npts_line[ic] = iend_line[ic]-istrt_line[ic];
1425  ic++;
1426  istrt_line[ic] = i;
1427  }//endfor
1428  }//endfor
1429  iend_line[ic] = nPacked;
1430  npts_line[ic] = iend_line[ic]-istrt_line[ic];
1431  ic++;
1432  if(ic!=nline_tot){
1433  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1434  CkPrintf("Toasty Line Flip-lines!\n");
1435  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1436  CkExit();
1437  }//endif
1438 
1439  complex *arrCPt;
1440  if(iopt==1){
1441  arrCPt = new complex[nPacked];
1442  CmiMemcpy(arrCPt,arrCP,(nPacked*sizeof(complex)));
1443  }//endif
1444 
1445  int *kxt = new int[nPacked];
1446  int *kyt = new int[nPacked];
1447  int *kzt = new int[nPacked];
1448  CmiMemcpy(kxt,kx,(nPacked*sizeof(int)));
1449  CmiMemcpy(kyt,ky,(nPacked*sizeof(int)));
1450  CmiMemcpy(kzt,kz,(nPacked*sizeof(int)));
1451 
1452  int jc = 0;
1453  int lc = 0;
1454  for(int i=0;i<nchareG; i++){
1455  for(int j=i;j<nline_tot;j+=nchareG){
1456  for(int lt=istrt_line[j],l=jc;lt<iend_line[j];lt++,l++){
1457  kx[l] = kxt[lt];
1458  ky[l] = kyt[lt];
1459  kz[l] = kzt[lt];
1460  if(iopt==1){arrCP[l] = arrCPt[lt];}
1461  }//endfor
1462  jc+=npts_line[j];
1463  lc++;
1464  }//endfor
1465  }//endfor
1466  if(jc!=nPacked) {
1467  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1468  CkPrintf("Toasty Line Flip-pts!\n");
1469  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1470  CkExit();
1471  }//endif
1472  if(lc!=nline_tot){
1473  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1474  CkPrintf("Toasty Line Flip-lines!\n");
1475  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1476  CkExit();
1477  }//endif
1478 
1479  ic = 0;
1480  istrt_line[0] = 0;
1481  int kxmax = kx[0];
1482  kx_line[0] = kx[0];
1483  ky_line[0] = ky[0];
1484  for(int i = 1;i<nPacked;i++){
1485  kxmax = MAX(kx[i],kxmax);
1486  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[(i-1)]){
1487  iend_line[ic] = i;
1488  npts_line[ic] = iend_line[ic]-istrt_line[ic];
1489  ic++;
1490  istrt_line[ic] = i;
1491  kx_line[ic] = kx[i];
1492  ky_line[ic] = ky[i];
1493  }//endfor
1494  }//endfor
1495  iend_line[ic] = nPacked;
1496  npts_line[ic] = iend_line[ic]-istrt_line[ic];
1497  ic++;
1498  if(ic!=nline_tot){
1499  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1500  CkPrintf("Toasty Line Flip-lines.b!\n");
1501  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1502  CkExit();
1503  }//endif
1504 
1505 ///////////////////////////////////////////////////////////////////////////////////==
1506 /// Decompose if necessary : note istrt_lgrp, iend_lgrp only define when iget_decomp==1
1507 
1508  if(iget_decomp==1){
1509  if(istrt_lgrp==NULL || iend_lgrp == NULL){
1510  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1511  CkPrintf("Toasty Line Flip memory!\n");
1512  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1513  CkExit();
1514  }//endif
1515  ParaGrpParse::get_chareG_line_prms(nPacked,nchareG,nline_tot,npts_line,
1516  istrt_lgrp,iend_lgrp,npts_lgrp,nline_lgrp,true);
1517  }//endif
1518 
1519 ///////////////////////////////////////////////////////////////////////////////////==
1520 /// For each decomposd chunk : sort on kx
1521 /// note istrt_lgrp, iend_lgrp only define when iget_decomp==1
1522 
1523  int mal_size = MAX(nline_tot,(kxmax+1));
1524  int *kx_tmp = new int[nline_tot];
1525  int *ky_tmp = new int[nline_tot];
1526  int *k_tmp = new int[mal_size];
1527  if(iopt==1){CmiMemcpy(arrCPt,arrCP,(nPacked*sizeof(complex)));}
1528  CmiMemcpy(kxt,kx,(nPacked*sizeof(int)));
1529  CmiMemcpy(kyt,ky,(nPacked*sizeof(int)));
1530  CmiMemcpy(kzt,kz,(nPacked*sizeof(int)));
1531 
1532  int loff = 0;
1533  int joff = 0;
1534  for(int i=0;i<nchareG;i++){
1535  for(int l=0;l<nline_lgrp[i];l++){
1536  kx_tmp[l] = kx_line[(l+loff)];
1537  ky_tmp[l] = ky_line[(l+loff)];
1538  kx_ind[l] = l;
1539  }//endfor
1540  if(nline_lgrp[i]>1){sort_kxky(nline_lgrp[i],kx_tmp,ky_tmp,kx_ind,k_tmp,ny);}
1541  for(int l=0;l<nline_lgrp[i];l++){
1542  int istrt = istrt_line[(kx_ind[l]+loff)];
1543  int iend = iend_line[(kx_ind[l]+loff)];
1544  for(int j=istrt,jk=joff;j<iend;j++,jk++){
1545  if(iopt==1){arrCP[jk] = arrCPt[j];}
1546  kx[jk] = kxt[j];
1547  ky[jk] = kyt[j];
1548  kz[jk] = kzt[j];
1549  }//endfor
1550  joff += npts_line[(kx_ind[l]+loff)];
1551  }//endfor
1552  loff += nline_lgrp[i];
1553  }//endfor
1554  if(joff!=nPacked) {
1555  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1556  CkPrintf("Toasty Line Flip-pts.2!\n");
1557  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1558  CkExit();
1559  }//endif
1560  if(loff!=nline_tot){
1561  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1562  CkPrintf("Toasty Line Flip-lines.2!\n");
1563  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1564  CkExit();
1565  }//endif
1566 
1567 ///////////////////////////////////////////////////////////////////////////////////==
1568 /// Debug output
1569 
1570 #ifdef _CP_DEBUG_LINE_
1571  if(iopt==1 && doublePack==1){
1572  double norm = 0;
1573  for(int i=1;i<nPacked;i++){
1574  double wght_now = 2.0;
1575  if(kx[i]==0 && ky[i]<0){wght_now=0.0;}
1576  if(kx[i]==0 && ky[i]==0 && kz[i]<0){wght_now=0.0;}
1577  if(kx[i]==0 && ky[i]==0 && kz[i]==0){wght_now=1.0;}
1578  norm += (wght_now)*arrCP[i].getMagSqr();
1579  }//endfor
1580  double normt = 0;
1581  for(int i=1;i<nPacked;i++){
1582  double wght_now = 2.0;
1583  if(kxt[i]==0 && kyt[i]<0){wght_now=0.0;}
1584  if(kxt[i]==0 && kyt[i]==0 && kzt[i]<0){wght_now=0.0;}
1585  if(kxt[i]==0 && kyt[i]==0 && kzt[i]==0){wght_now=1.0;}
1586  normt += (wght_now)*arrCPt[i].getMagSqr();
1587  }//endif
1588  CkPrintf("state : %g %g\n",norm,normt);
1589  }//endif
1590 
1591  if(iopt==1 && doublePack==0){
1592  double norm = 0;
1593  for(int i=1;i<nPacked;i++){
1594  norm += arrCP[i].getMagSqr();
1595  }//endfor
1596  double normt = 0;
1597  for(int i=1;i<nPacked;i++){
1598  normt += arrCPt[i].getMagSqr();
1599  }//endif
1600  CkPrintf("state : %g %g\n",norm,normt);
1601  }//endif
1602 #endif
1603 
1604 ///////////////////////////////////////////////////////////////////////////////////==
1605 /// Fix !double pack
1606 
1607 #ifdef _INPUT_FOR_KPTS_MESSES_UP_
1608  if(!config.doublePack){
1609  CkPrintf("@@@@@@@@@@@@@@@@@@@@_warning_@@@@@@@@@@@@@@@@@@@@\n");
1610  CkPrintf("The rundescriptor needs some love for the non-double pack\n");
1611  CkPrintf("It is not consistent with new FFT logic due to input data order\n");
1612  CkPrintf("If the data is just reordered all should be well, %s\n",fromFile);
1613  CkPrintf("Raz is on the job.\n");
1614  CkPrintf("@@@@@@@@@@@@@@@@@@@@_warning_@@@@@@@@@@@@@@@@@@@@\n");
1615  }//endif
1616 #endif
1617 
1618 ///////////////////////////////////////////////////////////////////////////////////==
1619 /// A little output to the screen!
1620 
1621  *nplane_ret = nplane;
1622  *nline_tot_ret = nline_tot;
1623 
1624  if(iopt==1){delete [] arrCPt;}
1625  delete [] kxt;
1626  delete [] kyt;
1627  delete [] kzt;
1628  delete [] istrt_line;
1629  delete [] iend_line;
1630  delete [] npts_line;
1631  delete [] kx_line;
1632  delete [] ky_line;
1633  delete [] kx_tmp;
1634  delete [] ky_tmp;
1635  delete [] k_tmp;
1636  delete [] kx_ind;
1637 
1638 #ifdef _CP_DEBUG_UTIL_VERBOSE_
1639  CkPrintf("Done proceesing %s state %s\n",stuff,fromFile);
1640 #endif
1641 
1642 //----------------------------------------------------------------------------------
1643  }//end routine
1644 ///////////////////////////////////////////////////////////////////////////////////==
1645 
1646 
1647 //////////////////////////////////////////////////////////////////////////////
1648 /// /////////////////////////////////////////////////////////////////////////cc
1649 //////////////////////////////////////////////////////////////////////////////
1650 
1651 void create_line_decomp_descriptor(CPcharmParaInfo *sim)
1652 
1653 //////////////////////////////////////////////////////////////////////////////
1654  { //begin rotunie
1655 //////////////////////////////////////////////////////////////////////////////
1656 /// Set the file name and data points
1657 
1658  char fname[1000];
1659  sprintf (fname, "%s/Spin.0_Kpt.0_Bead.0_Temper.0/state1.out",config.dataPath);
1660 
1661  int numData = config.numData;
1662  int ibinary_opt = sim->ibinary_opt;
1663  int sizeY = sim->sizeY;
1664  int sizeZ = sim->sizeZ;
1665  int nchareG = sim->nchareG;
1666  int sizeXNL = sim->ngrid_nloc_a;
1667  int nxNL = sim->ngrid_nloc_a;
1668  int nyNL = sim->ngrid_nloc_b;
1669  int gen_wave = sim->gen_wave;
1670  int ncoef = sim->ncoef;
1671  int doublePack = config.doublePack;
1672  double gExpandFact = config.gExpandFact;
1673 
1674 
1675 //////////////////////////////////////////////////////////////////////////////
1676 /// Get the complex data, Psi(g) and the run descriptor (z-lines in g-space)
1677 
1678  complex *complexPoints = new complex[numData];
1679  CkVec<RunDescriptor> runDescriptorVec;
1680  int nline_tot;
1681  int *istrt_lgrp = new int [nchareG];
1682  int *iend_lgrp = new int [nchareG];
1683  int *npts_lgrp = new int [nchareG];
1684  int *nline_lgrp = new int [nchareG];
1685  int *kx_line = NULL;
1686  int *ky_line = NULL;
1687  int *kx = NULL;
1688  int *ky = NULL;
1689  int *kz = NULL;
1690 
1691  readStateIntoRuns(numData,ncoef,complexPoints,runDescriptorVec,fname,ibinary_opt,
1692  &nline_tot,&(sim->nplane_x),istrt_lgrp,iend_lgrp,
1693  npts_lgrp,nline_lgrp,&kx_line,&ky_line,&kx,&ky,&kz,1,gen_wave,
1694  sizeX,sizeY,sizeZ);
1695  int nplane = sim->nplane_x;
1696 
1697  if(config.nGplane_x != nplane && config.doublePack){
1698  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1699  CkPrintf("Mismatch in allowed gspace planes %d %d\n",config.nGplane_x,nplane);
1700  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1701  CkExit();
1702  }//endif
1703 
1704  double temp = ((double)nplane)*gExpandFact;
1705  int mychareG = (int)temp;
1706  if(mychareG!=nchareG || mychareG!=config.nchareG){
1707  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1708  CkPrintf("Mismatch in allowed gspace chare arrays %d %d %d\n",
1709  mychareG,nchareG,config.nchareG);
1710  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1711  CkExit();
1712  }//endif
1713 
1714  int kx_min = kx[0];
1715  int kx_max = kx[0];
1716  int ky_min = ky[0];
1717  int ky_max = ky[0];
1718  int kz_min = kz[0];
1719  int kz_max = kz[0];
1720  for(int i = 1;i<numData;i++){
1721  kx_min = (kx_min <= kx[i] ? kx_min : kx[i]);
1722  ky_min = (ky_min <= ky[i] ? ky_min : ky[i]);
1723  kz_min = (kz_min <= kz[i] ? kz_min : kz[i]);
1724  kx_max = (kx_max >= kx[i] ? kx_max : kx[i]);
1725  ky_max = (ky_max >= ky[i] ? ky_max : ky[i]);
1726  kz_max = (kz_max >= kz[i] ? kz_max : kz[i]);
1727  }//endfor
1728  kx_max = (kx_max >= -kx_min ? kx_max : -kx_min);
1729  ky_max = (ky_max >= -ky_min ? ky_max : -ky_min);
1730  kz_max = (kz_max >= -kz_min ? kz_max : -kz_min);
1731 
1732  sim->kx_max = kx_max;
1733  sim->ky_max = ky_max;
1734  sim->kz_max = kz_max;
1735 
1736  if(2*kx_max >= sizeX/2 || 2*ky_max >= sizeY/2 || 2*kz_max >= sizeZ/2){
1737  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1738  CkPrintf("Bad FFT sizes : (%d,%d), (%d,%d), (%d,%d)\n",
1739  2*kx_max,sizeX/2,2*ky_max,sizeY/2,2*kz_max,sizeZ/2);
1740  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1741  CkExit();
1742  }//endif
1743 
1744 //////////////////////////////////////////////////////////////////////////////
1745 /// setup the non-local stuff
1746 
1747  int *kx_lineNL = new int[nline_tot];
1748  int *ky_lineNL = new int[nline_tot];
1749 
1750  int ic = 0;
1751  kx_lineNL[ic] = kx[ic];
1752  ky_lineNL[ic] = ky[ic];
1753  if(kx_lineNL[ic]<0){kx_lineNL[ic]+=nxNL;}
1754  if(ky_lineNL[ic]<0){ky_lineNL[ic]+=nyNL;}
1755  for(int i = 1;i<numData;i++){
1756  if(kx[i]!=kx[(i-1)] || ky[i]!=ky[(i-1)]){
1757  ic++;
1758  kx_lineNL[ic] = kx[i];
1759  ky_lineNL[ic] = ky[i];
1760  if(kx_lineNL[ic]<0){kx_lineNL[ic]+=nxNL;}
1761  if(ky_lineNL[ic]<0){ky_lineNL[ic]+=nyNL;}
1762  }//endfor
1763  }//endfor
1764  ic++;
1765 
1766  if(ic!=nline_tot){
1767  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1768  CkPrintf("Incorrect number of lines in util.C\n");
1769  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1770  CkExit();
1771  }//endif
1772 
1773 //////////////////////////////////////////////////////////////////////////////
1774 /// Create the line decomposition and a sorted run descriptor
1775 /// There are two rundescriptors per line : Noah's arc sort
1776 
1777  int nlines_max=0;
1778  for(int i=0;i<nchareG;i++){nlines_max=MAX(nlines_max,nline_lgrp[i]);}
1779  int **index_tran_upack = cmall_int_mat(0,nchareG,0,nlines_max,"util.C");
1780  int **index_tran_upackNL = cmall_int_mat(0,nchareG,0,nlines_max,"util.C");
1781 
1782  int yspace = sizeX;
1783  int yspaceNL = sizeXNL;
1784  if(doublePack){yspace =sizeX/2+1;}
1785  if(doublePack){yspaceNL=sizeXNL/2+1;}
1786  for(int igrp=0;igrp<nchareG;igrp++){
1787  for(int i=istrt_lgrp[igrp],j=0;i<iend_lgrp[igrp];i++,j++){
1788  index_tran_upack[igrp][j] = kx_line[i] + ky_line[i]*yspace;
1789  index_tran_upackNL[igrp][j] = kx_lineNL[i] + ky_lineNL[i]*yspaceNL;
1790  }//endfor
1791  }//endfor
1792 
1793  CkVec<RunDescriptor> *sortedRunDescriptors;
1794  sortedRunDescriptors = new CkVec<RunDescriptor> [nchareG];
1795  for(int igrp = 0; igrp < nchareG; igrp++){
1796  for(int i=istrt_lgrp[igrp];i<iend_lgrp[igrp];i++){
1797  int j = 2*i;
1798  int j1 = 2*i+1;
1799  sortedRunDescriptors[igrp].push_back(runDescriptorVec[j]);
1800  sortedRunDescriptors[igrp].push_back(runDescriptorVec[j1]);
1801  }//endfor
1802  }//endfor
1803 
1804  int *index_output_off = new int[nchareG];
1805  int numPoints = 0;
1806  for(int x = 0; x < nchareG; x ++) {
1807  index_output_off[x] = numPoints;
1808  int nnn = 0;
1809  int runsToBeSent = sortedRunDescriptors[x].size();
1810  for (int j = 0; j < runsToBeSent; j++){
1811  numPoints += sortedRunDescriptors[x][j].length;
1812  nnn += sortedRunDescriptors[x][j].length;
1813  }//endfor
1814  if(nnn != npts_lgrp[x]){
1815  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1816  CkPrintf("Incorrect number of points in gspace chare\n");
1817  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1818  CkExit();
1819  }//endif
1820  }//endfor
1821 
1822  if(numPoints!=numData){
1823  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1824  CkPrintf("Incorrect number of total g-space points\n");
1825  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1826  CkExit();
1827  }//endif
1828 
1829 //////////////////////////////////////////////////////////////////////////////
1830 /// Figure out the unique/redundent communication :
1831 /// Send out the unique : Receive and over write the redudant.
1832 
1833  // find the unique and redundent g's on each chare
1834  int *num_uni = new int [nchareG];
1835  int *num_red = new int [nchareG];
1836  int nk0_max=0;
1837  for(int i = 0; i < nchareG; i ++) {
1838  num_uni[i]=0;
1839  num_red[i]=0;
1840  if(doublePack==1){
1841  for(int j=0;j<npts_lgrp[i];j++){
1842  int iii = j+index_output_off[i];
1843  if(kx[iii]==0 && ky[iii]>0){num_uni[i]++;}
1844  if(kx[iii]==0 && ky[iii]<0){num_red[i]++;}
1845  if(kx[iii]==0 && ky[iii]==0 && kz[iii]>=0){num_uni[i]++;}
1846  if(kx[iii]==0 && ky[iii]==0 && kz[iii]<0){num_red[i]++;}
1847  }//endfor
1848  nk0_max=MAX(num_uni[i],nk0_max);
1849  nk0_max=MAX(num_red[i],nk0_max);
1850  }//endif
1851  }//endif
1852 
1853  // Find where my unique guys go and make a list so I can send them.
1854  // Make a list of where the unique guys arrive so I can receive them.
1855 
1856  RedundantCommPkg *RCommPkg = new RedundantCommPkg [nchareG];
1857  for(int i=0;i<nchareG;i++){RCommPkg[i].Init(nk0_max,nchareG);} // mallocs and zeros
1858 
1859  if(doublePack==1){
1860 
1861  for(int i=0;i<nchareG;i++){
1862  int *num_send = RCommPkg[i].num_send;
1863  int **lst_send = RCommPkg[i].lst_send;
1864  for(int j=0;j<nchareG;j++){
1865  int *num_recv = RCommPkg[j].num_recv;
1866  int **lst_recv = RCommPkg[j].lst_recv;
1867  for(int ip=0;ip<num_uni[i];ip++){
1868  int iii = ip+index_output_off[i] + num_red[i];
1869  for(int jp=0;jp<num_red[j];jp++){
1870  int jjj = jp+index_output_off[j];
1871  if(ky[iii]==-ky[jjj] && kz[iii]==-kz[jjj]){
1872  lst_send[j][num_send[j]] = ip+num_red[i];
1873  lst_recv[i][num_recv[i]] = jp;
1874  num_send[j]++;
1875  num_recv[i]++;
1876  }//endif
1877  }//endfor
1878  }//endfor
1879  }//endfor
1880  }//endfor
1881 
1882  for(int i=0;i<nchareG;i++){
1883  int *num_send = RCommPkg[i].num_send;
1884  int *num_recv = RCommPkg[i].num_recv;
1885  int num_recv_tot = 0;
1886  int num_send_tot = 0;
1887  for(int j=0;j<nchareG;j++){
1888  if(num_send[j]>0){num_send_tot++;}
1889  if(num_recv[j]>0){num_recv_tot++;}
1890  }//endif
1891  RCommPkg[i].num_recv_tot = num_recv_tot;
1892  RCommPkg[i].num_send_tot = num_send_tot;
1893  }//endfor
1894 
1895  for(int i=0;i<nchareG;i++){
1896  int *num_send = RCommPkg[i].num_send;
1897  int **lst_send = RCommPkg[i].lst_send;
1898  for(int j=0;j<nchareG;j++){
1899  int *num_recv = RCommPkg[j].num_recv;
1900  int **lst_recv = RCommPkg[j].lst_recv;
1901  if(num_send[j]!=num_recv[i]){
1902  CkPrintf("Bad Num of Redundent g-vectors %d %d\n",i,j);
1903  CkExit();
1904  }//endif
1905  for(int ip=0;ip<num_send[j];ip++){
1906  int iii = index_output_off[i]+lst_send[j][ip];
1907  int jjj = index_output_off[j]+lst_recv[i][ip];
1908  if(kx[iii]!=0 || kx[jjj]!=0 || ky[iii]!=-ky[jjj] || kz[iii]!=-kz[jjj] ||
1909  ky[iii]<0 || (kz[iii]<0 && ky[iii]==0)){
1910  CkPrintf("Bad Num of Redundent g-vectors %d %d %d %d %d\n",iii,jjj,ip,i,j);
1911  CkPrintf("%d %d %d : %d %d %d\n",kx[iii],ky[iii],kz[iii],
1912  kx[jjj],ky[jjj],kz[jjj]);
1913  CkExit();
1914  }//endif
1915  }//endfor
1916  }//endfor
1917  }//endfor
1918 
1919  }//endif :: Do redundancy when doublepack is off
1920 
1921 //////////////////////////////////////////////////////////////////////////////
1922 /// Pack up the stuff, clean up the memory and exit
1923 
1924  sim->npts_per_chareG = npts_lgrp;
1925  sim->index_tran_upack = index_tran_upack;
1926  sim->index_tran_upackNL = index_tran_upackNL;
1927  sim->nlines_max = nlines_max;
1928  sim->nlines_per_chareG = nline_lgrp;
1929  sim->sortedRunDescriptors = sortedRunDescriptors;
1930  sim->npts_tot = numData;
1931  sim->nlines_tot = nline_tot;
1932  sim->index_output_off = index_output_off;
1933  sim->RCommPkg = RCommPkg;
1934 
1935  delete [] istrt_lgrp;
1936  delete [] iend_lgrp;
1937  delete [] complexPoints;
1938  delete [] kx_line;
1939  delete [] ky_line;
1940  delete [] kx_lineNL;
1941  delete [] ky_lineNL;
1942  delete [] kx;
1943  delete [] ky;
1944  delete [] kz;
1945  delete [] num_uni;
1946  delete [] num_red;
1947 
1948 //////////////////////////////////////////////////////////////////////////////
1949  }//end routine
1950 //////////////////////////////////////////////////////////////////////////////
1951 
1952 
1953 ///////////////////////////////////////////////////////////////////////////////
1954 ///////////////////////////////////////////////////////////////////////////////
1955 ///////////////////////////////////////////////////////////////////////////////
1956 void writeStateFile(int ncoef,complex *psi,complex *vpsi,
1957  int *k_x,int *k_y,int *k_z,int cp_min_opt,
1958  int sizeX,int sizeY,int sizeZ,char *psiName,char *vpsiName,
1959  int ibinary_write_opt,int iteration, int istate,
1960  int ispin, int ikpt, int ibead, int itemper)
1961 ///////////////////////////////////////////////////////////////////////////////
1962  { //begin rotunie
1963 ///////////////////////////////////////////////////////////////////////////////
1964 
1965  int *index = new int [ncoef];
1966  int *ktemp = new int [ncoef];
1967  int istrt=0;
1968  sort_psi_output(ncoef,k_x,k_y,k_z,index,ktemp,&istrt);
1969  int ncoef_true=ncoef-istrt;
1970 
1971  if(istate==1){
1972  char fname[1000];
1973  sprintf (fname, "%s/Spin.%d_Kpt.%d_Bead.%d_Temper.%d/TimeStamp",config.dataPathOut,ispin,ikpt,ibead,itemper);
1974  FILE *fp = fopen(fname,"w");
1975  fprintf(fp,"time step = %d\n",iteration);
1976  fclose(fp);
1977  }//endif
1978 
1979  if(ibinary_write_opt==0){
1980  FILE *fp = fopen(psiName,"w");
1981  assert(fp!=NULL);
1982  fprintf(fp,"%d %d %d %d\n",ncoef_true,sizeX,sizeY,sizeZ);
1983  for(int i=istrt+1;i<ncoef;i++){
1984  fprintf(fp,"%g %g %d %d %d \n",psi[index[i]].re,psi[index[i]].im,
1985  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
1986  }//endfor
1987  int i = istrt;
1988  fprintf(fp,"%g %g %d %d %d \n",psi[index[i]].re,psi[index[i]].im,
1989  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
1990  fclose(fp);
1991  }
1992 #if CMK_PROJECTIONS_USE_ZLIB
1993  else if(ibinary_write_opt==2){
1994  strcat(psiName,".gz");
1995  gzFile zfp = gzopen(psiName,"w");
1996  assert(zfp!=NULL);
1997  gzprintf(zfp,"%d %d %d %d\n",ncoef_true,sizeX,sizeY,sizeZ);
1998  for(int i=istrt+1;i<ncoef;i++){
1999  gzprintf(zfp,"%g %g %d %d %d \n",psi[index[i]].re,psi[index[i]].im,
2000  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
2001  }//endfor
2002  int i = istrt;
2003  gzprintf(zfp,"%g %g %d %d %d \n",psi[index[i]].re,psi[index[i]].im,
2004  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
2005  gzclose(zfp);
2006  }
2007  else if(ibinary_write_opt==3){
2008  strcat(psiName,".gz");
2009  gzFile zfp = gzopen(psiName,"w");
2010  assert(zfp!=NULL);
2011  int n=1;
2012  gzwrite(zfp,&ncoef_true,sizeof(int));
2013  gzwrite(zfp,&sizeX,sizeof(int));
2014  gzwrite(zfp,&sizeY,sizeof(int));
2015  gzwrite(zfp,&sizeZ,sizeof(int));
2016  for(int i=istrt+1;i<ncoef;i++){
2017  gzwrite(zfp,&psi[index[i]].re,sizeof(double));
2018  gzwrite(zfp,&psi[index[i]].im,sizeof(double));
2019  gzwrite(zfp,&k_x[index[i]],sizeof(int));
2020  gzwrite(zfp,&k_y[index[i]],sizeof(int));
2021  gzwrite(zfp,&k_z[index[i]],sizeof(int));
2022  }//endfor
2023  int i = istrt;
2024  gzwrite(zfp,&psi[index[i]].re,sizeof(double));
2025  gzwrite(zfp,&psi[index[i]].im,sizeof(double));
2026  gzwrite(zfp,&k_x[index[i]],sizeof(int));
2027  gzwrite(zfp,&k_y[index[i]],sizeof(int));
2028  gzwrite(zfp,&k_z[index[i]],sizeof(int));
2029  gzclose(zfp);
2030  }
2031 #endif
2032  else{
2033  FILE *fp = fopen(psiName,"w");
2034  assert(fp!=NULL);
2035  int n=1;
2036  fwrite(&ncoef_true,sizeof(int),n,fp);
2037  fwrite(&sizeX,sizeof(int),n,fp);
2038  fwrite(&sizeY,sizeof(int),n,fp);
2039  fwrite(&sizeZ,sizeof(int),n,fp);
2040  for(int i=istrt+1;i<ncoef;i++){
2041  fwrite(&psi[index[i]].re,sizeof(double),n,fp);
2042  fwrite(&psi[index[i]].im,sizeof(double),n,fp);
2043  fwrite(&k_x[index[i]],sizeof(int),n,fp);
2044  fwrite(&k_y[index[i]],sizeof(int),n,fp);
2045  fwrite(&k_z[index[i]],sizeof(int),n,fp);
2046  }//endfor
2047  int i = istrt;
2048  fwrite(&psi[index[i]].re,sizeof(double),n,fp);
2049  fwrite(&psi[index[i]].im,sizeof(double),n,fp);
2050  fwrite(&k_x[index[i]],sizeof(int),n,fp);
2051  fwrite(&k_y[index[i]],sizeof(int),n,fp);
2052  fwrite(&k_z[index[i]],sizeof(int),n,fp);
2053  fclose(fp);
2054  }//endif
2055 
2056 
2057  if(cp_min_opt==0){
2058 
2059  if(ibinary_write_opt==0){
2060  FILE *fp = fopen(vpsiName,"w");
2061  assert(fp!=NULL);
2062  fprintf(fp,"%d %d %d %d\n",ncoef_true,sizeX,sizeY,sizeZ);
2063  for(int i=istrt+1;i<ncoef;i++){
2064  fprintf(fp,"%g %g %d %d %d \n",vpsi[index[i]].re,vpsi[index[i]].im,
2065  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
2066  }//endfor
2067  int i = istrt;
2068  fprintf(fp,"%g %g %d %d %d \n",vpsi[index[i]].re,vpsi[index[i]].im,
2069  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
2070  fclose(fp);
2071  }
2072 #if CMK_PROJECTIONS_USE_ZLIB
2073  else if(ibinary_write_opt==2){
2074  strcat(vpsiName,".gz");
2075  gzFile zfp=gzopen(vpsiName,"w");
2076  assert(zfp!=NULL);
2077  gzprintf(zfp,"%d %d %d %d\n",ncoef_true,sizeX,sizeY,sizeZ);
2078  for(int i=istrt+1;i<ncoef;i++){
2079  gzprintf(zfp,"%g %g %d %d %d \n",vpsi[index[i]].re,vpsi[index[i]].im,
2080  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
2081  }//endfor
2082  int i = istrt;
2083  gzprintf(zfp,"%g %g %d %d %d \n",vpsi[index[i]].re,vpsi[index[i]].im,
2084  k_x[index[i]],k_y[index[i]],k_z[index[i]]);
2085  gzclose(zfp);
2086 
2087  }
2088  else if(ibinary_write_opt==3){
2089  strcat(vpsiName,".gz");
2090  gzFile zfp=gzopen(vpsiName,"w");
2091  assert(zfp!=NULL);
2092  int n=1;
2093  gzwrite(zfp,&ncoef_true,sizeof(int));
2094  gzwrite(zfp,&sizeX,sizeof(int));
2095  gzwrite(zfp,&sizeY,sizeof(int));
2096  gzwrite(zfp,&sizeZ,sizeof(int));
2097  for(int i=istrt+1;i<ncoef;i++){
2098  gzwrite(zfp,&vpsi[index[i]].re,sizeof(double));
2099  gzwrite(zfp,&vpsi[index[i]].im,sizeof(double));
2100  gzwrite(zfp,&k_x[index[i]],sizeof(int));
2101  gzwrite(zfp,&k_y[index[i]],sizeof(int));
2102  gzwrite(zfp,&k_z[index[i]],sizeof(int));
2103  }//endfor
2104  int i = istrt;
2105  gzwrite(zfp,&vpsi[index[i]].re,sizeof(double));
2106  gzwrite(zfp,&vpsi[index[i]].im,sizeof(double));
2107  gzwrite(zfp,&k_x[index[i]],sizeof(int));
2108  gzwrite(zfp,&k_y[index[i]],sizeof(int));
2109  gzwrite(zfp,&k_z[index[i]],sizeof(int));
2110  gzclose(zfp);
2111  }//endif
2112 #endif
2113 else{
2114  FILE *fp = fopen(vpsiName,"w");
2115  assert(fp!=NULL);
2116  int n=1;
2117  fwrite(&ncoef_true,sizeof(int),n,fp);
2118  fwrite(&sizeX,sizeof(int),n,fp);
2119  fwrite(&sizeY,sizeof(int),n,fp);
2120  fwrite(&sizeZ,sizeof(int),n,fp);
2121  for(int i=istrt+1;i<ncoef;i++){
2122  fwrite(&vpsi[index[i]].re,sizeof(double),n,fp);
2123  fwrite(&vpsi[index[i]].im,sizeof(double),n,fp);
2124  fwrite(&k_x[index[i]],sizeof(int),n,fp);
2125  fwrite(&k_y[index[i]],sizeof(int),n,fp);
2126  fwrite(&k_z[index[i]],sizeof(int),n,fp);
2127  }//endfor
2128  int i = istrt;
2129  fwrite(&vpsi[index[i]].re,sizeof(double),n,fp);
2130  fwrite(&vpsi[index[i]].im,sizeof(double),n,fp);
2131  fwrite(&k_x[index[i]],sizeof(int),n,fp);
2132  fwrite(&k_y[index[i]],sizeof(int),n,fp);
2133  fwrite(&k_z[index[i]],sizeof(int),n,fp);
2134  fclose(fp);
2135  }//endif
2136  }//endif
2137 
2138  delete [] index;
2139  delete [] ktemp;
2140 
2141 //////////////////////////////////////////////////////////////////////////////
2142  }//end routine
2143 //////////////////////////////////////////////////////////////////////////////
2144 
2145 
2146 ///////////////////////////////////////////////////////////////////////////////
2147 ///////////////////////////////////////////////////////////////////////////////
2148 ///////////////////////////////////////////////////////////////////////////////
2149 void sort_psi_output(int n,int *kx,int *ky,int *kz,int *index,int *ktemp,int *istrt_ret){
2150 ///////////////////////////////////////////////////////////////////////////////
2151 /// Find max and min
2152 
2153  int kxmax = kx[0]; int kxmin=kx[0];
2154  int kymax = ky[0]; int kymin=ky[0];
2155  int kzmax = kz[0]; int kzmin=kz[0];
2156  for(int i=0;i<n;i++){
2157  kxmax = MAX(kxmax,kx[i]); kxmin = MIN(kxmin,kx[i]);
2158  kymax = MAX(kymax,ky[i]); kymin = MIN(kymin,ky[i]);
2159  kzmax = MAX(kzmax,kz[i]); kzmin = MIN(kzmin,kz[i]);
2160  }//endfor
2161  int nx = kxmax-kxmin+1;
2162  int ny = kymax-kymin+1;
2163  int nz = kzmax-kzmin+1;
2164 
2165 ///////////////////////////////////////////////////////////////////////////////
2166 /// Create a sortable 1d array and sort
2167 
2168 
2169  for(int i=0;i<n;i++){index[i]=i;}
2170  for(int i=0;i<n;i++){ktemp[i]=(kz[i]-kzmin)+(ky[i]-kymin)*nz+(kx[i]-kxmin)*nz*ny;}
2171  sort_commence(n,ktemp,index);
2172 
2173 ///////////////////////////////////////////////////////////////////////////////
2174 /// Find g=0 term
2175 
2176  int istrt=0;
2177  for(int i=0;i<n;i++){
2178  if(kx[index[i]]==0&&ky[index[i]]==0&&kz[index[i]]==0){
2179  istrt=i; break;
2180  }//endif
2181  }//endfor
2182  (*istrt_ret) = istrt;
2183 
2184 //-----------------------------------------------------------------------------
2185  }// end routine
2186 ///////////////////////////////////////////////////////////////////////////////
2187 
2188 
2189 ///////////////////////////////////////////////////////////////////////////////
2190 ///////////////////////////////////////////////////////////////////////////////
2191 ///////////////////////////////////////////////////////////////////////////////
2192 void sort_kxky(int n,int *kx,int *ky,int *index,int *ktemp,int sizeY){
2193 ///////////////////////////////////////////////////////////////////////////////
2194 /// Sort on kx
2195 
2196  for(int i=0;i<n;i++){index[i]=i;}
2197  sort_commence(n,kx,index);
2198  for(int i=0;i<n;i++){ktemp[i]=ky[i];}
2199  for(int i=0;i<n;i++){ky[i]=ktemp[index[i]];}
2200 
2201 ///////////////////////////////////////////////////////////////////////////////
2202 /// Sort on ky for each kx
2203 
2204  int kmin = kx[0];
2205  int kmax = kx[(n-1)];
2206  for(int i=0;i<=kmax-kmin;i++){ktemp[i]=0;}
2207  for(int i=0;i<n;i++){ktemp[kx[i]-kmin]++;}
2208 
2209  int ioff=0;
2210  for(int i=0;i<=kmax-kmin;i++){
2211  if(ktemp[i]>1){sort_commence(ktemp[i],&ky[ioff],&index[ioff]);}
2212  ioff += ktemp[i];
2213  }//endfor
2214 
2215 //////////////////////////////////////////////////////////////////////////////
2216  }//end routine
2217 //////////////////////////////////////////////////////////////////////////////
2218 
2219 ///////////////////////////////////////////////////////////////////////////////
2220 ///////////////////////////////////////////////////////////////////////////////
2221 ///////////////////////////////////////////////////////////////////////////////
2222 void sort_kxky_old(int n,int *kx,int *ky,int *index,int *kyt){
2223 ///////////////////////////////////////////////////////////////////////////////
2224 
2225  int nk0=0;
2226  for(int i=0;i<n;i++){index[i]=i;}
2227  for(int i=0;i<n;i++){
2228  if(kx[i]==0){
2229  int itmp = index[nk0];
2230  index[nk0] = index[i];
2231  index[i] = itmp;
2232  nk0++;
2233  }//endif
2234  }//endfor
2235 
2236  for(int i=0;i<nk0;i++){kyt[i]=ky[index[i]];}
2237  sort_commence(nk0,kyt,index);
2238 
2239 //////////////////////////////////////////////////////////////////////////////
2240  }//end routine
2241 //////////////////////////////////////////////////////////////////////////////
2242 
2243 
2244 
2245 //////////////////////////////////////////////////////////////////////////////
2246 /// Initialization Function
2247 //////////////////////////////////////////////////////////////////////////////
2248 //////////////////////////////////////////////////////////////////////////////
2249 //////////////////////////////////////////////////////////////////////////////
2250 void getSplitDecomp(int *istrt_ret,int *iend_ret,int *n_ret,
2251  int ntot, int ndiv,int idiv)
2252 //////////////////////////////////////////////////////////////////////////////
2253  {//begin routine
2254 //////////////////////////////////////////////////////////////////////////////
2255 
2256  if(idiv>=ndiv || ntot< ndiv){
2257  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
2258  CkPrintf("Incorrect input to RhoGHart collection creator.\n");
2259  CkPrintf("idiv %d ndiv %d, ntot %d.\n");
2260  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
2261  CkExit();
2262  }//endif
2263 
2264  int n = (ntot/ndiv);
2265  int r = (ntot%ndiv);
2266 
2267  int istrt = n*idiv;
2268  if(idiv>=r){istrt += r;}
2269  if(idiv<r) {istrt += idiv;}
2270  if(idiv<r) {n++;}
2271  int iend = n+istrt;
2272 
2273  if(n==0){
2274  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
2275  CkPrintf("No lines in a RhoGHart collection!!\n");
2276  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
2277  CkExit();
2278  }//endif
2279 
2280  (*n_ret) = n;
2281  (*istrt_ret) = istrt;
2282  (*iend_ret) = iend;
2283 
2284 //---------------------------------------------------------------------------
2285  }//end routine
2286 //////////////////////////////////////////////////////////////////////////////
2287 
2288 
2289 
2290 //////////////////////////////////////////////////////////////////////////////
2291 /// Create some decompositions and find the best one.
2292 //////////////////////////////////////////////////////////////////////////////
2293 //////////////////////////////////////////////////////////////////////////////
2294 //////////////////////////////////////////////////////////////////////////////
2295 void create_subPlane_decomp(int nplane_x,int *listGx,int *mapGrpGx,
2296  int nchareRhoGEext,int *numSubGx,
2297  int *nline_lgrp_eext,int *kx_line,
2298  int **nline_send_eext_y, int rhoRsubplanes){
2299 //////////////////////////////////////////////////////////////////////////////
2300 /// try out some several decomps and find the best
2301 
2302  int *list1 = new int [nplane_x];
2303  long seed = 1145;
2304 
2305  int score_max = 0;
2306  for(int j=0;j<nplane_x;j++){list1[j]=j;}
2307  for(int ntry=1;ntry<100;ntry++){
2308 
2309  //------------------------------------------------------
2310  // Mix up the kx
2311  if(ntry>1){
2312  double stuff = altRandom(&seed);
2313  int index = (int) (((double)nplane_x)*stuff);
2314  index = MIN(index,nplane_x-1);
2315  stuff = altRandom(&seed);
2316  int jndex = (int) (((double)nplane_x)*stuff);
2317  jndex = MIN(jndex,nplane_x-1);
2318  int itemp = list1[index];
2319  list1[index] = list1[jndex];
2320  list1[jndex] = itemp;
2321  }//endif
2322 
2323  //------------------------------------------------------
2324  // Create the decomp
2325  int iii = 0;
2326  for(int igrp=0;igrp<rhoRsubplanes;igrp++){
2327  int num = numSubGx[igrp];
2328  for(int jc=0,ic=iii;ic<iii+num;ic++,jc++){
2329  mapGrpGx[list1[ic]] = igrp;
2330  }//endfor
2331  iii += num;
2332  }//endfor
2333 
2334  //------------------------------------------------------
2335  // Score the decomp
2336  int score;
2337  score_subPlane_decomp(nchareRhoGEext,rhoRsubplanes,nline_lgrp_eext,
2338  mapGrpGx,kx_line,nline_send_eext_y,&score);
2339 
2340  //------------------------------------------------------
2341  // Keep the best decomp
2342  // CkPrintf("try %d score %d score_max %d : ",ntry,score,score_max);
2343  // for(int i=0;i<nplane_x;i++){CkPrintf("%d ",list1[i]);} CkPrintf("\n");
2344  if(score<score_max || ntry==1){
2345  score_max = score;
2346  for(int i=0;i<nplane_x;i++){listGx[i]=list1[i];}
2347  }else{
2348  for(int i=0;i<nplane_x;i++){list1[i]=listGx[i];}
2349  }//endif
2350 
2351  }//endfor
2352 
2353  delete [] list1;
2354 
2355 //////////////////////////////////////////////////////////////////////////////
2356  }//end routine
2357 //////////////////////////////////////////////////////////////////////////////
2358 
2359 
2360 //////////////////////////////////////////////////////////////////////////////
2361 /// Score a decomposition
2362 //////////////////////////////////////////////////////////////////////////////
2363 //////////////////////////////////////////////////////////////////////////////
2364 //////////////////////////////////////////////////////////////////////////////
2365 void score_subPlane_decomp(int nchareRhoGEext,int rhoRsubplanes, int *nline_lgrp_eext,
2366  int *mapGrpGx, int *kx_line,int **nline_send_eext_y,
2367  int *Rhart_max){
2368 //////////////////////////////////////////////////////////////////////////////
2369 /// Count the sends
2370 
2371  for(int igrp=0,i=0;igrp<nchareRhoGEext;igrp++){
2372  for(int ic=0;ic<rhoRsubplanes;ic++){
2373  nline_send_eext_y[igrp][ic]=0;
2374  }}//endfor
2375 
2376  for(int igrp=0,i=0;igrp<nchareRhoGEext;igrp++){
2377  for(int j=0;j<nline_lgrp_eext[igrp];j++,i++){
2378  int ic = mapGrpGx[kx_line[i]]; //subPlane index
2379  nline_send_eext_y[igrp][ic]++;
2380  }//endfor
2381  }//endfor
2382 
2383 //////////////////////////////////////////////////////////////////////////////
2384 /// Find the total number of send/recvs
2385 
2386  Rhart_max[0] = 0;
2387  for( int j=0; j< nchareRhoGEext;j++){
2388  int recvCountFromRHartExt = 0;
2389  for(int i=0;i<rhoRsubplanes;i++){
2390  if(nline_send_eext_y[j][i]>0){recvCountFromRHartExt++;}
2391  }//endfor
2392  Rhart_max[0]+=recvCountFromRHartExt;
2393  }//endfor
2394 
2395 //////////////////////////////////////////////////////////////////////////////
2396  }//end routine
2397 //////////////////////////////////////////////////////////////////////////////
2398 
2399 
2400 //////////////////////////////////////////////////////////////////////////////
2401 //////////////////////////////////////////////////////////////////////////////
2402 //////////////////////////////////////////////////////////////////////////////
2403 void create_gx_decomp(int nktot, int nline, int *kx_line, int *mapl,
2404  int nchare,int *nline_grp){
2405 //////////////////////////////////////////////////////////////////////////////
2406 /// count the pts in the lines
2407 
2408  int *npts = new int [nline];
2409  for(int i =0;i<nline;i++){npts[i]=0;}
2410  for(int i =0;i<nktot;i++){npts[mapl[kx_line[i]]]++;} // ky pts per kx line
2411 
2412 //////////////////////////////////////////////////////////////////////////////
2413 /// Find the best division
2414 
2415  double dev_min = 0.0;
2416  int nbal_min = 0;
2417  int ibal_min = 0;
2418  int ifirst = 0;
2419  int mmm = (nktot/nchare);
2420  for(int ibal=0;ibal<=mmm-1;ibal++){
2421  int nbal = ibal;
2422  int ntarg = (nktot/nchare);
2423  if(ntarg > nbal){ntarg -= nbal;}
2424  int nmax = 0;
2425  int nmin = nktot;
2426  int nnow = 0;
2427  int ic = 0;
2428  for(int i=0;i<nline;i++){
2429  nnow += npts[i];
2430  if( (nnow>=ntarg) && (ic<(nchare-1)) ){
2431  ic+=1;
2432  nmin = MIN(nmin,nnow);
2433  nmax = MAX(nmax,nnow);
2434  nnow = 0;
2435  }//endif
2436  }//endfor
2437  nmin = MIN(nmin,nnow);
2438  nmax = MAX(nmax,nnow);
2439  double dev = 100.0*((double)(nmax-nmin))/((double)MAX(nmin,1));
2440  if(ic==nchare-1){
2441  if(dev<dev_min || ifirst==0){
2442  ifirst = 1;
2443  dev_min = dev;
2444  nbal_min = nbal;
2445  ibal_min = ibal;
2446  }//endif
2447  }//endif
2448  }//endfor
2449 
2450 ///////////////////////////////////////////////////////////////////////////=
2451 /// Store the good decomposition
2452 
2453  int ntarg = (nktot/nchare);
2454  if(ntarg > nbal_min){ntarg = ntarg-nbal_min;}
2455 
2456  int ic = 0;
2457  int nnow = 0;
2458  for(int i=0;i<nchare;i++){nline_grp[i]=0;}
2459  for(int i=0;i<nline;i++){
2460  nline_grp[ic]++;
2461  nnow += npts[i];
2462  if( (nnow>=ntarg) && (ic<(nchare-1)) ){
2463  ic +=1;
2464  nnow = 0;
2465  }//endif
2466  }//endfor
2467 
2468  delete [] npts;
2469 
2470 //////////////////////////////////////////////////////////////////////////////
2471  }//end routine
2472 //////////////////////////////////////////////////////////////////////////////
2473 
2474 
2475 
2476 //////////////////////////////////////////////////////////////////////////////
2477 //////////////////////////////////////////////////////////////////////////////
2478 //////////////////////////////////////////////////////////////////////////////
2479 FILE *openScreenfWrite(const char *dirnameBase, const char *fname, int temper, int bead, bool beadfile)
2480 {
2481  char dirPath[1024];
2482  char subdir[1024];
2483  char lfname[1024];
2484  char beadname[1024];
2485  memset(dirPath, 0 , 1024);
2486  memset(lfname, 0 , 1024);
2487  snprintf(subdir,1023,"/Temper.%d/",temper);
2488  strncpy(dirPath, dirnameBase, 1023);
2489  strncpy(lfname, fname, 1023);
2490  strncat(dirPath,subdir,1023);
2491  if(beadfile){
2492  snprintf(beadname,1023,".Bead.%d",bead);
2493  strncat(lfname,beadname,1023);
2494  }
2495  strncat(dirPath,lfname,1023);
2496  FILE *file=fopen(dirPath,"a");
2497  // CkPrintf("opened temper file %s for screen output \n",dirPath);
2498  assert(file!=NULL);
2499  return file;
2500 //////////////////////////////////////////////////////////////////////////////
2501 }//end routine
2502 //////////////////////////////////////////////////////////////////////////////
openScreenfWrite
FILE * openScreenfWrite(const char *dirnameBase, const char *fname, int temper, int bead, bool beadfile)
this stuff is in a C file because multiple includes of .def will give you too many implementations of...
Definition: util.C:2479
make_rho_runs
void make_rho_runs(CPcharmParaInfo *sim)
== rho space run descriptors for spherical cutoff fft support
Definition: util.C:34
get_rho_kvectors
void get_rho_kvectors(double ecut4, double *hmati, int **kx_ret, int **ky_ret, int **kz_ret, int *nline_tot_ret, int *nPacked_ret, int ka_max, int kb_max, int kc_max)
Definition: util.C:777
sort_psi_output
void sort_psi_output(int n, int *kx, int *ky, int *kz, int *index, int *ktemp, int *istrt_ret)
Definition: util.C:2149
sort_kxky
void sort_kxky(int n, int *kx, int *ky, int *index, int *ktemp, int sizeY)
Definition: util.C:2192
readStateIntoRuns
void readStateIntoRuns(int nPacked, int ncoef, complex *arrCP, CkVec< RunDescriptor > &runs, const char *fromFile, int ibinary_opt, int *nline_tot_ret, int *nplane_ret, int *istrt_lgrp, int *iend_lgrp, int *npts_lgrp, int *nline_lgrp, int **kx_line_ret, int **ky_line_ret, int **kx_ret, int **ky_ret, int **kz_ret, int iget_decomp, int gen_wave, int nx_in, int ny_in, int nz_in)
Definition: util.C:862
RunDescriptor
== Index logic for lines of constant x,y in gspace.
Definition: RunDescriptor.h:22
config
Config config
addtogroup Uber
Definition: cpaimd.ci:233
ParaGrpParse::get_chareG_line_prms
static void get_chareG_line_prms(int, int, int, int *, int *, int *, int *, int *, bool)
functions
Definition: para_grp_parse.C:21
readState
void readState(int nPacked, complex *arrCP, const char *fromFile, int ibinary_opt, int *nline_tot_ret, int *nplane_ret, int *kx, int *ky, int *kz, int *nx_ret, int *ny_ret, int *nz_ret, int *istrt_lgrp, int *iend_lgrp, int *npts_lgrp, int *nline_lgrp, int iget_decomp, int iget_vstate)
Definition: util.C:1108
getSplitDecomp
void getSplitDecomp(int *istrt_ret, int *iend_ret, int *n_ret, int ntot, int ndiv, int idiv)
Initialization Function /////////////////////////////////////////////////////////////////////////// /...
Definition: util.C:2250
ParaGrpParse::flip_data_set
static void flip_data_set(int, int *, int *, int *, int *)
Definition: LINE_EQUAL_DECOMP_V.2/para_grp_parse.C:135
create_gx_decomp
void create_gx_decomp(int nktot, int nline, int *kx_line, int *mapl, int nchare, int *nline_grp)
Definition: util.C:2403
score_subPlane_decomp
void score_subPlane_decomp(int nchareRhoGEext, int rhoRsubplanes, int *nline_lgrp_eext, int *mapGrpGx, int *kx_line, int **nline_send_eext_y, int *Rhart_max)
Score a decomposition /////////////////////////////////////////////////////////////////////////// ///...
Definition: util.C:2365
create_subPlane_decomp
void create_subPlane_decomp(int nplane_x, int *listGx, int *mapGrpGx, int nchareRhoGEext, int *numSubGx, int *nline_lgrp_eext, int *kx_line, int **nline_send_eext_y, int rhoRsubplanes)
Create some decompositions and find the best one. ///////////////////////////////////////////////////...
Definition: util.C:2295
create_line_decomp_descriptor
void create_line_decomp_descriptor(CPcharmParaInfo *sim)
/////////////////////////////////////////////////////////////////////////cc
Definition: util.C:1651
processState
void processState(int nPacked, int nktot, complex *arrCP, const char *fromFile, int ibinary_opt, int *nline_tot_ret, int *nplane_ret, int *kx, int *ky, int *kz, int *istrt_lgrp, int *iend_lgrp, int *npts_lgrp, int *nline_lgrp, int iget_decomp, int iget_vstate, int iopt, int ny)
Definition: util.C:1332