Blockstructured Adaptive Mesh Refinement in object-oriented C++
Blockstructured Adaptive Mesh Refinement in object-oriented C++
00001 00006 #include "DAGHDistribution.h" 00007 00008 #ifdef DEBUG_PRINT 00009 #include "CommServer.h" 00010 #endif 00011 00012 /*************************************************************************/ 00013 /* Partitioning functions based on Paul Walker's implementation */ 00014 /*************************************************************************/ 00015 void partition_one(const BBox& wholebbox, BBox* boxes, const int np, 00016 const int dim); 00017 void partition_two(const BBox& wholebbox, BBox* boxes, const int np, 00018 const int dim1, const int dim2); 00019 void partition_all(const BBox& wholebbox, BBox* boxes, const int np); 00020 void pwalker_partition(const BBox& wholebbox, BBox* boxes, const int np); 00021 /*************************************************************************/ 00022 00023 DAGHDistribution::DAGHDistribution(const int disttype, BBoxList& bbl) 00024 : type(disttype), boxes(0), wfunc(0) 00025 { 00026 boxes = new BBox[bbl.number()]; 00027 register BBox* bb = 0; 00028 register int cnt = 0; 00029 for (bb=bbl.first();bb;bb=bbl.next()) { 00030 boxes[cnt++] = *bb; 00031 } 00032 } 00033 00034 void DAGHDistribution::set_dist(const int disttype, BBoxList& bbl) 00035 { 00036 type = disttype; 00037 if (boxes) delete [] boxes; 00038 boxes = new BBox[bbl.number()]; 00039 register BBox* bb = 0; 00040 register int cnt = 0; 00041 for (bb=bbl.first();bb;bb=bbl.next()) { 00042 boxes[cnt++] = *bb; 00043 } 00044 } 00045 00046 unsigned long DAGHDistribution::get_load(GridUnitList& ggul, 00047 const short* olap) 00048 { 00049 if (wfunc) { 00050 unsigned long load = 0, n = 0; 00051 int levels = ggul.levels(); 00052 for (int l=0;l<levels;l++) 00053 load += (*wfunc)(&(n=ggul.numelems(l,olap)),&l); 00054 return load; 00055 } 00056 else return( ggul.load(olap) ); 00057 } 00058 00059 void DAGHDistribution::partition(GridUnitList*& ggul, 00060 GridUnitList*& lgul, 00061 const int np, const int me, 00062 const int minw, const short* olap) 00063 { 00064 if (type == DAGHCompositeDistribution) 00065 partitionComposite(ggul, lgul, np, me, minw, olap); 00066 00067 else if (type == DAGHBlockXDistribution) 00068 partitionBlockOne(DAGH_X, ggul, lgul, np, me, minw, olap); 00069 else if (type == DAGHBlockYDistribution) 00070 partitionBlockOne(DAGH_Y, ggul, lgul, np, me, minw, olap); 00071 else if (type == DAGHBlockZDistribution) 00072 partitionBlockOne(DAGH_Z, ggul, lgul, np, me, minw, olap); 00073 00074 else if (type == DAGHBlockXYDistribution) 00075 partitionBlockTwo(DAGH_X, DAGH_Y, ggul, lgul, np, me, minw, olap); 00076 else if (type == DAGHBlockYZDistribution) 00077 partitionBlockTwo(DAGH_Y, DAGH_Z, ggul, lgul, np, me, minw, olap); 00078 else if (type == DAGHBlockXZDistribution) 00079 partitionBlockTwo(DAGH_X, DAGH_Z, ggul, lgul, np, me, minw, olap); 00080 00081 else if (type == DAGHBlockAllDistribution) 00082 partitionBlockAll(ggul, lgul, np, me, minw, olap); 00083 00084 else if (type == DAGHUserDefDistribution) 00085 partitionUserDef(ggul, lgul, np, me, minw, olap); 00086 00087 else 00088 assert(0); 00089 } 00090 00091 void DAGHDistribution::partition(const BBox& wbbox, 00092 BBox& mybbox, 00093 const int np, 00094 const int me) 00095 { 00096 if (type == DAGHCompositeDistribution) 00097 type = DAGHBlockXDistribution; 00098 00099 if (me >= np) { 00100 mybbox = BBox::_empty_bbox; 00101 return; 00102 } 00103 00104 if (type != DAGHUserDefDistribution) { 00105 if (boxes) delete [] boxes; 00106 boxes = new BBox[np]; 00107 00108 if (type == DAGHBlockXDistribution) 00109 partition_one(wbbox, boxes, np, DAGH_X); 00110 else if (type == DAGHBlockYDistribution) 00111 partition_one(wbbox, boxes, np, DAGH_Y); 00112 else if (type == DAGHBlockZDistribution) 00113 partition_one(wbbox, boxes, np, DAGH_Z); 00114 00115 else if (type == DAGHBlockXYDistribution) 00116 partition_two(wbbox, boxes, np, DAGH_X, DAGH_Y); 00117 else if (type == DAGHBlockYZDistribution) 00118 partition_two(wbbox, boxes, np, DAGH_Y, DAGH_Z); 00119 else if (type == DAGHBlockXZDistribution) 00120 partition_two(wbbox, boxes, np, DAGH_X, DAGH_Z); 00121 00122 else if (type == DAGHBlockAllDistribution) 00123 partition_all(wbbox, boxes, np); 00124 00125 else 00126 assert(0); 00127 } 00128 mybbox = boxes[me]; 00129 } 00130 00131 /* This function is slightly modified. Perhaps, segmentation faults are created 00132 by the complex GridUnit-modifications. No errors have been observed after 00133 lgul->add() and operator= have been introduced. RD */ 00134 00135 void DAGHDistribution::partitionComposite(GridUnitList*& ggul, 00136 GridUnitList*& lgul, 00137 const int np, const int me, 00138 const int minw, const short* olap) 00139 { 00140 #ifdef DEBUG_PRINT 00141 /* Cannot partition an empty list */ 00142 assert(ggul->number() > 0); 00143 00144 /* Zero procs... ? */ 00145 assert(np > 0); 00146 #endif 00147 00148 //if (!localspans.isempty()) localspans.empty(); 00149 00150 if (me < np) { 00151 if (lgul) lgul->empty(); 00152 else { 00153 lgul = new GridUnitList; 00154 } 00155 } 00156 00157 // dMapIndex *lmin = 0; 00158 // dMapIndex lmax; 00159 00160 unsigned long work = get_load(*ggul,olap); 00161 00162 unsigned long thresh = (unsigned long) (1.0*(work/np)+0.5); 00163 register unsigned long curthresh = thresh; 00164 00165 int p = 0, pcnt = 0; 00166 unsigned long tmpw = 0, w = 0; 00167 00168 register GridUnit *tmpgu = ggul->first(); 00169 00170 #ifdef DEBUG_PRINT_DIST 00171 ( comm_service::log() << "\n{ Composite Distribution: " 00172 << np << " " << work << " " << thresh 00173 << "} \n" 00174 ).flush(); 00175 #endif 00176 00177 for (int cnt=0;tmpgu;cnt++,tmpgu=ggul->next()) { 00178 w += (tmpw = tmpgu->guWork(olap)); 00179 00180 if (p>=np) { 00181 tmpgu->guSetOwner(np-1); 00182 if (np-1 == me) { 00183 GridUnit* newgu = lgul->add(); 00184 *newgu = *tmpgu; 00185 } 00186 continue; 00187 } 00188 00189 #ifdef DEBUG_PRINT_DIST 00190 ( comm_service::log() << "\n{ Partitioning: " 00191 << cnt << " " << p << " " 00192 << curthresh << " " << w << " " << tmpw << " " 00193 << tmpgu->guBaseIndex() << " " 00194 << *ggul->currec() 00195 << "} \n" 00196 00197 ).flush(); 00198 #endif 00199 00200 if (w <= curthresh) { 00201 00202 #ifdef DEBUG_PRINT_DIST 00203 if (p>=np) 00204 ( comm_service::log() << "\n************* Current Composite List *************\n" 00205 << *ggul 00206 << "\n************* ***************** *************\n" 00207 << "\n************* Current Bucket *************\n" 00208 << *((SimpleBucketVoid *) ggul) 00209 << "\n************* ***************** *************\n" 00210 ).flush(); 00211 #endif 00212 #ifdef DEBUG_PRINT 00213 assert (p<np); 00214 #endif 00215 tmpgu->guSetOwner(p); pcnt++; 00216 if (p == me) { 00217 GridUnit* newgu = lgul->add(); 00218 *newgu = *tmpgu; 00219 } 00220 continue; 00221 // if (!lmin) lmin = new dMapIndex(tmpgu->guBaseIndex()); 00222 // lmax = tmpgu->guBaseIndex(); 00223 } 00224 else 00225 if (w-curthresh < tmpw/4) { 00226 00227 #ifdef DEBUG_PRINT_DIST 00228 if (p>=np) 00229 ( comm_service::log() << "\n************* Current Composite List *************\n" 00230 << *ggul 00231 << "\n************* ***************** *************\n" 00232 << "\n************* Current Bucket *************\n" 00233 << *((SimpleBucketVoid *) ggul) 00234 << "\n************* ***************** *************\n" 00235 ).flush(); 00236 #endif 00237 #ifdef DEBUG_PRINT 00238 assert(p<np); 00239 #endif 00240 tmpgu->guSetOwner(p); pcnt++; 00241 if (p == me) { 00242 GridUnit* newgu = lgul->add(); 00243 *newgu = *tmpgu; 00244 } 00245 // lmax = tmpgu->guBaseIndex(); 00246 p++; pcnt = 0; 00247 if (p<np) 00248 curthresh = w + (unsigned long) (1.0*((work>w ? work-w : 0)/(np-p))+0.5); 00249 continue; 00250 } 00251 else 00252 if (w >= curthresh && tmpgu->guExtent(tmpgu->guCrsLev()) > minw) { 00253 w -= tmpw; 00254 cnt--; 00255 unsigned long oldwork = work; 00256 ggul->decompose(); 00257 work = get_load(*ggul,olap); 00258 00259 // if (work == oldwork) assert(0); 00260 00261 #ifdef DEBUG_PRINT_DIST 00262 ( comm_service::log() << "\n************* After Decompose *************\n" 00263 << "Cur Proc:" << p << " " 00264 << *ggul->current() << "\n" 00265 << "Cur GU:" << ggul->current() << "\n" 00266 << "\n************* ***************** *************\n" 00267 ).flush(); 00268 #endif 00269 thresh = (unsigned long) (1.0*(work/np)+0.5); 00270 curthresh = curthresh+(unsigned long)(1.0*((work>oldwork ? 00271 work-oldwork : 0)/(np-p))+0.5); 00272 continue; 00273 } 00274 else 00275 if (pcnt == 0) { 00276 00277 #ifdef DEBUG_PRINT_DIST 00278 if (p>=np) 00279 ( comm_service::log() << "\n************* Current Composite List *************\n" 00280 << *ggul 00281 << "\n************* ***************** *************\n" 00282 << "\n************* Current Bucket *************\n" 00283 << *((SimpleBucketVoid *) ggul) 00284 << "\n************* ***************** *************\n" 00285 ).flush(); 00286 #endif 00287 #ifdef DEBUG_PRINT 00288 assert(p<np); 00289 #endif 00290 tmpgu->guSetOwner(p); pcnt++; 00291 if (p == me) { 00292 GridUnit* newgu = lgul->add(); 00293 *newgu = *tmpgu; 00294 } 00295 // lmax = tmpgu->guBaseIndex(); 00296 p++; pcnt = 0; 00297 if (p<np) 00298 curthresh = w + (unsigned long) (1.0*((work>w ? 00299 work-w : 0)/(np-p))+0.5); 00300 continue; 00301 } 00302 else { 00303 p++; pcnt = 0; 00304 00305 #ifdef DEBUG_PRINT_DIST 00306 if (p>=np) 00307 ( comm_service::log() << "\n************* Current Composite List *************\n" 00308 << *ggul 00309 << "\n************* ***************** *************\n" 00310 << "\n************* Current Bucket *************\n" 00311 << *((SimpleBucketVoid *) ggul) 00312 << "\n************* ***************** *************\n" 00313 ).flush(); 00314 #endif 00315 #ifdef DEBUG_PRINT 00316 assert(p<np); 00317 #endif 00318 if (p>=np) continue; 00319 00320 tmpgu->guSetOwner(p); 00321 if (p == me) { 00322 GridUnit* newgu = lgul->add(); 00323 *newgu = *tmpgu; 00324 } 00325 // if (!lmin) lmin = new dMapIndex(tmpgu->guBaseIndex()); 00326 // lmax = tmpgu->guBaseIndex(); 00327 curthresh = (w-tmpw) + (unsigned long) (1.0*((work+tmpw>w ? 00328 work+tmpw-w : 0)/(np-p))+0.5); 00329 if (w >= curthresh) { 00330 p++; pcnt = 0; 00331 if (p<np) 00332 curthresh = w + (unsigned long) (1.0*((work>w ? 00333 work-w : 0)/(np-p))+0.5); 00334 } 00335 continue; 00336 } 00337 } 00338 if (p<me && me<np) { 00339 cerr << "Cannot distribute anything to processor " << me 00340 << ". Refine base grid!\n"; 00341 } 00342 } 00343 00344 void DAGHDistribution::partitionBlockOne(const int axis, 00345 GridUnitList*& ggul, 00346 GridUnitList*& lgul, 00347 const int np, const int me, 00348 const int minw, const short* olap) 00349 { 00350 #ifdef DEBUG_PRINT 00351 /* Cannot partition an empty list */ 00352 assert(ggul->number() > 0); 00353 00354 /* Zero procs... ? */ 00355 assert(np > 0); 00356 #endif 00357 00358 if (boxes) delete [] boxes; 00359 boxes = new BBox[np]; 00360 00361 const int flev = (ggul->first())->guFineLev(); 00362 BBoxList tmpbbl; 00363 ggul->bboxlist(tmpbbl, flev, 0, 0); 00364 BBox wbbox = tmpbbl.reduce(); 00365 partition_one(wbbox, boxes, np, axis); 00366 00367 GridUnitList tmpgul; 00368 for (register int p=0;p<np;p++) { 00369 ggul->intersect(boxes[p], flev, tmpgul, 0); 00370 GridUnitList tmpgul(*ggul); 00371 tmpgul *= boxes[p]; 00372 ggul->setowner(tmpgul,p); 00373 if (p == me) { 00374 if (lgul) { 00375 lgul->empty(); 00376 *lgul = tmpgul; 00377 } 00378 else { 00379 lgul = new GridUnitList(tmpgul); 00380 } 00381 lgul->setowner(me); 00382 } 00383 } 00384 } 00385 00386 void DAGHDistribution::partitionBlockTwo(const int axis1, const int axis2, 00387 GridUnitList*& ggul, 00388 GridUnitList*& lgul, 00389 const int np, const int me, 00390 const int minw, const short* olap) 00391 { 00392 #ifdef DEBUG_PRINT 00393 /* Cannot partition an empty list */ 00394 assert(ggul->number() > 0); 00395 00396 /* Zero procs... ? */ 00397 assert(np > 0); 00398 #endif 00399 00400 if (boxes) delete [] boxes; 00401 boxes = new BBox[np]; 00402 00403 const int flev = ggul->first()->guFineLev(); 00404 BBoxList tmpbbl; 00405 00406 ggul->bboxlist(tmpbbl, flev, 0, 0); 00407 BBox wbbox = tmpbbl.reduce(); 00408 partition_two(wbbox, boxes, np, axis1, axis2); 00409 GridUnitList tmpgul; 00410 for (register int p=0;p<np;p++) { 00411 ggul->intersect(boxes[p], flev, tmpgul, 0); 00412 GridUnitList tmpgul(*ggul); 00413 tmpgul *= boxes[p]; 00414 ggul->setowner(tmpgul,p); 00415 if (p == me) { 00416 if (lgul) { 00417 lgul->empty(); 00418 *lgul = tmpgul; 00419 } 00420 else { 00421 lgul = new GridUnitList(tmpgul); 00422 } 00423 lgul->setowner(me); 00424 } 00425 } 00426 } 00427 00428 void DAGHDistribution::partitionBlockAll(GridUnitList*& ggul, 00429 GridUnitList*& lgul, 00430 const int np, const int me, 00431 const int minw, const short* olap) 00432 { 00433 #ifdef DEBUG_PRINT 00434 /* Cannot partition an empty list */ 00435 assert(ggul->number() > 0); 00436 00437 /* Zero procs... ? */ 00438 assert(np > 0); 00439 #endif 00440 00441 if (boxes) delete [] boxes; 00442 boxes = new BBox[np]; 00443 00444 const int flev = (ggul->first())->guFineLev(); 00445 BBoxList tmpbbl; 00446 ggul->bboxlist(tmpbbl, flev, 0, 0); 00447 BBox wbbox = tmpbbl.reduce(); 00448 partition_all(wbbox, boxes, np); 00449 00450 GridUnitList tmpgul; 00451 for (register int p=0;p<np;p++) { 00452 ggul->intersect(boxes[p], flev, tmpgul, 0); 00453 GridUnitList tmpgul(*ggul); 00454 tmpgul *= boxes[p]; 00455 ggul->setowner(tmpgul,p); 00456 if (p == me) { 00457 if (lgul) { 00458 lgul->empty(); 00459 *lgul = tmpgul; 00460 } 00461 else { 00462 lgul = new GridUnitList(tmpgul); 00463 } 00464 lgul->setowner(me); 00465 } 00466 } 00467 } 00468 00469 void DAGHDistribution::partitionUserDef(GridUnitList*& ggul, 00470 GridUnitList*& lgul, 00471 const int np, const int me, 00472 const int minw, const short* olap) 00473 { 00474 #ifdef DEBUG_PRINT 00475 /* The list boxes has to be defined */ 00476 assert(boxes); 00477 00478 /* Cannot partition an empty list */ 00479 assert(ggul->number() > 0); 00480 00481 /* Zero procs... ? */ 00482 assert(np > 0); 00483 #endif 00484 00485 const int flev = (ggul->first())->guFineLev(); 00486 BBoxList tmpbbl; 00487 ggul->bboxlist(tmpbbl, flev, 0, 0); 00488 BBox wbbox = tmpbbl.reduce(); 00489 partition_all(wbbox, boxes, np); 00490 00491 GridUnitList tmpgul; 00492 for (register int p=0;p<np;p++) { 00493 ggul->intersect(boxes[p], flev, tmpgul, 0); 00494 GridUnitList tmpgul(*ggul); 00495 tmpgul *= boxes[p]; 00496 ggul->setowner(tmpgul,p); 00497 if (p == me) { 00498 if (lgul) { 00499 lgul->empty(); 00500 *lgul = tmpgul; 00501 } 00502 else { 00503 lgul = new GridUnitList(tmpgul); 00504 } 00505 lgul->setowner(me); 00506 } 00507 } 00508 } 00509 00510 /*$int DAGHDistribution::islocal(dMapIndex const &idx) 00511 { 00512 DAGHSpan *span = 0; 00513 00514 DAGHListLoop(localspans,span,DAGHSpan) 00515 if (span->contains(idx)) return 1; 00516 DAGHEndLoop 00517 00518 return 0; 00519 }$*/ 00520 00521 void partition_one(const BBox& wholebbox, BBox* boxes, const int np, 00522 const int dim) 00523 { 00524 00525 const int extent = wholebbox.extents(dim); 00526 00527 #ifdef DEBUG_PRINT 00528 assert (extent >= np); // I should have atleast extent points. */ 00529 #endif 00530 00531 const int dx = extent/np; 00532 int Xtra = extent-(dx*np); 00533 00534 for (register int p=0;p<np;p++) { 00535 const int dlower = p*dx + Min(Xtra,p); 00536 const int dupper = (np-p-1)*dx + Max((Xtra-p-1),0); 00537 boxes[p] = wholebbox; 00538 boxes[p].growlower(dim,dlower); 00539 boxes[p].growupper(dim,-dupper); 00540 } 00541 } 00542 00543 /* A modification of pwalker_partition by Paul Walker */ 00544 void partition_two(const BBox& wholebbox, BBox* boxes, const int np, 00545 const int dim1, const int dim2) 00546 { 00547 /* Great now figure out the decomposition in n1,n2 */ 00548 int n1,n2; 00549 00550 /* Initizialize */ 00551 n1 = np; 00552 n2 = 1; 00553 if (n1 % 2 == 0) { 00554 while (n1 % 2 == 0 && n1 > n2 00555 && (abs(n1-n2) >= abs(n1/2-n2*2))) { 00556 n1 /= 2; 00557 n2 *= 2; 00558 } 00559 } else { 00560 while (n1 % 3 == 0 && n1 > n2 00561 && (abs(n1-n2) >= abs(n1/3-n2*3))) { 00562 n1 /= 3; 00563 n2 *= 3; 00564 } 00565 } 00566 00567 const int rank = wholebbox.rank; 00568 const Coords& glb(wholebbox.lower()); 00569 const Coords& gub(wholebbox.upper()); 00570 const Coords& step(wholebbox.stepsize()); 00571 00572 /* Set up deltas */ 00573 Coords delt(rank,0); 00574 delt(dim1) = (gub(dim1)-glb(dim1)+step(dim1))/n1; 00575 delt(dim2) = (gub(dim2)-glb(dim2)+step(dim2))/n2; 00576 00577 /* Set up correction vector */ 00578 Coords cor(rank,0); 00579 cor(0) = (gub(0)-glb(0)+step(0))%n1; 00580 cor(1) = (gub(1)-glb(1)+step(1))%n2; 00581 00582 /* Initialize the bboxes */ 00583 for (register int p=0;p<np;p++) boxes[p] = wholebbox; 00584 00585 for (int ii=0;ii<n1;ii++) { 00586 for (int jj=0;jj<n2;jj++) { 00587 /* We may care about topology here later... */ 00588 int whichp = ii + jj*n1; 00589 00590 Coords c1(rank,0), c2 (rank,0); 00591 00592 c1(0) = (cor(0) > ii) ? 1 : 0; 00593 c1(1) = (cor(1) > jj) ? 1 : 0; 00594 00595 /*$c2(0) = n1 - (cor(0) - ii) - 1; 00596 c2(1) = n2 - (cor(1) - jj) - 1;$*/ 00597 00598 c2(0) = Min(cor(0),ii); 00599 c2(1) = Min(cor(1),jj); 00600 00601 Coords& lb(boxes[whichp].lower()); 00602 Coords& ub(boxes[whichp].upper()); 00603 00604 lb(0) = glb(0) + (ii * (delt(0)) + c2(0)); 00605 lb(1) = glb(1) + (jj * (delt(1)) + c2(1)); 00606 00607 ub(0) = lb(0) + delt(0) + c1(0) - 1; 00608 ub(1) = lb(1) + delt(1) + c1(1) - 1; 00609 } 00610 } 00611 } 00612 00613 /* A modification of pwalker_partition by Paul Walker */ 00614 void partition_all(const BBox& wholebbox, BBox* boxes, const int np) 00615 { 00616 /* Great now figure out the decomposition in n1,n2,n3 */ 00617 int n1,n2,n3; 00618 00619 /* Initizialize */ 00620 n1 = np; 00621 n2 = 1; 00622 n3 = 1; 00623 for (register int i=0;i<2;i++) { 00624 int twof = 0; 00625 if (n1 % 2 == 0) twof = 1; 00626 if (twof) { 00627 while (n1 % 2 == 0 && n1 > n2 00628 && (abs(n1-n2) >= abs(n1/2-n2*2))) { 00629 n1 /= 2; 00630 n2 *= 2; 00631 } 00632 } else { 00633 while (n1 % 3 == 0 && n1 > n2 00634 && (abs(n1-n2) >= abs(n1/3-n2*3))) { 00635 n1 /= 3; 00636 n2 *= 3; 00637 } 00638 } 00639 twof = 0; 00640 if (n2 % 2 == 0) twof = 1; 00641 if (twof) { 00642 while (n2 % 2 == 0 && n2 > n3 00643 && (abs(n2-n3) >= abs(n2/2-n3*2))) { 00644 n2 /= 2; 00645 n3 *= 2; 00646 } 00647 } else { 00648 while (n2 % 3 == 0 && n2 > n3 00649 && (abs(n2-n3) >= abs(n2/3-n3*3))) { 00650 n2 /= 3; 00651 n3 *= 3; 00652 } 00653 } 00654 } 00655 /* Sort based on global box (later)... */ 00656 00657 const int rank = wholebbox.rank; 00658 const Coords& glb(wholebbox.lower()); 00659 const Coords& gub(wholebbox.upper()); 00660 const Coords& step(wholebbox.stepsize()); 00661 00662 /* Set up deltas */ 00663 Coords delt(rank,0); 00664 delt(0) = (gub(0)-glb(0)+step(0))/n1; 00665 delt(1) = (gub(1)-glb(1)+step(1))/n2; 00666 delt(2) = (gub(2)-glb(2)+step(2))/n3; 00667 00668 /* Set up correction vector */ 00669 Coords cor(rank,0); 00670 cor(0) = (gub(0)-glb(0)+step(0))%n1; 00671 cor(1) = (gub(1)-glb(1)+step(1))%n2; 00672 cor(2) = (gub(2)-glb(2)+step(2))%n3; 00673 00674 /*$comm_service::log() << "glb " << glb << endl; 00675 comm_service::log() << "gub " << gub << endl; 00676 comm_service::log() << "deltas " << delt << endl; 00677 comm_service::log() << "correction " << cor << endl;$*/ 00678 00679 /* Initialize the bboxes */ 00680 for (register int p=0;p<np;p++) boxes[p] = wholebbox; 00681 00682 for (int ii=0;ii<n1;ii++) { 00683 for (int jj=0;jj<n2;jj++) { 00684 for (int kk=0;kk<n3;kk++) { 00685 /* We may care about topology here later... */ 00686 int whichp = ii + jj*n1 + kk*n1*n2; 00687 00688 /*$comm_service::log() << whichp << " " 00689 << ii << " " 00690 << jj << " " 00691 << kk << endl;$*/ 00692 00693 Coords c1(rank,0), c2 (rank,0); 00694 00695 c1(0) = (cor(0) > ii) ? 1 : 0; 00696 c1(1) = (cor(1) > jj) ? 1 : 0; 00697 c1(2) = (cor(2) > kk) ? 1 : 0; 00698 00699 /*$c2(0) = n1 - (cor(0) - ii) - 1; 00700 c2(1) = n2 - (cor(1) - jj) - 1; 00701 c2(2) = n3 - (cor(2) - kk) - 1;$*/ 00702 00703 c2(0) = Min(cor(0),ii); 00704 c2(1) = Min(cor(1),jj); 00705 c2(2) = Min(cor(2),kk); 00706 00707 Coords& lb(boxes[whichp].lower()); 00708 Coords& ub(boxes[whichp].upper()); 00709 00710 lb(0) = glb(0) + (ii * (delt(0)) + c2(0)); 00711 lb(1) = glb(1) + (jj * (delt(1)) + c2(1)); 00712 lb(2) = glb(2) + (kk * (delt(2)) + c2(2)); 00713 00714 ub(0) = lb(0) + delt(0) + c1(0) - 1; 00715 ub(1) = lb(1) + delt(1) + c1(1) - 1; 00716 ub(2) = lb(2) + delt(2) + c1(2) - 1; 00717 00718 /*$comm_service::log() << "lb " << lb << endl; 00719 comm_service::log() << "ub " << ub << endl;$*/ 00720 00721 } 00722 } 00723 } 00724 } 00725 00726 void pwalker_partition(const BBox& wholebbox, BBox* boxes, const int np) 00727 { 00728 Coords glb(wholebbox.lower()); 00729 Coords gub(wholebbox.upper()); 00730 /* Great now figure out the decomposition in nx,ny,nz */ 00731 int nx,ny,nz; 00732 00733 /* Initizialize */ 00734 nx = np; 00735 ny = 1; 00736 nz = 1; 00737 for (int i=0;i<2;i++) { 00738 int twof = 0; 00739 if (nx % 2 == 0) twof = 1; 00740 if (twof) { 00741 while (nx % 2 == 0 && nx > ny 00742 && (abs(nx-ny) >= abs(nx/2-ny*2))) { 00743 nx /= 2; 00744 ny *= 2; 00745 } 00746 } else { 00747 while (nx % 3 == 0 && nx > ny 00748 && (abs(nx-ny) >= abs(nx/3-ny*3))) { 00749 nx /= 3; 00750 ny *= 3; 00751 } 00752 } 00753 twof = 0; 00754 if (ny % 2 == 0) twof = 1; 00755 if (twof) { 00756 while (ny % 2 == 0 && ny > nz 00757 && (abs(ny-nz) >= abs(ny/2-nz*2))) { 00758 ny /= 2; 00759 nz *= 2; 00760 } 00761 } else { 00762 while (ny % 3 == 0 && ny > nz 00763 && (abs(ny-nz) >= abs(ny/3-nz*3))) { 00764 ny /= 3; 00765 nz *= 3; 00766 } 00767 } 00768 } 00769 /* Sort based on global box (later)... */ 00770 00771 /* Now make the coordinates for each box ... */ 00772 BBox mt(3,0); 00773 for (int p=0;p<np;p++) boxes[p] = mt; 00774 00775 int ii,jj,kk; 00776 int *delt = new int[3], *layout = new int[3]; 00777 00778 /* Set up layout and deltas */ 00779 layout[0] = nx; layout[1] = ny; layout[2] = nz; 00780 for (int mm=0;mm<3;mm++) { 00781 delt[mm] = (gub(mm)-glb(mm))/layout[mm]; 00782 } 00783 00784 Coords mylb(3,0),myub(3,0); 00785 for (ii=0;ii<nx;ii++) 00786 for (jj=0;jj<ny;jj++) 00787 for (kk=0;kk<nz;kk++) { 00788 /* We may care about topology here later... */ 00789 int whichp = ii + jj*nx + kk*nx*ny; 00790 mylb(0) = delt[0]*ii; 00791 mylb(1) = delt[1]*jj; 00792 mylb(2) = delt[2]*kk; 00793 00794 myub(0) = 00795 (ii != nx-1 ? delt[0]*(ii+1)-1 : gub(0)); 00796 myub(1) = 00797 (jj != ny-1 ? delt[1]*(jj+1)-1 : gub(1)); 00798 myub(2) = 00799 (kk != nz-1 ? delt[2]*(kk+1)-1 : gub(2)); 00800 00801 boxes[whichp].setlower(mylb); 00802 boxes[whichp].setupper(myub); 00803 boxes[whichp].setstepsize(0,1); 00804 boxes[whichp].setstepsize(1,1); 00805 boxes[whichp].setstepsize(2,1); 00806 } 00807 }
Quickstart Users Guide Programmers Reference Installation Examples Download
AMROC Main Home Contactlast update: 06/01/04