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