examples/charm++/barnes-charm/src/load.C

   1 /*************************************************************************/
   2 /*                                                                       */
   3 /*  Copyright (c) 1994 Stanford University                               */
   4 /*                                                                       */
   5 /*  All rights reserved.                                                 */
   6 /*                                                                       */
   7 /*  Permission is given to use, copy, and modify this software for any   */
   8 /*  non-commercial purpose as long as this copyright notice is not       */
   9 /*  removed.  All other uses, including redistribution in whole or in    */
  10 /*  part, are forbidden without prior written permission.                */
  11 /*                                                                       */
  12 /*  This software is provided with absolutely no warranty and no         */
  13 /*  support.                                                             */
  14 /*                                                                       */
  15 /*************************************************************************/
  16
  17 EXTERN_ENV
  18 #define global extern
  19
  20 #include "code.h"
  21 #include "defs.h"
  22
  23 bool intcoord();
  24 cellptr makecell(unsigned int ProcessId);
  25 leafptr makeleaf(unsigned int ProcessId);
  26 cellptr SubdivideLeaf(leafptr le, cellptr parent, unsigned int l,
  27                       unsigned int ProcessId);
  28
  29 cellptr InitCell(cellptr parent, unsigned int ProcessId);
  30 leafptr InitLeaf(cellptr parent, unsigned int ProcessId);
  31 nodeptr loadtree(bodyptr p, cellptr root, unsigned int ProcessId);
  32
  33 /*
  34  * MAKETREE: initialize tree structure for hack force calculation.
  35  */
  36
  37 maketree(ProcessId)
  38    unsigned ProcessId;
  39 {
  40    bodyptr p, *pp;
  41
  42    Local[ProcessId].myncell = 0;
  43    Local[ProcessId].mynleaf = 0;
  44    if (ProcessId == 0) {
  45       Local[ProcessId].mycelltab[Local[ProcessId].myncell++] = Global->G_root;
  46    }
  47    Local[ProcessId].Current_Root = (nodeptr) Global->G_root;
  48    for (pp = Local[ProcessId].mybodytab;
  49         pp < Local[ProcessId].mybodytab+Local[ProcessId].mynbody; pp++) {
  50       p = *pp;
  51       if (Mass(p) != 0.0) {
  52          Local[ProcessId].Current_Root
  53             = (nodeptr) loadtree(p, (cellptr) Local[ProcessId].Current_Root,
  54                                  ProcessId);
  55       }
  56       else {
  57          LOCK(Global->io_lock);
  58          fprintf(stderr, "Process %d found body %d to have zero mass\n",
  59                  ProcessId, (int) p);
  60          UNLOCK(Global->io_lock);
  61       }
  62    }
  63    BARRIER(Global->Bartree,NPROC);
  64    hackcofm( 0, ProcessId );
  65    BARRIER(Global->Barcom,NPROC);
  66 }
  67
  68 cellptr InitCell(parent, ProcessId)
  69    cellptr parent;
  70    unsigned ProcessId;
  71 {
  72    cellptr c;
  73    int i, Mycell;
  74
  75    c = makecell(ProcessId);
  76    c->processor = ProcessId;
  77    c->next = NULL;
  78    c->prev = NULL;
  79    if (parent == NULL)
  80       Level(c) = IMAX >> 1;
  81    else
  82       Level(c) = Level(parent) >> 1;
  83    Parent(c) = (nodeptr) parent;
  84    ChildNum(c) = 0;
  85    return (c);
  86 }
  87
  88 leafptr InitLeaf(parent, ProcessId)
  89    cellptr parent;
  90    unsigned ProcessId;
  91 {
  92    leafptr l;
  93    int i, Mycell;
  94
  95    l = makeleaf(ProcessId);
  96    l->processor = ProcessId;
  97    l->next = NULL;
  98    l->prev = NULL;
  99    if (parent==NULL)
 100       Level(l) = IMAX >> 1;
 101    else
 102       Level(l) = Level(parent) >> 1;
 103    Parent(l) = (nodeptr) parent;
 104    ChildNum(l) = 0;
 105    return (l);
 106 }
 107
 108 printtree (n)
 109    nodeptr n;
 110 {
 111    int k;
 112    cellptr c;
 113    leafptr l;
 114    bodyptr p;
 115    nodeptr tmp;
 116    unsigned long nseq;
 117    int xp[NDIM];
 118
 119    switch (Type(n)) {
 120     case CELL:
 121       c = (cellptr) n;
 122       nseq = c->seqnum;
 123       printf("Cell : Cost = %d, ", Cost(c));
 124       PRTV("Pos", Pos(n));
 125       printf("\n");
 126       for (k = 0; k < NSUB; k++) {
 127          printf("Child #%d: ", k);
 128          if (Subp(c)[k] == NULL) {
 129             printf("NONE");
 130          }
 131          else {
 132             if (Type(Subp(c)[k]) == CELL) {
 133                nseq = ((cellptr) Subp(c)[k])->seqnum;
 134                printf("C: Cost = %d, ", Cost(Subp(c)[k]));
 135             }
 136             else {
 137                nseq = ((leafptr) Subp(c)[k])->seqnum;
 138                printf("L: # Bodies = %2d, Cost = %d, ",
 139                       ((leafptr) Subp(c)[k])->num_bodies, Cost(Subp(c)[k]));
 140             }
 141             tmp = Subp(c)[k];
 142             PRTV("Pos", Pos(tmp));
 143          }
 144          printf("\n");
 145       }
 146       for (k=0;k<NSUB;k++) {
 147          if (Subp(c)[k] != NULL) {
 148             printtree(Subp(c)[k]);
 149          }
 150       }
 151       break;
 152     case LEAF:
 153       l = (leafptr) n;
 154       nseq = l->seqnum;
 155       printf("Leaf : # Bodies = %2d, Cost = %d, ", l->num_bodies, Cost(l));
 156       PRTV("Pos", Pos(n));
 157       printf("\n");
 158       for (k = 0; k < l->num_bodies; k++) {
 159          p = Bodyp(l)[k];
 160          printf("Body #%2d: Num = %2d, Level = %o, ",
 161                 p - bodytab, k, Level(p));
 162          PRTV("Pos",Pos(p));
 163          printf("\n");
 164       }
 165       break;
 166     default:
 167       fprintf(stderr, "Bad type\n");
 168       exit(-1);
 169       break;
 170    }
 171    fflush(stdout);
 172 }
 173
 174 /*
 175  * LOADTREE: descend tree and insert particle.
 176  */
 177
 178 nodeptr
 179 loadtree(p, root, ProcessId)
 180    bodyptr p;                        /* body to load into tree */
 181    cellptr root;
 182    unsigned ProcessId;
 183 {
 184    int l, xq[NDIM], xp[NDIM], xor[NDIM], subindex(), flag;
 185    int i, j, root_level;
 186    bool valid_root;
 187    int kidIndex;
 188    volatile nodeptr *volatile qptr, mynode;
 189    cellptr c;
 190    leafptr le;
 191
 192    intcoord(xp, Pos(p));
 193    valid_root = TRUE;
 194    for (i = 0; i < NDIM; i++) {
 195       xor[i] = xp[i] ^ Local[ProcessId].Root_Coords[i];
 196    }
 197    for (i = IMAX >> 1; i > Level(root); i >>= 1) {
 198       for (j = 0; j < NDIM; j++) {
 199          if (xor[j] & i) {
 200             valid_root = FALSE;
 201             break;
 202          }
 203       }
 204       if (!valid_root) {
 205          break;
 206       }
 207    }
 208    if (!valid_root) {
 209       if (root != Global->G_root) {
 210          root_level = Level(root);
 211          for (j = i; j > root_level; j >>= 1) {
 212             root = (cellptr) Parent(root);
 213          }
 214          valid_root = TRUE;
 215          for (i = IMAX >> 1; i > Level(root); i >>= 1) {
 216             for (j = 0; j < NDIM; j++) {
 217                if (xor[j] & i) {
 218                   valid_root = FALSE;
 219                   break;
 220                }
 221             }
 222             if (!valid_root) {
 223                printf("P%d body %d\n", ProcessId, p - bodytab);
 224                root = Global->G_root;
 225             }
 226          }
 227       }
 228    }
 229    root = Global->G_root;
 230    mynode = (nodeptr) root;
 231    kidIndex = subindex(xp, Level(mynode));
 232    qptr = &Subp(mynode)[kidIndex];
 233
 234    l = Level(mynode) >> 1;
 235
 236    flag = TRUE;
 237    while (flag) {                           /* loop descending tree     */
 238       if (l == 0) {
 239          error("not enough levels in tree\n");
 240       }
 241       if (*qptr == NULL) {
 242          /* lock the parent cell */
 243          ALOCK(CellLock->CL, ((cellptr) mynode)->seqnum % MAXLOCK);
 244          if (*qptr == NULL) {
 245             le = InitLeaf((cellptr) mynode, ProcessId);
 246             Parent(p) = (nodeptr) le;
 247             Level(p) = l;
 248             ChildNum(p) = le->num_bodies;
 249             ChildNum(le) = kidIndex;
 250             Bodyp(le)[le->num_bodies++] = p;
 251             *qptr = (nodeptr) le;
 252             flag = FALSE;
 253          }
 254          AULOCK(CellLock->CL, ((cellptr) mynode)->seqnum % MAXLOCK);
 255          /* unlock the parent cell */
 256       }
 257       if (flag && *qptr && (Type(*qptr) == LEAF)) {
 258          /*   reached a "leaf"?      */
 259          ALOCK(CellLock->CL, ((cellptr) mynode)->seqnum % MAXLOCK);
 260          /* lock the parent cell */
 261          if (Type(*qptr) == LEAF) {             /* still a "leaf"?      */
 262             le = (leafptr) *qptr;
 263             if (le->num_bodies == MAX_BODIES_PER_LEAF) {
 264                *qptr = (nodeptr) SubdivideLeaf(le, (cellptr) mynode, l,
 265                                                   ProcessId);
 266             }
 267             else {
 268                Parent(p) = (nodeptr) le;
 269                Level(p) = l;
 270                ChildNum(p) = le->num_bodies;
 271                Bodyp(le)[le->num_bodies++] = p;
 272                flag = FALSE;
 273             }
 274          }
 275          AULOCK(CellLock->CL, ((cellptr) mynode)->seqnum % MAXLOCK);
 276          /* unlock the node           */
 277       }
 278       if (flag) {
 279          mynode = *qptr;
 280          kidIndex = subindex(xp, l);
 281          qptr = &Subp(*qptr)[kidIndex];  /* move down one level  */
 282          l = l >> 1;                            /* and test next bit    */
 283       }
 284    }
 285    SETV(Local[ProcessId].Root_Coords, xp);
 286    return Parent((leafptr) *qptr);
 287 }
 288 \f
 289
 290 /* * INTCOORD: compute integerized coordinates.  * Returns: TRUE
 291 unless rp was out of bounds.  */
 292
 293 bool intcoord(xp, rp)
 294   int xp[NDIM];                  /* integerized coordinate vector [0,IMAX) */
 295   vector rp;                     /* real coordinate vector (system coords) */
 296 {
 297    int k;
 298    bool inb;
 299    double xsc;
 300    double tmp;
 301
 302    inb = TRUE;
 303    for (k = 0; k < NDIM; k++) {
 304       xsc = (rp[k] - Global->rmin[k]) / Global->rsize;
 305       if (0.0 <= xsc && xsc < 1.0) {
 306         tmp = IMAX * xsc;
 307          xp[k] = (int)tmp;
 308       }
 309       else {
 310          inb = FALSE;
 311       }
 312    }
 313    return (inb);
 314 }
 315
 316 /*
 317  * SUBINDEX: determine which subcell to select.
 318  */
 319
 320 int subindex(x, l)
 321   int x[NDIM];                       /* integerized coordinates of particle */
 322   int l;                             /* current level of tree */
 323 {
 324    int i, k;
 325    int yes;
 326
 327    i = 0;
 328    yes = FALSE;
 329    if (x[0] & l) {
 330       i += NSUB >> 1;
 331       yes = TRUE;
 332    }
 333    for (k = 1; k < NDIM; k++) {
 334       if (((x[k] & l) && !yes) || (!(x[k] & l) && yes)) {
 335          i += NSUB >> (k + 1);
 336          yes = TRUE;
 337       }
 338       else yes = FALSE;
 339    }
 340
 341    return (i);
 342 }
 343
 344 \f
 345
 346 /*
 347  * HACKCOFM: descend tree finding center-of-mass coordinates.
 348  */
 349
 350 hackcofm(nc, ProcessId)
 351   int nc;
 352   unsigned ProcessId;
 353 {
 354    int i,Myindex;
 355    nodeptr r;
 356    leafptr l;
 357    leafptr* ll;
 358    bodyptr p;
 359    cellptr q;
 360    cellptr *cc;
 361    vector tmpv, dr;
 362    real drsq;
 363    matrix drdr, Idrsq, tmpm;
 364
 365    /* get a cell using get*sub.  Cells are got in reverse of the order in */
 366    /* the cell array; i.e. reverse of the order in which they were created */
 367    /* this way, we look at child cells before parents                    */
 368
 369    for (ll = Local[ProcessId].myleaftab + Local[ProcessId].mynleaf - 1;
 370         ll >= Local[ProcessId].myleaftab; ll--) {
 371       l = *ll;
 372       Mass(l) = 0.0;
 373       Cost(l) = 0;
 374       CLRV(Pos(l));
 375       for (i = 0; i < l->num_bodies; i++) {
 376          p = Bodyp(l)[i];
 377          Mass(l) += Mass(p);
 378          Cost(l) += Cost(p);
 379          MULVS(tmpv, Pos(p), Mass(p));
 380          ADDV(Pos(l), Pos(l), tmpv);
 381       }
 382       DIVVS(Pos(l), Pos(l), Mass(l));
 383 #ifdef QUADPOLE
 384       CLRM(Quad(l));
 385       for (i = 0; i < l->num_bodies; i++) {
 386          p = Bodyp(l)[i];
 387          SUBV(dr, Pos(p), Pos(l));
 388          OUTVP(drdr, dr, dr);
 389          DOTVP(drsq, dr, dr);
 390          SETMI(Idrsq);
 391          MULMS(Idrsq, Idrsq, drsq);
 392          MULMS(tmpm, drdr, 3.0);
 393          SUBM(tmpm, tmpm, Idrsq);
 394          MULMS(tmpm, tmpm, Mass(p));
 395          ADDM(Quad(l), Quad(l), tmpm);
 396       }
 397 #endif
 398       Done(l)=TRUE;
 399    }
 400    for (cc = Local[ProcessId].mycelltab+Local[ProcessId].myncell-1;
 401         cc >= Local[ProcessId].mycelltab; cc--) {
 402       q = *cc;
 403       Mass(q) = 0.0;
 404       Cost(q) = 0;
 405       CLRV(Pos(q));
 406       for (i = 0; i < NSUB; i++) {
 407          r = Subp(q)[i];
 408          if (r != NULL) {
 409             while(!Done(r)) {
 410                /* wait */
 411             }
 412             Mass(q) += Mass(r);
 413             Cost(q) += Cost(r);
 414             MULVS(tmpv, Pos(r), Mass(r));
 415             ADDV(Pos(q), Pos(q), tmpv);
 416             Done(r) = FALSE;
 417          }
 418       }
 419       DIVVS(Pos(q), Pos(q), Mass(q));
 420 #ifdef QUADPOLE
 421       CLRM(Quad(q));
 422       for (i = 0; i < NSUB; i++) {
 423          r = Subp(q)[i];
 424          if (r != NULL) {
 425             SUBV(dr, Pos(r), Pos(q));
 426             OUTVP(drdr, dr, dr);
 427             DOTVP(drsq, dr, dr);
 428             SETMI(Idrsq);
 429             MULMS(Idrsq, Idrsq, drsq);
 430             MULMS(tmpm, drdr, 3.0);
 431             SUBM(tmpm, tmpm, Idrsq);
 432             MULMS(tmpm, tmpm, Mass(r));
 433             ADDM(tmpm, tmpm, Quad(r));
 434             ADDM(Quad(q), Quad(q), tmpm);
 435          }
 436       }
 437 #endif
 438       Done(q)=TRUE;
 439    }
 440 }
 441 \f
 442 cellptr
 443 SubdivideLeaf (le, parent, l, ProcessId)
 444    leafptr le;
 445    cellptr parent;
 446    unsigned int l;
 447    unsigned int ProcessId;
 448 {
 449    cellptr c;
 450    int i, index;
 451    int xp[NDIM];
 452    bodyptr bodies[MAX_BODIES_PER_LEAF];
 453    int num_bodies;
 454    bodyptr p;
 455
 456    /* first copy leaf's bodies to temp array, so we can reuse the leaf */
 457    num_bodies = le->num_bodies;
 458    for (i = 0; i < num_bodies; i++) {
 459       bodies[i] = Bodyp(le)[i];
 460       Bodyp(le)[i] = NULL;
 461    }
 462    le->num_bodies = 0;
 463    /* create the parent cell for this subtree */
 464    c = InitCell(parent, ProcessId);
 465    ChildNum(c) = ChildNum(le);
 466    /* do first particle separately, so we can reuse le */
 467    p = bodies[0];
 468    intcoord(xp, Pos(p));
 469    index = subindex(xp, l);
 470    Subp(c)[index] = (nodeptr) le;
 471    ChildNum(le) = index;
 472    Parent(le) = (nodeptr) c;
 473    Level(le) = l >> 1;
 474    /* set stuff for body */
 475    Parent(p) = (nodeptr) le;
 476    ChildNum(p) = le->num_bodies;
 477    Level(p) = l >> 1;
 478    /* insert the body */
 479    Bodyp(le)[le->num_bodies++] = p;
 480    /* now handle the rest */
 481    for (i = 1; i < num_bodies; i++) {
 482       p = bodies[i];
 483       intcoord(xp, Pos(p));
 484       index = subindex(xp, l);
 485       if (!Subp(c)[index]) {
 486          le = InitLeaf(c, ProcessId);
 487          ChildNum(le) = index;
 488          Subp(c)[index] = (nodeptr) le;
 489       }
 490       else {
 491          le = (leafptr) Subp(c)[index];
 492       }
 493       Parent(p) = (nodeptr) le;
 494       ChildNum(p) = le->num_bodies;
 495       Level(p) = l >> 1;
 496       Bodyp(le)[le->num_bodies++] = p;
 497    }
 498    return c;
 499 }
 500
 501 /*
 502  * MAKECELL: allocation routine for cells.
 503  */
 504
 505 cellptr makecell(ProcessId)
 506    unsigned ProcessId;
 507 {
 508    cellptr c;
 509    int i, Mycell;
 510
 511    if (Local[ProcessId].mynumcell == maxmycell) {
 512       error("makecell: Proc %d needs more than %d cells; increase fcells\n",
 513             ProcessId,maxmycell);
 514    }
 515    Mycell = Local[ProcessId].mynumcell++;
 516    c = Local[ProcessId].ctab + Mycell;
 517    c->seqnum = ProcessId*maxmycell+Mycell;
 518    Type(c) = CELL;
 519    Done(c) = FALSE;
 520    Mass(c) = 0.0;
 521    for (i = 0; i < NSUB; i++) {
 522       Subp(c)[i] = NULL;
 523    }
 524    Local[ProcessId].mycelltab[Local[ProcessId].myncell++] = c;
 525    return (c);
 526 }
 527
 528 /*
 529  * MAKELEAF: allocation routine for leaves.
 530  */
 531
 532 leafptr makeleaf(ProcessId)
 533    unsigned ProcessId;
 534 {
 535    leafptr le;
 536    int i, Myleaf;
 537
 538    if (Local[ProcessId].mynumleaf == maxmyleaf) {
 539       error("makeleaf: Proc %d needs more than %d leaves; increase fleaves\n",
 540             ProcessId,maxmyleaf);
 541    }
 542    Myleaf = Local[ProcessId].mynumleaf++;
 543    le = Local[ProcessId].ltab + Myleaf;
 544    le->seqnum = ProcessId * maxmyleaf + Myleaf;
 545    Type(le) = LEAF;
 546    Done(le) = FALSE;
 547    Mass(le) = 0.0;
 548    le->num_bodies = 0;
 549    for (i = 0; i < MAX_BODIES_PER_LEAF; i++) {
 550       Bodyp(le)[i] = NULL;
 551    }
 552    Local[ProcessId].myleaftab[Local[ProcessId].mynleaf++] = le;
 553    return (le);
 554 }
 555
 556