de43deee3b4c9706329c5cd6449311aad2c7c64d
[charm.git] / examples / charm++ / barnes-charm / TreePiece.cpp
1 #include "barnes.h"
2
3 TreePiece::TreePiece(CmiUInt8 p_, int which, int level_,  real rx, real ry, real rz, real rs, CkArrayIndex1D parentIdx){
4   nodeptr p;
5   isTopLevel = false;
6
7   usesAtSync = CmiFalse;
8   setMigratable(false);
9   
10   rmin[0] = rx;
11   rmin[1] = ry;
12   rmin[2] = rz;
13
14   rsize = rs;
15
16   myLevel = level_;
17   whichChildAmI = which;
18   parentIndex = parentIdx;
19   // save parent
20   parent = (nodeptr)p_;
21   haveParent = true;
22
23   mycelltab.reserve(maxmycell);
24   myleaftab.reserve(maxmyleaf);
25   ctab = new cell [maxmycell];
26   ltab = new leaf [maxmyleaf];
27
28   myRoot = NULL;
29
30   numTotalMsgs = -1;
31   numRecvdMsgs = 0;
32   haveChildren = false;
33   myNumParticles = 0;
34
35   myncell = 0;
36   mynleaf = 0;
37
38   haveCounts = false;
39   pendingChildren = 0;
40   partsToChild.resize(NSUB);
41
42   for(int i = 0; i < NSUB; i++){
43     sentTo[i] = 0;
44     childrenTreePieces[i] = -1;
45     isPending[i] = false;
46     partsToChild[i].reserve(INIT_PARTS_PER_CHILD);
47   }
48   
49   wantToSplit = false;
50
51 #ifdef MEMCHECK
52   CkPrintf("piece %d after construction\n", thisIndex);
53   CmiMemoryCheck();
54 #endif
55 }
56
57 TreePiece::TreePiece(CmiUInt8 p_, int which, int level_, CkArrayIndex1D parentIdx){
58   isTopLevel = true;
59   nodeptr p;
60   myLevel = level_;
61   // don't save parent if top-level tp. parent will be set by
62   // acceptroots
63   whichChildAmI = thisIndex;
64
65   usesAtSync = CmiFalse;
66   setMigratable(false);
67
68   mycelltab.reserve(maxmycell);
69   myleaftab.reserve(maxmyleaf);
70   ctab = new cell [maxmycell];
71   ltab = new leaf [maxmyleaf];
72
73   myRoot = NULL;
74
75   numTotalMsgs = -1;
76   numRecvdMsgs = 0;
77   haveChildren = false;
78   myNumParticles = 0;
79
80   myncell = 0;
81   mynleaf = 0;
82
83   haveCounts = false;
84   pendingChildren = 0;
85   partsToChild.resize(NSUB);
86
87   for(int i = 0; i < NSUB; i++){
88     sentTo[i] = 0;
89     childrenTreePieces[i] = -1;
90     isPending[i] = false;
91     partsToChild[i].reserve(INIT_PARTS_PER_CHILD);
92   }
93   
94   wantToSplit = false;
95
96 #ifdef MEMCHECK
97   CkPrintf("piece %d after construction\n", thisIndex);
98   CmiMemoryCheck();
99 #endif
100 }
101
102 void TreePiece::processParticles(bodyptr *particles, int num){
103   // get contents of this message
104   int xp[NDIM];
105   bodyptr p;
106
107   for(int i = 0; i < num; i++){
108     int c; // part i goes to child c
109     int relc; // index of child relative to this node (0..NSUB)
110     p = particles[i]; 
111     intcoord(xp,Pos(p),rmin,rsize);
112     relc = subindex(xp,Level(myRoot));
113     partsToChild[relc].push_back(particles[i]);
114   }
115 }
116
117 void TreePiece::sendParticlesToChildren(){
118   for(int i = 0; i < NSUB; i++){
119     int len = partsToChild[i].length();
120     if(childrenTreePieces[i] < 0 && len > 0){
121       CkAssert(!isPending[i]);
122       int child = NSUB*thisIndex+numTreePieces+i;
123       // TODO
124       // common code outside
125       // chares or chare arrays?
126       // remote method instead of message sends
127       pieces[child].insert((CmiUInt8)myRoot, i, myLevel >> 1, rmin[0], rmin[1], rmin[2], rsize, thisIndex);
128       childrenTreePieces[i] = child;
129       pendingChildren++;
130       isPending[i] = true;
131 #ifdef VERBOSE_PIECES
132       CkPrintf("piece %d inserting child %d, pending children: %d\n", thisIndex, child, pendingChildren);
133 #endif
134       // create msg from partsToChild[i], send
135       ParticleMsg *amsg = new (len) ParticleMsg;
136       amsg->num = len;
137       memcpy(amsg->particles, partsToChild[i].getVec(), len*sizeof(bodyptr));
138       pieces[child].recvParticles(amsg);
139 #ifdef VERBOSE_PIECES
140       CkPrintf("piece %d sent %d particles to child %d. sentTo[%d] = %d\n", thisIndex, len, child, i, sentTo[i]);
141 #endif
142       partsToChild[i].length() = 0;
143     }
144     else if(len > 0){
145       int child = childrenTreePieces[i];
146       if(!isPending[i]){
147         pendingChildren++;
148 #ifdef VERBOSE_PIECES
149         CkPrintf("piece %d has child %d, pending children: %d\n", thisIndex, child, pendingChildren);
150 #endif
151         isPending[i] = true;
152       }
153       pieces[child].recvRootFromParent((CmiUInt8)myRoot, rmin[0], rmin[1], rmin[2], rsize);
154       // create msg from partsToChild[i], send
155       ParticleMsg *amsg = new (len) ParticleMsg;
156       amsg->num = len;
157       memcpy(amsg->particles, partsToChild[i].getVec(), len*sizeof(bodyptr));
158       pieces[child].recvParticles(amsg);
159 #ifdef VERBOSE_PIECES
160       CkPrintf("piece %d sent %d particles to child %d. sentTo[%d] = %d\n", thisIndex, len, child, i, sentTo[i]);
161 #endif
162       partsToChild[i].length() = 0;
163     }
164   }
165 }
166
167 void TreePiece::resetPartsToChild(){
168   for(int i = 0; i < NSUB; i++){
169     partsToChild[i].length() = 0;
170   }
171 }
172
173 void TreePiece::recvRootFromParent(CmiUInt8 r, real rx, real ry, real rz, real rs){
174   parent = (nodeptr)r;
175   haveParent = true;
176   rmin[0] = rx;
177   rmin[1] = ry;
178   rmin[2] = rz;
179   rsize = rs;
180
181 #ifdef VERBOSE_PIECES
182       CkPrintf("piece %d recvd root from parent\n", thisIndex);
183 #endif
184
185   if(wantToSplit){
186 #ifdef VERBOSE_PIECES
187       CkPrintf("piece %d wanted to split\n", thisIndex);
188 #endif
189     // create own root
190     myRoot = (nodeptr) InitCell((cellptr)parent, thisIndex);
191     Subp(parent)[whichChildAmI] = myRoot;
192     Mass(myRoot) = 0.0;
193     Cost(myRoot) = 0;
194     CLRV(Pos(myRoot));
195     NodeKey(myRoot) = (NodeKey(parent) << NDIM)+whichChildAmI;
196
197     // process own particles
198     processParticles(myParticles.getVec(), myNumParticles);
199     myNumParticles = 0;
200     myParticles.length() = 0; 
201
202     // process messages
203     for(int i = 0; i < bufferedMsgs.length(); i++){
204       processParticles(bufferedMsgs[i]->particles, bufferedMsgs[i]->num);
205       delete bufferedMsgs[i];
206     }
207
208     // send out all particles
209     sendParticlesToChildren();
210   }
211 #ifdef MEMCHECK
212   CkPrintf("piece %d after recvParticles (recvRootFromParent)\n", thisIndex);
213   CmiMemoryCheck();
214 #endif
215 }
216
217 void TreePiece::acceptRoots(CmiUInt8 roots_, real rsize_, real rmx, real rmy, real rmz, CkCallback &cb){
218   rsize = rsize_;
219   rmin[0] = rmx;
220   rmin[1] = rmy;
221   rmin[2] = rmz;
222
223 #ifdef VERBOSE_PIECES
224   CkPrintf("piece %d acceptRoot rmin: (%f,%f,%f) rsize: %f\n", thisIndex, rmin[0], rmin[1], rmin[2], rsize);
225 #endif
226
227   if(isTopLevel){
228 #ifdef MEMCHECK
229   CkPrintf("piece %d before acceptRoots \n", thisIndex);
230   CmiMemoryCheck();
231 #endif
232     nodeptr *pp = (nodeptr *)roots_;
233     nodeptr p = pp[thisIndex/NSUB];
234     parent = p;
235     haveParent = true;
236 #ifdef MEMCHECK
237   CkPrintf("piece %d after acceptRoots\n", thisIndex);
238   CmiMemoryCheck();
239 #endif
240
241 #ifdef VERBOSE_PIECES
242     CkPrintf("piece [%d] acceptRoot parent 0x%x\n", thisIndex, parent);
243 #endif
244   }
245   contribute(0,0,CkReduction::concat,cb);
246 }
247
248 void TreePiece::childDone(int which){
249   pendingChildren--;
250   
251   // 'which' child just finished building its tree (and hence calculating
252   // moments. add these to your own. 
253 #ifdef VERBOSE_PIECES
254   CkPrintf("piece %d child %d done, updateMoments, pendingChildren: %d\n", thisIndex, which, pendingChildren);
255 #endif
256   updateMoments(which);
257
258   if(pendingChildren == 0){
259     DIVVS(Pos(myRoot), Pos(myRoot), Mass(myRoot));
260 #ifdef MEMCHECK
261     CkPrintf("piece %d before childDone \n", thisIndex);
262     CmiMemoryCheck();
263 #endif
264 #ifdef VERBOSE_PIECES
265     CkPrintf("piece %d all children done\n", thisIndex);
266 #endif
267     if(!isTopLevel){
268       // talk to parent
269 #ifdef VERBOSE_PIECES
270       CkPrintf("piece %d (whichChildAmI: %d) -> parent %d, i'm done\n", thisIndex, whichChildAmI, parentIndex.index);
271 #endif
272       pieces[parentIndex].childDone(whichChildAmI);
273     }
274     CkCallback cb(CkIndex_ParticleChunk::doneTreeBuild(), CkArrayIndex1D(0), chunks);
275     contribute(0,0,CkReduction::concat,cb);
276 #ifdef MEMCHECK
277     CkPrintf("piece %d after childDone\n", thisIndex);
278     CmiMemoryCheck();
279 #endif
280   }
281   else if(pendingChildren < 0){
282     CkAbort("pendingChildren < 0\n");
283   }
284 }
285
286 void TreePiece::recvParticles(ParticleMsg *msg){ 
287   bodyptr *particles = msg->particles;
288
289 #ifdef MEMCHECK
290   CkPrintf("piece %d before recvParticles \n", thisIndex);
291   CmiMemoryCheck();
292 #endif
293
294   int newtotal = myNumParticles+msg->num;
295
296   if(newtotal > maxPartsPerTp && !haveChildren){
297     // insert children into pieces array
298 #ifdef VERBOSE_PIECES
299     CkPrintf("piece %d has too many (%d+%d) particles; creating children\n", thisIndex, myNumParticles, msg->num);
300 #endif
301     if(haveParent){
302 #ifdef VERBOSE_PIECES
303       CkPrintf("piece %d has parent, process\n", thisIndex);
304 #endif
305       // first create own root
306       // then process message and own particles, to see where they should go
307       myRoot = (nodeptr) InitCell((cellptr)parent, thisIndex);
308       Subp(parent)[whichChildAmI] = myRoot;
309       Mass(myRoot) = 0.0;
310       Cost(myRoot) = 0;
311       CLRV(Pos(myRoot));
312       NodeKey(myRoot) = (NodeKey(parent) << NDIM)+whichChildAmI;
313
314       // process own particles
315       processParticles(myParticles.getVec(), myNumParticles);
316       myNumParticles = 0;
317       myParticles.length() = 0; 
318
319       // process message 
320       processParticles(msg->particles, msg->num);
321       haveChildren = true;
322
323       // send out all particles
324       sendParticlesToChildren();
325     }
326     else{
327 #ifdef VERBOSE_PIECES
328       CkPrintf("piece %d doesn't have parent, buffer msg\n", thisIndex);
329 #endif
330       // buffer message and return
331       bufferedMsgs.push_back(msg);
332       wantToSplit = true;
333       return;
334     }
335   }
336   else if(!haveChildren){
337     // stuff msg into own parts
338 #ifdef VERBOSE_PIECES
339     CkPrintf("piece %d adding %d particles to self (%d), total: %d\n", thisIndex, msg->num, myNumParticles, myNumParticles+msg->num);
340 #endif
341
342     // this is how many particles we had before expanding
343     int savedpos = myNumParticles;
344     // we now have msg->num additional particles
345     myNumParticles += msg->num;
346     // expand array of particles to include new ones
347     myParticles.resize(myNumParticles);
348     // this is where we start copying new particles to 
349     bodyptr *savedstart = myParticles.getVec()+savedpos;
350
351     memcpy(savedstart, msg->particles, (msg->num)*sizeof(bodyptr));
352   }
353   else{
354     // have no particles of own, send to children
355 #ifdef VERBOSE_PIECES
356     CkPrintf("piece %d has children\n", thisIndex);
357 #endif
358     // at this point, we have a list of particles 
359     // destined for each child
360     processParticles(msg->particles, msg->num);
361     sendParticlesToChildren();
362     
363   }
364
365   delete msg;
366 #ifdef MEMCHECK
367   CkPrintf("piece %d after recvParticles\n", thisIndex);
368   CmiMemoryCheck();
369 #endif
370 }
371
372 void TreePiece::doBuildTree(){
373   if(haveChildren){
374 #ifdef VERBOSE_PIECES
375     CkPrintf("piece %d fake doneTreeBuild()\n", thisIndex);
376 #endif
377   }
378   else{
379     // don't have children, build own tree
380 #ifdef VERBOSE_PIECES
381     CkPrintf("piece %d doesn't have children, building tree\n", thisIndex);
382 #endif
383     buildTree();
384     hackcofm(0,thisIndex);
385     if(!isTopLevel && haveParent){
386       // once you've built your own tree, 
387       // you must notify your parent that you're done
388 #ifdef VERBOSE_PIECES
389       CkPrintf("piece %d real !topLevel doneTreeBuild(), haveParent: %d\n", thisIndex, haveParent);
390 #endif
391       pieces[parentIndex].childDone(whichChildAmI);
392     }
393 #ifdef VERBOSE_PIECES
394     else{
395       CkPrintf("piece %d real topLevel doneTreeBuild()\n", thisIndex);
396     }
397 #endif
398     CkCallback cb(CkIndex_ParticleChunk::doneTreeBuild(), CkArrayIndex1D(0), chunks);
399     contribute(0,0,CkReduction::concat,cb);
400   }
401 #ifdef MEMCHECK
402   CkPrintf("piece %d after doBuildTree()\n", thisIndex);
403   CmiMemoryCheck();
404 #endif
405 }
406
407 void TreePiece::buildTree(){
408   bodyptr p, *pp;
409   int ProcessId = thisIndex;
410   nodeptr Current_Root;
411
412 #ifdef MEMCHECK
413   CkPrintf("piece %d before buildTree\n", thisIndex);
414   CmiMemoryCheck();
415 #endif
416   // start from parent 
417   // this way, particles can be added to Root until it
418   // needs to be split.
419   Current_Root = (nodeptr) parent;
420   done = 0;
421   for(int i = 0; i < myParticles.length(); i++){
422     //CkPrintf("[%d] inserting particle %d, current_root: %ld\n", thisIndex, myParticles[i]->num, NodeKey(Current_Root));
423     (nodeptr) loadtree(myParticles[i], (cellptr) Current_Root, ProcessId);
424     //Current_Root = (nodeptr) loadtree(myParticles[i], (cellptr) Current_Root, ProcessId);
425     done++;
426   }
427
428 #ifdef MEMCHECK
429   CkPrintf("piece %d after buildTree\n", thisIndex);
430   CmiMemoryCheck();
431 #endif
432 }
433
434 nodeptr TreePiece::loadtree(bodyptr p, cellptr root, unsigned int ProcessId){
435   int l, xq[NDIM], xp[NDIM], flag;
436   int i, j, root_level;
437   bool valid_root;
438   int kidIndex;
439   nodeptr *qptr, mynode;
440   cellptr c;
441   leafptr le;
442
443 #ifdef MEMCHECK
444   CkPrintf("piece %d before loadtree\n", thisIndex);
445   CmiMemoryCheck();
446 #endif
447   intcoord(xp, Pos(p),rmin,rsize);
448   valid_root = TRUE;
449   mynode = (nodeptr) root;
450   kidIndex = subindex(xp, Level(mynode));
451   qptr = &Subp(mynode)[kidIndex];
452
453   l = Level(mynode) >> 1;
454
455   flag = TRUE;
456   while (flag) {                           /* loop descending tree     */
457     if (l == 0) {
458       //CkPrintf("[%d] not enough levels in tree, %d done\n", thisIndex, done);
459       CkAbort("tree depth\n");
460     }
461     if (*qptr == NULL) { 
462         le = InitLeaf((cellptr) mynode, ProcessId);
463         ParentOf(p) = (nodeptr) le;
464         Level(p) = l;
465         ChildNum(p) = le->num_bodies;
466         ChildNum(le) = kidIndex;
467         NodeKey(le) = (NodeKey(mynode) << NDIM) + kidIndex;
468         //CkPrintf("[%d] new leaf at kidindex: %d, key: %ld\n", thisIndex, kidIndex, NodeKey(le));
469         Bodyp(le)[le->num_bodies++] = p;
470
471         *qptr = (nodeptr) le;
472         flag = FALSE;
473     }
474     if (flag && *qptr && (Type(*qptr) == LEAF)) {
475       /*   reached a "leaf"?      */
476         le = (leafptr) *qptr;
477         //CkPrintf("[%d] reached existing leaf at kidindex %d key: %ld\n", thisIndex, kidIndex, NodeKey(le));
478         if (le->num_bodies == MAX_BODIES_PER_LEAF) {
479           //CkPrintf("[%d] too many particles. splitting\n", thisIndex);
480           *qptr = (nodeptr) SubdivideLeaf(le, (cellptr) mynode, l, ProcessId);
481         }
482         else {
483           //CkPrintf("[%d] enough space\n", thisIndex);
484           ParentOf(p) = (nodeptr) le;
485           Level(p) = l;
486           ChildNum(p) = le->num_bodies;
487           Bodyp(le)[le->num_bodies++] = p;
488           flag = FALSE;
489         }
490     }
491     if (flag) {
492       mynode = *qptr;
493       kidIndex = subindex(xp, l);
494       qptr = &Subp(*qptr)[kidIndex];  /* move down one level  */
495       l = l >> 1;                            /* and test next bit    */
496     }
497   }
498 #ifdef MEMCHECK
499   CkPrintf("piece %d after loadtree\n", thisIndex);
500   CmiMemoryCheck();
501 #endif
502   return ParentOf((leafptr) *qptr);
503
504 }
505
506 leafptr TreePiece::InitLeaf(cellptr parent, unsigned ProcessId)
507 {
508    leafptr l;
509    int i, Mycell;
510
511    l = makeleaf(ProcessId);
512    l->next = NULL;
513    l->prev = NULL;
514    if (parent==NULL)
515       Level(l) = IMAX >> 1;
516    else
517       Level(l) = Level(parent) >> 1;
518    ParentOf(l) = (nodeptr) parent;
519    ChildNum(l) = 0;
520
521    return (l);
522 }
523
524 /*
525  * MAKECELL: allocation routine for cells.
526  */
527
528 cellptr TreePiece::makecell(unsigned ProcessId)
529 {
530    cellptr c;
531    int i, Mycell;
532     
533    if (myncell == maxmycell) {
534       CkPrintf("makecell: Proc %d needs more than %d cells; increase fcells\n", 
535             ProcessId,maxmycell);
536       CkAbort("makecell\n");
537    }
538    Mycell = myncell++;
539    c = ctab + Mycell;
540    Type(c) = CELL;
541    Done(c) = FALSE;
542    Mass(c) = 0.0;
543    for (i = 0; i < NSUB; i++) {
544       Subp(c)[i] = NULL;
545    }
546    mycelltab.push_back(c);
547    return (c);
548 }
549
550 /*
551  * MAKELEAF: allocation routine for leaves.
552  */
553
554 leafptr TreePiece::makeleaf(unsigned ProcessId)
555 {
556    leafptr le;
557    int i, Myleaf;
558     
559    if (mynleaf == maxmyleaf) {
560       CkPrintf("makeleaf: Proc %d needs more than %d leaves; increase fleaves\n",
561             ProcessId,maxmyleaf);
562       CkAbort("makeleaf\n");
563    }
564    Myleaf = mynleaf++;
565    le = ltab + Myleaf;
566    le->seqnum = ProcessId * maxmyleaf + Myleaf;
567    Type(le) = LEAF;
568    Done(le) = FALSE;
569    Mass(le) = 0.0;
570    le->num_bodies = 0;
571    for (i = 0; i < MAX_BODIES_PER_LEAF; i++) {
572       Bodyp(le)[i] = NULL;
573    }
574    myleaftab.push_back(le);
575    return (le);
576 }
577
578 cellptr
579 TreePiece::SubdivideLeaf (leafptr le, cellptr parent_, unsigned int l, unsigned int ProcessId)
580 {
581    cellptr c;
582    int i, index;
583    int xp[NDIM];
584    bodyptr bodies[MAX_BODIES_PER_LEAF];
585    int num_bodies;
586    bodyptr p;
587
588    /* first copy leaf's bodies to temp array, so we can reuse the leaf */
589    num_bodies = le->num_bodies;
590    for (i = 0; i < num_bodies; i++) {
591       bodies[i] = Bodyp(le)[i];
592       Bodyp(le)[i] = NULL;
593    }
594    le->num_bodies = 0;
595    /* create the parent cell for this subtree */
596    c = InitCell(parent_, ProcessId);
597    ChildNum(c) = ChildNum(le);
598    NodeKey(c) = (NodeKey(parent_) << NDIM)+ChildNum(c);
599    //CkPrintf("new cell key: %ld\n", NodeKey(c));
600    /* do first particle separately, so we can reuse le */
601    p = bodies[0];
602    intcoord(xp, Pos(p),rmin,rsize);
603    index = subindex(xp, l);
604    Subp(c)[index] = (nodeptr) le;
605    ChildNum(le) = index;
606    ParentOf(le) = (nodeptr) c;
607    NodeKey(le) = (NodeKey(c) << NDIM)+index; 
608    //CkPrintf("existing leaf (gets particle %d) key: %ld\n", p->num, NodeKey(le));
609
610    Level(le) = l >> 1;
611    /* set stuff for body */
612    ParentOf(p) = (nodeptr) le;
613    ChildNum(p) = le->num_bodies;
614    Level(p) = l >> 1;
615    /* insert the body */
616    Bodyp(le)[le->num_bodies++] = p;
617    /* now handle the rest */
618    for (i = 1; i < num_bodies; i++) {
619       p = bodies[i];
620       intcoord(xp, Pos(p),rmin,rsize);
621       index = subindex(xp, l);
622       if (!Subp(c)[index]) {
623          le = InitLeaf(c, ProcessId);
624          ChildNum(le) = index;
625          Subp(c)[index] = (nodeptr) le;
626          NodeKey(le) = (NodeKey(c) << NDIM) + index;
627          //CkPrintf("i: %d, created leaf index %d (gets particle %d) key: %ld\n", i, index, p->num, NodeKey(le));
628       }
629       else {
630          le = (leafptr) Subp(c)[index];
631          //CkPrintf("i: %d, existing leaf index %d (gets particle %d) key: %ld\n", i, index, p->num, NodeKey(le));
632       }
633       ParentOf(p) = (nodeptr) le;
634       ChildNum(p) = le->num_bodies;
635       Level(p) = l >> 1;
636       Bodyp(le)[le->num_bodies++] = p;
637    }
638    return c;
639 }
640
641 cellptr TreePiece::InitCell(cellptr parent_, unsigned ProcessId)
642 {
643   cellptr c;
644   int i, Mycell;
645
646   c = makecell(ProcessId);
647   c->next = NULL;
648   c->prev = NULL;
649   if (parent_ == NULL)
650     Level(c) = IMAX >> 1;
651   else
652     Level(c) = Level(parent_) >> 1;
653   ParentOf(c) = (nodeptr) parent_;
654   ChildNum(c) = 0;
655    
656   return (c);
657 }
658
659 /*
660  * SUBINDEX: determine which subcell to select.
661  */
662
663 int subindex(int x[NDIM], int l)
664 {
665    int i, k;
666    int yes;
667     
668    i = 0;
669    yes = FALSE;
670    if (x[0] & l) {
671       i += NSUB >> 1;
672       yes = TRUE;
673    }
674    for (k = 1; k < NDIM; k++) {
675       if (((x[k] & l) && !yes) || (!(x[k] & l) && yes)) { 
676          i += NSUB >> (k + 1);
677          yes = TRUE;
678       }
679       else yes = FALSE;
680    }
681
682    return (i);
683 }
684
685 /* * INTCOORD: compute integerized coordinates.  * Returns: TRUE
686 unless rp was out of bounds.  */
687
688 bool intcoord(int *xp, vector rp, vector rmin, real rsize)
689 {
690    int k;
691    bool inb;
692    double xsc;
693    double tmp;
694     
695    inb = TRUE;
696    for (k = 0; k < NDIM; k++) {
697       xsc = (rp[k] - rmin[k]) / rsize; 
698       if (0.0 <= xsc && xsc < 1.0) {
699         tmp = IMAX * xsc;
700          xp[k] = (int)tmp;
701       }
702       else {
703          inb = FALSE;
704       }
705    }
706    return (inb);
707 }
708
709 /*
710  * HACKCOFM: descend tree finding center-of-mass coordinates.
711  */
712
713 void TreePiece::hackcofm(int nc, unsigned ProcessId)
714 {
715    int i,Myindex;
716    nodeptr r;
717    leafptr l;
718    leafptr* ll;
719    bodyptr p;
720    cellptr q;
721    cellptr *cc;
722    vector tmpv, dr;
723    real drsq;
724    matrix drdr, Idrsq, tmpm;
725
726    /* get a cell using get*sub.  Cells are got in reverse of the order in */
727    /* the cell array; i.e. reverse of the order in which they were created */
728    /* this way, we look at child cells before parents                    */
729     
730    for (ll = myleaftab.getVec() + mynleaf - 1; ll >= myleaftab.getVec(); ll--) {
731       l = *ll;
732       Mass(l) = 0.0;
733       Cost(l) = 0;
734       CLRV(Pos(l));
735       for (i = 0; i < l->num_bodies; i++) {
736          p = Bodyp(l)[i];
737          Mass(l) += Mass(p);
738          Cost(l) += Cost(p);
739          MULVS(tmpv, Pos(p), Mass(p));
740          ADDV(Pos(l), Pos(l), tmpv);
741       }
742       DIVVS(Pos(l), Pos(l), Mass(l));
743 #ifdef QUADPOLE
744       CLRM(Quad(l));
745       for (i = 0; i < l->num_bodies; i++) {
746          p = Bodyp(l)[i];
747          SUBV(dr, Pos(p), Pos(l));
748          OUTVP(drdr, dr, dr);
749          DOTVP(drsq, dr, dr);
750          SETMI(Idrsq);
751          MULMS(Idrsq, Idrsq, drsq);
752          MULMS(tmpm, drdr, 3.0);
753          SUBM(tmpm, tmpm, Idrsq);
754          MULMS(tmpm, tmpm, Mass(p));
755          ADDM(Quad(l), Quad(l), tmpm);
756       }
757 #endif
758       Done(l)=TRUE;
759    }
760    for (cc = mycelltab.getVec()+myncell-1; cc >= mycelltab.getVec(); cc--) {
761       q = *cc;
762       Mass(q) = 0.0;
763       Cost(q) = 0;
764       CLRV(Pos(q));
765       for (i = 0; i < NSUB; i++) {
766          r = Subp(q)[i];
767          if (r != NULL) {
768             Mass(q) += Mass(r);
769             Cost(q) += Cost(r);
770             MULVS(tmpv, Pos(r), Mass(r));
771             ADDV(Pos(q), Pos(q), tmpv);
772             Done(r) = FALSE;
773          }
774       }
775       DIVVS(Pos(q), Pos(q), Mass(q));
776 #ifdef QUADPOLE
777       CLRM(Quad(q));
778       for (i = 0; i < NSUB; i++) {
779          r = Subp(q)[i];
780          if (r != NULL) {
781             SUBV(dr, Pos(r), Pos(q));
782             OUTVP(drdr, dr, dr);
783             DOTVP(drsq, dr, dr);
784             SETMI(Idrsq);
785             MULMS(Idrsq, Idrsq, drsq);
786             MULMS(tmpm, drdr, 3.0);
787             SUBM(tmpm, tmpm, Idrsq);
788             MULMS(tmpm, tmpm, Mass(r));
789             ADDM(tmpm, tmpm, Quad(r));
790             ADDM(Quad(q), Quad(q), tmpm);
791          }
792       }
793 #endif
794    }
795 }
796
797 // called only treepiece has children, so that myRoot will
798 // be valid.
799 void TreePiece::updateMoments(int which){
800
801   int i,Myindex;
802   nodeptr r;
803   leafptr l;
804   leafptr* ll;
805   bodyptr p;
806   cellptr q;
807   cellptr *cc;
808   vector tmpv, dr;
809   real drsq;
810   matrix drdr, Idrsq, tmpm;
811
812   // add moments of child 'which' to meRoot
813   q = (cellptr) myRoot;
814   r = Subp(q)[which];
815   if (r != NULL){
816     Mass(q) += Mass(r);
817     Cost(q) += Cost(r);
818     MULVS(tmpv, Pos(r), Mass(r));
819     ADDV(Pos(q), Pos(q), tmpv);
820
821 #ifdef QUADPOLE
822     CLRM(Quad(q));
823     r = Subp(q)[which];
824     CkAssert (r != NULL);
825     SUBV(dr, Pos(r), Pos(q));
826     OUTVP(drdr, dr, dr);
827     DOTVP(drsq, dr, dr);
828     SETMI(Idrsq);
829     MULMS(Idrsq, Idrsq, drsq);
830     MULMS(tmpm, drdr, 3.0);
831     SUBM(tmpm, tmpm, Idrsq);
832     MULMS(tmpm, tmpm, Mass(r));
833     ADDM(tmpm, tmpm, Quad(r));
834     ADDM(Quad(q), Quad(q), tmpm);
835 #endif
836   }
837 }
838
839 void TreePiece::cleanup(CkCallback &cb_){
840
841   // reset your parent's pointer to your topmost 
842   // cell/leaf. every treepiece has a parent node
843
844   Subp(parent)[whichChildAmI] = NULL;
845
846   numTotalMsgs = -1;
847   numRecvdMsgs = 0;
848   haveChildren = false;
849   haveCounts = false;
850   for(int i = 0; i < NSUB; i++){
851     sentTo[i] = 0;
852     partsToChild[i].length() = 0;
853     isPending[i] = false;
854   }
855   // treepieces need not have myRoots
856   // since we begin inserting particles as 
857   // children of the parent node. therefore,
858   // unless more than one bucket was created 
859   // during the treebuild process, or the node
860   // was split because it was too fat, myRoot
861   // will be NULL
862   if(myRoot != NULL){
863     for(int i = 0; i < NSUB; i++){
864       Subp(myRoot)[i] = NULL;
865     }
866   }
867   myParticles.length() = 0;
868   myNumParticles = 0;
869   pendingChildren = 0;
870   mycelltab.length() = 0;
871   myleaftab.length() = 0;
872   myncell = 0;
873   mynleaf = 0;
874   haveParent = false;
875   bufferedMsgs.length() = 0;
876   wantToSplit = false;
877
878   contribute(0,0,CkReduction::concat,cb_);
879 }
880