4d595fe1657a8aefa33cc770aacd72c28a0a897b
[charm.git] / examples / charm++ / barnes-charm / Main.cpp
1 #include "barnes.h"
2 #include "defs.h"
3
4 CProxy_Main mainChare;
5 CProxy_TreePiece pieces;
6 CProxy_ParticleChunk chunks;
7 int maxmybody;
8
9 int numTreePieces;
10 int maxPartsPerTp;
11 int numParticleChunks;
12
13 // number of particle chunks 
14 int NPROC;
15 real fcells;
16 real fleaves;
17 real tstop;
18 int nbody;
19 real dtime;
20 real eps;
21 real epssq;
22 real tol;
23 real tolsq;
24 real dtout;
25 real dthf;
26
27 int maxmycell;
28 int maxmyleaf;
29
30 CkReduction::reducerType minmax_RealVectorType;
31
32 int log8floor(int arg){
33   int ret = -1;
34   for(; arg > 0; arg >>= 3){
35     ret++;
36   }
37   return (ret);
38 }
39
40 /*
41  * INITOUTPUT: initialize output routines.
42  */
43
44
45 void Main::initoutput()
46 {
47    printf("\n\t\t%s\n\n", headline.c_str());
48    printf("%10s%10s%10s%10s%10s%10s%10s%10s%10s\n",
49           "nbody", "dtime", "eps", "tol", "dtout", "tstop","fcells","fleaves","NPROC");
50    printf("%10d%10.5f%10.4f%10.2f%10.3f%10.3f%10.2f%10.2f%10d\n\n",
51           nbody, dtime, eps, tol, dtout, tstop, fcells, fleaves, NPROC);
52 }
53
54
55 Main::Main(CkArgMsg *m){
56   int c;
57
58   mainChare = thisProxy;
59
60   defaults.reserve(32);
61   for(int i = 0; i < m->argc; i++){
62     defaults.push_back(m->argv[i]);
63   }
64   startrun();
65   initoutput();
66   maxleaf = (int) ((double) fleaves * nbody);
67   CkPrintf("[main] maxleaf: %d, fleaves: %f, nbody: %d\n", maxleaf, fleaves, nbody);
68   tab_init();
69
70   maxPartsPerTp = getiparam("fat"); 
71   if(maxPartsPerTp < 0){
72     maxPartsPerTp = MAX_PARTS_PER_TP;
73   }
74   else if(maxPartsPerTp < MAX_BODIES_PER_LEAF){
75     maxPartsPerTp = MAX_BODIES_PER_LEAF;
76   }
77
78   
79   numParticleChunks = NPROC;
80   iterations = getiparam("it");
81   if(iterations < 0){
82     iterations = DEFAULT_NUM_ITERATIONS;
83   }
84
85   /*
86   numTreePieces = getiparam("pieces");
87   if(numTreePieces < 0){
88     numTreePieces = 8*numParticleChunks; 
89   }
90   */
91   int depth = getiparam("depth");
92   if(depth < 1){
93     depth = 1;
94   }
95
96   numTreePieces = 1 << (3*depth);
97   ckout << "pieces: " << numTreePieces << endl;
98   // various maximum count parameters are:
99   // nbody, maxmybody, maxcell, maxmycell, maxleaf, maxmyleaf
100   // nbody has already been set
101   maxmyleaf = maxleaf/NPROC;
102   maxcell = fcells * maxleaf;
103   maxmycell = maxcell/NPROC;
104
105
106   // create chunks, treepieces
107   CProxy_BlockMap myMap=CProxy_BlockMap::ckNew(); 
108   CkArrayOptions opts(numTreePieces); 
109   opts.setMap(myMap);
110   //int depth = log8floor(numTreePieces);
111   ckout << "top-level pieces depth: " << (depth+1) << ", " << (IMAX >> (depth+1)) << endl;
112   CProxy_TreePiece treeProxy = CProxy_TreePiece::ckNew((CmiUInt8)0,-1,(IMAX >> (depth+1)), CkArrayIndex1D(0), opts);
113   pieces = treeProxy;
114
115   myMap=CProxy_BlockMap::ckNew(); 
116   CkArrayOptions optss(numParticleChunks); 
117   optss.setMap(myMap);
118   CProxy_ParticleChunk chunkProxy = CProxy_ParticleChunk::ckNew(maxleaf, maxcell, optss);
119   chunks = chunkProxy;
120
121   topLevelRoots.reserve(numTreePieces);
122   // startup split into two so that global readonlys
123   // are initialized before we send start signal to 
124   // particle chunks
125   ckout << "Starting simulation" << endl;
126   thisProxy.startSimulation();
127 }
128
129 /*
130  * MAKECELL: allocation routine for cells.
131  */
132
133 cellptr Main::makecell(unsigned ProcessId)
134 {
135    cellptr c;
136    int i, Mycell;
137     
138    c = new cell;
139    Type(c) = CELL;
140    Done(c) = FALSE;
141    Mass(c) = 0.0;
142    for (i = 0; i < NSUB; i++) {
143       Subp(c)[i] = NULL;
144    }
145    return (c);
146 }
147
148 /*
149  * INIT_ROOT: Processor 0 reinitialize the global root at each time step
150  */
151 void Main::init_root (unsigned int ProcessId)
152 {
153
154   // create top portion of global tree
155   int depth = log8floor(numTreePieces);
156   G_root = makecell(ProcessId);
157   ParentOf(G_root) = NULL;
158   Mass(G_root) = 0.0;
159   Cost(G_root) = 0;
160   CLRV(Pos(G_root));
161   NodeKey(G_root) = nodekey(1); 
162
163
164   Level(G_root) = IMAX >> 1;
165   
166   CkPrintf("[main] Creating top-level tree, depth: %d, root: 0x%x\n", depth, G_root);
167   int totalNumCellsMade = createTopLevelTree(G_root, depth);
168   CkPrintf("totalNumCellsMade: %d\n", totalNumCellsMade+1);
169 }
170
171 int Main::createTopLevelTree(cellptr node, int depth){
172   if(depth == 1){
173     // lowest level, no more children will be created - save self
174     int index = topLevelRoots.push_back_v((nodeptr) node);
175 #ifdef VERBOSE_MAIN
176     CkPrintf("saving root %d 0x%x\n", index, node);
177 #endif
178     return 0;
179   }
180   
181   int numCellsMade = 0;
182   // here, ancestors of treepieces' roots are created,
183   // not the roots themselves.
184   // these are all cells. the roots will start out as NULL pointers,
185   // but if bodies are sent to them, they are modified to become
186   // leaves.
187
188   for (int i = 0; i < NSUB; i++){
189     cellptr child = makecell(-1);
190     Subp(node)[i] = (nodeptr) child;
191     ParentOf(child) = (nodeptr) node;
192     ChildNum(child) = i;
193     Level(child) = Level(node) >> 1;
194     nodekey k = NodeKey(node) << NDIM;
195     k += i;
196     NodeKey(child) = k;
197
198     Mass(child) = 0.0;
199     Cost(child) = 0;
200     CLRV(Pos(child));
201
202     numCellsMade++;
203     numCellsMade += createTopLevelTree((cellptr) Subp(node)[i], depth-1);
204   }
205
206   return numCellsMade;
207 }
208
209 void Main::startSimulation(){
210   // slavestart for chunks
211   chunks.SlaveStart((CmiUInt8)mybodytab, (CmiUInt8)bodytab, CkCallbackResumeThread());
212
213   double start;
214   double end;
215   double totalStart;
216   double iterationStart;
217   double totalEnd;
218
219   /* main loop */
220   int i = 0;
221   while (tnow < tstop + 0.1 * dtime && i < iterations) {
222     // create top-level tree
223     CkPrintf("**********************************\n");
224     CkPrintf("[main] iteration: %d, tnow: %f\n", i, tnow);
225     CkPrintf("[main] rmin: (%f,%f,%f), rsize: %f\n", rmin[0], rmin[1], rmin[2], rsize);
226     CkPrintf("**********************************\n");
227 #ifndef NO_TIME
228     start = CmiWallTimer();
229     iterationStart = CmiWallTimer();
230 #endif
231     if(i == 2){
232       totalStart = CmiWallTimer();
233     }
234     CkCallback cb(CkIndex_TreePiece::doBuildTree(), pieces);
235     CkStartQD(cb);
236     init_root(-1);
237     // send roots to pieces
238     pieces.acceptRoots((CmiUInt8)topLevelRoots.getVec(), rsize, rmin[0], rmin[1], rmin[2], CkCallbackResumeThread());
239     // send root to chunk
240     chunks.acceptRoot((CmiUInt8)G_root, rmin[0], rmin[1], rmin[2], rsize, CkCallbackResumeThread());
241     
242     // update top-level nodes' moments here, since all treepieces have 
243     // completed calculating theirs
244     updateTopLevelMoments();
245 #ifndef NO_TIME
246     end = CmiWallTimer();
247 #endif
248     CkPrintf("[main] Tree building ... %f s\n", (end-start));
249 #ifdef PRINT_TREE
250     graph();
251 #endif
252 #ifndef NO_PARTITION
253 #ifndef NO_TIME
254     start = CmiWallTimer();
255 #endif
256     chunks.partition(CkCallbackResumeThread());
257 #ifndef NO_TIME
258     end = CmiWallTimer();
259 #endif
260     CkPrintf("[main] Partitioning ...  %f s\n", (end-start));
261 #endif
262 #ifndef NO_FORCES
263 #ifndef NO_TIME
264     start = CmiWallTimer();
265 #endif
266     chunks.ComputeForces(CkCallbackResumeThread());
267 #ifndef NO_TIME
268     end = CmiWallTimer();
269 #endif
270     CkPrintf("[main] Forces ...  %f s\n", (end-start));
271 #endif
272 #ifndef NO_ADVANCE
273 #ifndef NO_TIME
274     start = CmiWallTimer();
275 #endif
276     chunks.advance(CkCallbackResumeThread());
277 #ifndef NO_TIME
278     end = CmiWallTimer();
279 #endif
280     CkPrintf("[main] Advance ... %f s\n", (end-start));
281 #endif
282 #ifndef NO_CLEANUP
283 #ifndef NO_TIME
284     start = CmiWallTimer();
285 #endif
286     pieces.cleanup(CkCallbackResumeThread());
287 #ifndef NO_TIME
288     end = CmiWallTimer();
289 #endif
290     CkPrintf("[main] Clean up ... %f s\n", (end-start));
291 #endif
292     i++;
293 #ifndef NO_TIME
294     totalEnd = CmiWallTimer();
295 #endif
296     CkPrintf("[main] Total ... %f s\n", (totalEnd-iterationStart));
297 #ifndef NO_LB
298     CkPrintf("[main] starting LB\n");
299     chunks.doAtSync(CkCallbackResumeThread());
300 #endif
301
302     // must reset the vector of top level roots so that the same
303     // root isn't used over and over again by the treepieces:
304     topLevelRoots.length() = 0;
305     tnow = tnow + dtime;
306   }
307   totalEnd = CmiWallTimer();
308
309   CkPrintf("[main] Completed simulation: %f s\n", (totalEnd-totalStart));
310 #ifdef OUTPUT_ACC
311   chunks.outputAccelerations(CkCallbackResumeThread());
312 #endif
313   CkExit();
314 }
315
316
317
318 /*
319  * TAB_INIT : allocate body and cell data space
320  */
321 void 
322 Main::tab_init()
323 {
324
325   /*allocate space for personal lists of body pointers */
326   maxmybody = (real)(nbody+maxleaf*MAX_BODIES_PER_LEAF)/(real) NPROC; 
327   CkPrintf("[main] maxmybody: %d, nbody: %d, maxleaf: %d, MBPL: %d, fleaves: %f\n", maxmybody, nbody, maxleaf, MAX_BODIES_PER_LEAF, fleaves);
328   mybodytab = new bodyptr [NPROC*maxmybody]; 
329   CkAssert(mybodytab != NULL);
330   /* space is allocated so that every */
331   /* process can have a maximum of maxmybody pointers to bodies */ 
332   /* then there is an array of bodies called bodytab which is  */
333   /* allocated in the distribution generation or when the distr. */
334   /* file is read */
335   maxmycell = maxcell / NPROC;
336   maxmyleaf = maxleaf / NPROC;
337
338 }
339
340
341 void
342 Main::Help () 
343 {
344    printf("There are a total of twelve parameters, and all of them have default values.\n");
345    printf("\n");
346    printf("1) infile (char*) : The name of an input file that contains particle data.  \n");
347    printf("    The  of the file is:\n");
348    printf("\ta) An int representing the number of particles in the distribution\n");
349    printf("\tb) An int representing the dimensionality of the problem (3-D)\n");
350    printf("\tc) A double representing the current time of the simulation\n");
351    printf("\td) Doubles representing the masses of all the particles\n");
352    printf("\te) A vector (length equal to the dimensionality) of doubles\n");
353    printf("\t   representing the positions of all the particles\n");
354    printf("\tf) A vector (length equal to the dimensionality) of doubles\n");
355    printf("\t   representing the velocities of all the particles\n");
356    printf("\n");
357    printf("    Each of these numbers can be separated by any amount of whitespace.\n");
358    printf("\n");
359    printf("2) nbody (int) : If no input file is specified (the first line is blank), this\n");
360    printf("    number specifies the number of particles to generate under a plummer model.\n");
361    printf("    Default is 16384.\n");
362    printf("\n");
363    printf("3) seed (int) : The seed used by the random number generator.\n");
364    printf("    Default is 123.\n");
365    printf("\n");
366    printf("4) outfile (char*) : The name of the file that snapshots will be printed to. \n");
367    printf("    This feature has been disabled in the SPLASH release.\n");
368    printf("    Default is NULL.\n");
369    printf("\n");
370    printf("5) dtime (double) : The integration time-step.\n");
371    printf("    Default is 0.025.\n");
372    printf("\n");
373    printf("6) eps (double) : The usual potential softening\n");
374    printf("    Default is 0.05.\n");
375    printf("\n");
376    printf("7) tol (double) : The cell subdivision tolerance.\n");
377    printf("    Default is 1.0.\n");
378    printf("\n");
379    printf("8) fcells (double) : The total number of cells created is equal to \n");
380    printf("    fcells * number of leaves.\n");
381    printf("    Default is 2.0.\n");
382    printf("\n");
383    printf("9) fleaves (double) : The total number of leaves created is equal to  \n");
384    printf("    fleaves * nbody.\n");
385    printf("    Default is 0.5.\n");
386    printf("\n");
387    printf("10) tstop (double) : The time to stop integration.\n");
388    printf("    Default is 0.075.\n");
389    printf("\n");
390    printf("11) dtout (double) : The data-output interval.\n");
391    printf("    Default is 0.25.\n");
392    printf("\n");
393    printf("12) NPROC (int) : The number of processors.\n");
394    printf("    Default is 1.\n");
395 }
396
397 /*
398  * STARTRUN: startup hierarchical N-body code.
399  */
400
401 void Main::startrun()
402 {
403    int seed;
404
405    infile = getparam("in");
406    if (!infile.empty()) {
407      CkPrintf("[main] input file: %s\n", infile.c_str());
408      inputdata();
409    }
410    else {
411       nbody = getiparam("nbody");
412       if (nbody < 1) {
413          ckerr << "startrun: absurd nbody\n";
414          nbody = 16384;
415       }
416       seed = getiparam("seed");
417       if(seed < 0)
418         seed = 123;
419    }
420
421    outfile = getparam("out");
422    dtime = getdparam("dtime");
423    if(isnan(dtime))
424      dtime = 0.025;
425
426    dthf = 0.5 * dtime;
427    eps = getdparam("eps");
428    if(isnan(eps))
429      eps = 0.05;
430    epssq = eps * eps;
431
432    real epssq = eps*eps;
433    tol = getdparam("tol");
434    if(isnan(tol))
435      tol = 1.0;
436    tolsq = tol*tol;
437    
438    fcells = getdparam("fcells");
439    if(isnan(fcells))
440      fcells = 2.0;
441
442    fleaves = getdparam("fleaves");
443    if(isnan(fleaves))
444      fleaves = 15;
445
446    tstop = getdparam("tstop");
447    if(isnan(tstop))
448      tstop = 0.225;
449
450    dtout = getdparam("dtout");
451    if(isnan(dtout))
452      dtout = 0.25;
453
454    NPROC = getiparam("NPROC");
455    if(NPROC < 0)
456      NPROC = 1;
457
458    pranset(seed);
459
460    if(infile.empty()){
461     testdata();
462    }
463    setbound();
464    tout = tnow + dtout;
465 }
466
467 /*
468  * SETBOUND: Compute the initial size of the root of the tree; only done
469  * before first time step, and only processor 0 does it
470  */
471 void Main::setbound()
472 {
473    int i;
474    real side;
475    bodyptr p;
476    vector min, max;
477
478    SETVS(min,1E99);
479    SETVS(max,-1E99);
480    side=0;
481
482    for (p = bodytab; p < bodytab+nbody; p++) {
483       for (i=0; i<NDIM;i++) {
484          if (Pos(p)[i]<min[i]) min[i]=Pos(p)[i] ;
485          if (Pos(p)[i]>max[i])  max[i]=Pos(p)[i] ;
486       }
487    }
488     
489    SUBV(max,max,min);
490    for (i=0; i<NDIM;i++) if (side<max[i]) side=max[i];
491    ADDVS(rmin,min,-side/100000.0);
492    rsize = 1.00002*side;
493 }
494
495 /*
496  * TESTDATA: generate Plummer model initial conditions for test runs,
497  * scaled to units such that M = -4E = G = 1 (Henon, Hegge, etc).
498  * See Aarseth, SJ, Henon, M, & Wielen, R (1974) Astr & Ap, 37, 183.
499  */
500
501 #define MFRAC  0.999                /* mass cut off at MFRAC of total */
502
503 void Main::testdata()
504 {
505    real rsc, vsc, rsq, r, v, x, y;
506    vector cmr, cmv;
507    register bodyptr p;
508    int rejects = 0;
509    int k;
510    int halfnbody, i;
511    float offset;
512    register bodyptr cp;
513    double tmp;
514
515    tnow = 0.0;
516    bodytab = new body [nbody];
517    CkAssert(bodytab);
518    if (bodytab == NULL) {
519       ckerr << "testdata: not enuf memory\n";
520    }
521    rsc = 9 * PI / 16;
522    vsc = sqrt(1.0 / rsc);
523
524    CLRV(cmr);
525    CLRV(cmv);
526
527    halfnbody = nbody / 2;
528    if (nbody % 2 != 0) halfnbody++;
529    for (p = bodytab; p < bodytab+halfnbody; p++) {
530       Type(p) = BODY;
531       Mass(p) = 1.0 / nbody;
532       Cost(p) = 1;
533
534       r = 1 / sqrt(pow(xrand(0.0, MFRAC), -2.0/3.0) - 1);
535       /*   reject radii greater than 10 */
536       while (r > 9.0) {
537          rejects++;
538          r = 1 / sqrt(pow(xrand(0.0, MFRAC), -2.0/3.0) - 1);
539       }        
540       pickshell(Pos(p), rsc * r);
541       ADDV(cmr, cmr, Pos(p));
542       do {
543          x = xrand(0.0, 1.0);
544          y = xrand(0.0, 0.1);
545
546       } while (y > x*x * pow(1 - x*x, 3.5));
547
548       v = sqrt(2.0) * x / pow(1 + r*r, 0.25);
549       pickshell(Vel(p), vsc * v);
550       ADDV(cmv, cmv, Vel(p));
551    }
552
553    offset = 4.0;
554
555    for (p = bodytab + halfnbody; p < bodytab+nbody; p++) {
556       Type(p) = BODY;
557       Mass(p) = 1.0 / nbody;
558       Cost(p) = 1;
559
560       cp = p - halfnbody;
561       for (i = 0; i < NDIM; i++){
562          Pos(p)[i] = Pos(cp)[i] + offset; 
563          ADDV(cmr, cmr, Pos(p));
564          Vel(p)[i] = Vel(cp)[i];
565          ADDV(cmv, cmv, Vel(p));
566       }
567    }
568
569    DIVVS(cmr, cmr, (real) nbody);
570    DIVVS(cmv, cmv, (real) nbody);
571
572    for (p = bodytab; p < bodytab+nbody; p++) {
573       SUBV(Pos(p), Pos(p), cmr);
574       SUBV(Vel(p), Vel(p), cmv);
575    }
576 }
577
578 /*
579  * PICKSHELL: pick a random point on a sphere of specified radius.
580  */
581
582 void Main::pickshell(real vec[], real rad)
583 {
584    register int k;
585    double rsq, rsc;
586
587    do {
588       for (k = 0; k < NDIM; k++) {
589          vec[k] = xrand(-1.0, 1.0);
590       }
591       DOTVP(rsq, vec, vec);
592    } while (rsq > 1.0);
593
594    rsc = rad / sqrt(rsq);
595    MULVS(vec, vec, rsc);
596 }
597
598 /*
599  * GETPARAM: export version prompts user for value.
600  */
601
602 string Main::getparam(string name)
603 {
604    int i, leng;
605    string def;
606    char* temp;
607
608   for(int i = 0; i < defaults.length(); i++){
609     if(defaults[i].find(name) != string::npos){
610       int pos = defaults[i].find("=");
611       string value = defaults[i].substr(pos+1,defaults[i].length()-pos-1);
612       return value;
613     }
614   }
615   ckout << "getparam: " << name.c_str() << " unknown\n";
616   return string();
617 }
618
619 /*
620  * GETIPARAM, ..., GETDPARAM: get int, long, bool, or double parameters.
621  */
622
623 int Main::getiparam(string name)
624 {
625   string val;
626
627   val = getparam(name);
628   if(val.empty())
629     return -1;
630   else
631     return (atoi(val.c_str()));
632 }
633
634 long Main::getlparam(string name)
635 {
636   string val;
637
638   val = getparam(name);
639   if(val.empty())
640     return -1;
641   else 
642     return (atol(val.c_str()));
643 }
644
645 bool Main::getbparam(string name)
646 {
647   string val;
648
649   val = getparam(name);
650   if (strchr("tTyY1", *(val.c_str())) != 0) {
651     return (true);
652   }
653   if (strchr("fFnN0", *(val.c_str())) != 0) {
654     return (false);
655   }
656   CkPrintf("getbparam: %s=%s not bool\n", name.c_str(), val.c_str());
657 }
658
659 double Main::getdparam(string name)
660 {
661   string val;
662
663   val = getparam(name);
664   if(val.empty())
665     return NAN;
666   else 
667     return (atof(val.c_str()));
668 }
669
670 /*
671  * EXTRVALUE: extract value from name=value string.
672  */
673
674 string Main::extrvalue(string &arg)
675 {
676    char *ap;
677
678    ap = (char *) arg.c_str();
679    while (*ap != 0)
680       if (*ap++ == '=')
681          return (string(ap));
682    return (string());
683 }
684
685 void Main::updateTopLevelMoments(){
686   int depth = log8floor(numTreePieces);
687 #ifdef VERBOSE_MAIN
688   CkPrintf("[main]: updateTopLevelMoments(%d)\n", depth);
689 #endif
690   moments((nodeptr)G_root, depth);
691 #ifdef MEMCHECK
692   CkPrintf("[main] check after moments\n");
693   CmiMemoryCheck();
694 #endif
695 }
696
697 nodeptr Main::moments(nodeptr node, int depth){
698   vector tmpv;
699   if(depth == 0){
700     return node;
701   }
702   for(int i = 0; i < NSUB; i++){
703     nodeptr child = Subp(node)[i];
704     if(child != NULL){
705       nodeptr mom = moments(child, depth-1);
706 #ifdef VERBOSE_MAIN
707       CkPrintf("node 0x%x (%f,%f,%f) with node 0x%x (%d) (%f,%f,%f) mass: %f\n", node, Pos(node)[0], Pos(node)[1], Pos(node)[2], mom, i, Pos(mom)[0], Pos(mom)[1], Pos(mom)[2], Mass(mom));
708 #endif
709       Mass(node) += Mass(mom);
710       Cost(node) += Cost(mom);
711       MULVS(tmpv, Pos(mom), Mass(mom));
712       ADDV(Pos(node), Pos(node), tmpv);
713     }
714   }
715   DIVVS(Pos(node), Pos(node), Mass(node));
716 #ifdef VERBOSE_MAIN
717   CkPrintf("Pos(node): (%f,%f,%f)\n", Pos(node)[0], Pos(node)[1], Pos(node)[2]);
718 #endif
719   return node;
720 }
721
722 #ifdef PRINT_TREE
723 void Main::graph(){
724   ofstream myfile;
725   ostringstream ostr;
726
727   ostr << "tree." << nbody << "." << maxPartsPerTp << ".dot";
728   CkPrintf("[main] output file name: %s\n", ostr.str().c_str());
729   myfile.open(ostr.str().c_str());
730   myfile << "digraph tree_" << nbody << "_" << maxPartsPerTp <<" {" << endl;
731   CkQ<nodeptr> nodes(4096);
732   nodes.enq((nodeptr)G_root);
733   while(!nodes.isEmpty()){
734     nodeptr curnode = nodes.deq();
735     
736     myfile << NodeKey(curnode) << " [label=\"" << "("<< NodeKey(curnode) << ", " << Mass(curnode) << ")" << "\\n (" << Pos(curnode)[0] << "," << Pos(curnode)[1] << "," << Pos(curnode)[2] << ") " << "\"];"<< endl;
737     for(int i = 0; i < NSUB; i++){
738       nodeptr childnode = Subp(curnode)[i];
739       if(childnode != NULL){
740         if(Type(childnode) == CELL){
741           nodes.enq(childnode);
742           myfile << NodeKey(curnode) << "->" << NodeKey(childnode) << endl;
743         }
744         else if(Type(childnode) == LEAF){
745           myfile << NodeKey(curnode) << "->" << NodeKey(childnode) << endl;
746           myfile << NodeKey(childnode) << " [label=\"" << "("<< NodeKey(childnode) << ", " << Mass(childnode) << ")" << "\\n (" << Pos(childnode)[0] << "," << Pos(childnode)[1] << "," << Pos(childnode)[2] << ") num " << (((leafptr)childnode)->bodyp[0])->num << " \"];"<< endl;
747         }
748       }
749       else {
750         myfile << NodeKey(curnode) << "-> NULL_" << NodeKey(curnode) << "_" << i << endl;
751       }
752     }
753   }
754   myfile << "}" << endl;
755   myfile.close();
756 }
757 #endif
758
759 CkReductionMsg *minmax_RealVector(int nmsg, CkReductionMsg **msgs){
760   real minmax[NDIM*2];
761
762   SETVS(minmax,1E99);
763   SETVS((minmax+NDIM),-1E99);
764
765   for(int j = 0; j < nmsg; j++){
766     CkAssert(msgs[j]->getSize() == NDIM*2*sizeof(real));
767     real *tmp = (real *) msgs[j]->getData();
768     real *min = tmp;
769     real *max = tmp+NDIM;
770     for (int i = 0; i < NDIM; i++) {
771       if (minmax[i] > min[i]) {
772         minmax[i] = min[i];
773       }
774       if (minmax[NDIM+i] < max[i]) {
775         minmax[NDIM+i] = max[i];
776       }
777     }
778   }
779
780   return CkReductionMsg::buildNew(NDIM*2*sizeof(real), minmax);
781 }
782
783 void register_minmax_RealVector(){
784   minmax_RealVectorType = CkReduction::addReducer(minmax_RealVector); 
785 }
786
787 void Main::recvGlobalSizes(CkReductionMsg *msg){
788   real *tmp = (real *)msg->getData();
789   real *min = tmp;
790   real *max = tmp+NDIM;
791
792   rsize=0;
793   SUBV(max,max,min);
794   for (int i = 0; i < NDIM; i++) {
795     if (rsize < max[i]) {
796       rsize = max[i];
797     }
798   }
799   ADDVS(rmin,min,-rsize/100000.0);
800   rsize = 1.00002*rsize;
801   chunks.cleanup();
802
803   delete msg;
804 }
805 #include "barnes.def.h"