45 #ifndef GEOMETRICGENERATOR 46 #define GEOMETRICGENERATOR 48 #include <Teuchos_Comm.hpp> 49 #include <Teuchos_ParameterList.hpp> 50 #include <Teuchos_FilteredIterator.hpp> 51 #include <Teuchos_ParameterEntry.hpp> 60 #include <Tpetra_MultiVector_decl.hpp> 63 #include <Teuchos_ArrayViewDecl.hpp> 64 #include <Teuchos_RCP.hpp> 65 #include <Tpetra_Distributor.hpp> 76 using Teuchos::CommandLineProcessor;
77 using Teuchos::ArrayView;
79 using Teuchos::ArrayRCP;
82 #define CATCH_EXCEPTIONS(pp) \ 83 catch (std::runtime_error &e) { \ 84 std::cout << "Runtime exception returned from " << pp << ": " \ 85 << e.what() << " FAIL" << std::endl; \ 88 catch (std::logic_error &e) { \ 89 std::cout << "Logic exception returned from " << pp << ": " \ 90 << e.what() << " FAIL" << std::endl; \ 93 catch (std::bad_alloc &e) { \ 94 std::cout << "Bad_alloc exception returned from " << pp << ": " \ 95 << e.what() << " FAIL" << std::endl; \ 98 catch (std::exception &e) { \ 99 std::cout << "Unknown exception returned from " << pp << ": " \ 100 << e.what() << " FAIL" << std::endl; \ 107 template <
typename tMVector_t>
114 template <
typename tMVector_t>
122 template <
typename tMVector_t>
124 int numObj, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids,
125 int dim,
double *coords_,
int *ierr)
127 typedef typename tMVector_t::scalar_type scalar_t;
133 double *val = coords_;
134 const scalar_t *x = dots_->
coordinates->getData(0).getRawPtr();
135 const scalar_t *y = dots_->
coordinates->getData(1).getRawPtr();
136 const scalar_t *z = dots_->
coordinates->getData(2).getRawPtr();
137 for (
int i=0; i < numObj; i++){
138 *val++ =
static_cast<double>(x[i]);
139 *val++ =
static_cast<double>(y[i]);
140 *val++ =
static_cast<double>(z[i]);
146 double *val = coords_;
147 const scalar_t *x = dots_->
coordinates->getData(0).getRawPtr();
148 const scalar_t *y = dots_->
coordinates->getData(1).getRawPtr();
149 for (
int i=0; i < numObj; i++){
150 *val++ =
static_cast<double>(x[i]);
151 *val++ =
static_cast<double>(y[i]);
156 template <
typename tMVector_t>
167 template <
typename tMVector_t>
169 ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids,
170 int num_wgts,
float *obj_wgts,
int *ierr)
172 typedef typename tMVector_t::global_ordinal_type gno_t;
177 size_t localLen = dots_->
coordinates->getLocalLength();
179 dots_->
coordinates->getMap()->getNodeElementList().getRawPtr();
181 if (
sizeof(ZOLTAN_ID_TYPE) ==
sizeof(gno_t))
182 memcpy(gids, ids,
sizeof(ZOLTAN_ID_TYPE) * localLen);
184 for (
size_t i=0; i < localLen; i++)
185 gids[i] = static_cast<ZOLTAN_ID_TYPE>(ids[i]);
188 float *wgts = obj_wgts;
189 for (
size_t i=0; i < localLen; i++)
190 for (
int w=0; w < num_wgts; w++)
191 *wgts++ = dots_->
weights[w][i];
197 const std::string
shapes[] = {
"SQUARE",
"RECTANGLE",
"CIRCLE",
"CUBE",
"RECTANGULAR_PRISM",
"SPHERE"};
198 #define SHAPE_COUNT 6 202 #define DISTRIBUTION_COUNT 2 204 #define HOLE_ALLOC_STEP 10 205 #define MAX_WEIGHT_DIM 10 206 #define INVALID(STR) "Invalid argument at " + STR 207 #define INVALIDSHAPE(STR, DIM) "Invalid shape name " + STR + " for " + DIM + ".\nValid shapes are \"SQUARE\", \"RECTANGLE\", \"CIRCLE\" for 2D, and \"CUBE\", \"RECTANGULAR_PRISM\", \"SPHERE\" for 3D" 209 #define INVALID_SHAPE_ARG(SHAPE, REQUIRED) "Invalid argument count for shape " + SHAPE + ". Requires " + REQUIRED + " argument(s)." 210 #define MAX_ITER_ALLOWED 500 214 template <
typename T>
247 template <
typename T>
253 this->center.
x = center_.
x;
254 this->center.
y = center_.
y;
255 this->center.
z = center_.
z;
256 this->edge1 = edge1_;
257 this->edge2 = edge2_;
258 this->edge3 = edge3_;
264 template <
typename T>
270 return fabs(dot.
x - this->center.x) < this->edge1 / 2 && fabs(dot.
y - this->center.y) < this->edge1 / 2;
274 template <
typename T>
279 return fabs(dot.
x - this->center.x) < this->edge1 / 2 && fabs(dot.
y - this->center.y) < this->edge2 / 2;
284 template <
typename T>
290 return (dot.
x - this->center.x)*(dot.
x - this->center.x) + (dot.
y - this->center.y) * (dot.
y - this->center.y) < this->edge1 * this->edge1;
294 template <
typename T>
300 return fabs(dot.
x - this->center.x) < this->edge1 / 2 && fabs(dot.
y - this->center.y) < this->edge1 / 2 && fabs(dot.
z - this->center.z) < this->edge1 / 2;
304 template <
typename T>
310 return fabs(dot.
x - this->center.x) < this->edge1 / 2 && fabs(dot.
y - this->center.y) < this->edge2 / 2 && fabs(dot.
z - this->center.z) < this->edge3 / 2;
314 template <
typename T>
320 return fabs((dot.
x - this->center.x) * (dot.
x - this->center.x) * (dot.
x - this->center.x)) +
321 fabs((dot.
y - this->center.y) * (dot.
y - this->center.y) * (dot.
y - this->center.y)) +
322 fabs((dot.
z - this->center.z) * (dot.
z - this->center.z) * (dot.
z - this->center.z))
324 this->edge1 * this->edge1 * this->edge1;
328 template <
typename T,
typename weighttype>
332 virtual weighttype get1DWeight(T x)=0;
333 virtual weighttype get2DWeight(T x, T y)=0;
334 virtual weighttype get3DWeight(T x, T y, T z)=0;
357 template <
typename T,
typename weighttype>
368 SteppedEquation(T a1_, T a2_, T a3_, T b1_, T b2_, T b3_, T c_, T x1_, T y1_, T z1_, T *steps_, T *values_,
int stepCount_):
WeightDistribution<T,weighttype>(){
381 this->stepCount = stepCount_;
382 if(this->stepCount > 0){
383 this->steps =
new T[this->stepCount];
384 this->values =
new T[this->stepCount + 1];
386 for (
int i = 0; i < this->stepCount; ++i){
387 this->steps[i] = steps_[i];
388 this->values[i] = values_[i];
390 this->values[this->stepCount] = values_[this->stepCount];
396 if(this->stepCount > 0){
397 delete [] this->steps;
398 delete [] this->values;
404 T expressionRes = this->a1 * pow( (x - this->x1), b1) + c;
405 if(this->stepCount > 0){
406 for (
int i = 0; i < this->stepCount; ++i){
407 if (expressionRes < this->steps[i])
return this->values[i];
409 return this->values[this->stepCount];
412 return weighttype(expressionRes);
417 T expressionRes = this->a1 * pow( (x - this->x1), b1) + this->a2 * pow( (y - this->y1), b2) + c;
418 if(this->stepCount > 0){
419 for (
int i = 0; i < this->stepCount; ++i){
420 if (expressionRes < this->steps[i])
return this->values[i];
422 return this->values[this->stepCount];
425 return weighttype(expressionRes);
430 std::cout << this->a1 <<
"*" <<
"math.pow( (" << x <<
"- " << this->x1 <<
" ), " << b1 <<
")" <<
"+" << this->a2<<
"*"<<
"math.pow( (" << y <<
"-" << this->y1 <<
"), " <<
431 b2 <<
" ) + " << this->a3 <<
" * math.pow( (" << z <<
"-" << this->z1 <<
"), " << b3 <<
")+ " << c <<
" == " <<
432 this->a1 * pow( (x - this->x1), b1) + this->a2 * pow( (y - this->y1), b2) + this->a3 * pow( (z - this->z1), b3) + c << std::endl;
437 T expressionRes = this->a1 * pow( (x - this->x1), b1) + this->a2 * pow( (y - this->y1), b2) + this->a3 * pow( (z - this->z1), b3) + this->c;
440 if(this->stepCount > 0){
441 for (
int i = 0; i < this->stepCount; ++i){
442 if (expressionRes < this->steps[i]) {
444 return this->values[i];
449 return this->values[this->stepCount];
452 return weighttype(expressionRes);
457 T expressionRes = this->a1 * pow( (p.
x - this->x1), b1) + this->a2 * pow( (p.
y - this->y1), b2) + this->a3 * pow( (p.
z - this->z1), b3);
458 if(this->stepCount > 0){
459 for (
int i = 0; i < this->stepCount; ++i){
460 if (expressionRes < this->steps[i])
return this->values[i];
462 return this->values[this->stepCount];
465 return weighttype(expressionRes);
470 template <
typename T,
typename lno_t,
typename gno_t>
481 numPoints(np_), dimension(dim), requested(0), assignedPrevious(0),
485 virtual T getXCenter() = 0;
486 virtual T getXRadius() =0;
489 Hole<T> **holes, lno_t holeCount,
490 float *sharedRatios_,
int myRank){
492 for (
int i = 0; i < myRank; ++i){
494 this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
495 if (i < this->numPoints % this->worldSize){
496 this->assignedPrevious += 1;
500 this->requested = requestedPointcount;
502 unsigned int slice = UINT_MAX/(this->worldSize);
503 unsigned int stateBegin = myRank * slice;
508 #ifdef HAVE_ZOLTAN2_OMP 514 #ifdef HAVE_ZOLTAN2_OMP 515 me = omp_get_thread_num();
516 tsize = omp_get_num_threads();
518 unsigned int state = stateBegin + me * slice/(tsize);
520 #ifdef HAVE_ZOLTAN2_OMP 523 for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
527 throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
531 bool isInHole =
false;
532 for(lno_t i = 0; i < holeCount; ++i){
533 if(holes[i][0].isInArea(p)){
538 if(isInHole)
continue;
577 void GetPoints(lno_t requestedPointcount, T **coords, lno_t tindex,
578 Hole<T> **holes, lno_t holeCount,
579 float *sharedRatios_,
int myRank){
581 for (
int i = 0; i < myRank; ++i){
583 this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
584 if (gno_t(i) < this->numPoints % this->worldSize){
585 this->assignedPrevious += 1;
589 this->requested = requestedPointcount;
591 unsigned int slice = UINT_MAX/(this->worldSize);
592 unsigned int stateBegin = myRank * slice;
593 #ifdef HAVE_ZOLTAN2_OMP 599 #ifdef HAVE_ZOLTAN2_OMP 600 me = omp_get_thread_num();
601 tsize = omp_get_num_threads();
603 unsigned int state = stateBegin + me * (slice/(tsize));
611 #ifdef HAVE_ZOLTAN2_OMP 614 for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
618 throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
622 bool isInHole =
false;
623 for(lno_t i = 0; i < holeCount; ++i){
624 if(holes[i][0].isInArea(p)){
629 if(isInHole)
continue;
630 coords[0][cnt + tindex] = p.
x;
631 if(this->dimension > 1){
632 coords[1][cnt + tindex] = p.
y;
633 if(this->dimension > 2){
634 coords[2][cnt + tindex] = p.
z;
644 template <
typename T,
typename lno_t,
typename gno_t>
661 T sd_x, T sd_y, T sd_z,
int wSize) :
663 standartDevx(sd_x), standartDevy(sd_y), standartDevz(sd_z)
665 this->center.
x = center_.
x;
666 this->center.
y = center_.
y;
667 this->center.
z = center_.
z;
675 for(
int i = 0; i < this->dimension; ++i){
678 p.
x = normalDist(this->center.
x, this->standartDevx, state);
681 p.
y = normalDist(this->center.
y, this->standartDevy, state);
684 p.
z = normalDist(this->center.
z, this->standartDevz, state);
687 throw "unsupported dimension";
695 T normalDist(T center_, T sd,
unsigned int &state) {
696 static bool derived=
false;
697 static T storedDerivation;
698 T polarsqrt, normalsquared, normal1, normal2;
701 normal1=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
702 normal2=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
703 normalsquared=normal1*normal1+normal2*normal2;
704 }
while ( normalsquared>=1.0 || normalsquared == 0.0);
706 polarsqrt=sqrt(-2.0*log(normalsquared)/normalsquared);
707 storedDerivation=normal1*polarsqrt;
709 return normal2*polarsqrt*sd + center_;
713 return storedDerivation*sd + center_;
718 template <
typename T,
typename lno_t,
typename gno_t>
730 return (rightMostx - leftMostx)/2 + leftMostx;
733 return (rightMostx - leftMostx)/2;
738 T l_z, T r_z,
int wSize ) :
740 leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y),
741 leftMostz(l_z), rightMostz(r_z){}
749 for(
int i = 0; i < this->dimension; ++i){
752 p.
x = uniformDist(this->leftMostx, this->rightMostx, state);
755 p.
y = uniformDist(this->leftMosty, this->rightMosty, state);
758 p.
z = uniformDist(this->leftMostz, this->rightMostz, state);
761 throw "unsupported dimension";
769 T uniformDist(T a, T b,
unsigned int &state)
771 return (b-a)*(rand_r(&state) / double(RAND_MAX)) + a;
775 template <
typename T,
typename lno_t,
typename gno_t>
792 return (rightMostx - leftMostx)/2 + leftMostx;
795 return (rightMostx - leftMostx)/2;
800 T l_x, T r_x, T l_y, T r_y, T l_z, T r_z ,
801 int myRank_,
int wSize) :
803 leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y), leftMostz(l_z), rightMostz(r_z), myRank(myRank_){
805 this->processCnt = 0;
806 this->along_X = alongX; this->along_Y = alongY; this->along_Z = alongZ;
808 if(this->along_X > 1)
809 this->xstep = (rightMostx - leftMostx) / (alongX - 1);
812 if(this->along_Y > 1)
813 this->ystep = (rightMosty - leftMosty) / (alongY - 1);
816 if(this->along_Z > 1)
817 this->zstep = (rightMostz - leftMostz) / (alongZ - 1);
820 xshift = 0; yshift = 0; zshift = 0;
835 this->zshift = pindex / (along_X * along_Y);
836 this->yshift = (pindex % (along_X * along_Y)) / along_X;
837 this->xshift = (pindex % (along_X * along_Y)) % along_X;
844 p.
x = xshift * this->xstep + leftMostx;
845 p.
y = yshift * this->ystep + leftMosty;
846 p.
z = zshift * this->zstep + leftMostz;
901 template <
typename scalar_t,
typename lno_t,
typename gno_t,
typename node_t>
906 int coordinate_dimension;
907 gno_t numGlobalCoords;
908 lno_t numLocalCoords;
909 float *loadDistributions;
911 bool distinctCoordSet;
913 int distributionCount;
918 int numWeightsPerCoord;
920 RCP<const Teuchos::Comm<int> > comm;
932 float perturbation_ratio;
934 typedef Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t>
tMVector_t;
935 typedef Tpetra::Map<lno_t, gno_t, node_t>
tMap_t;
938 template <
typename tt>
939 tt getParamVal(
const Teuchos::ParameterEntry& pe,
const std::string ¶mname){
942 returnVal = Teuchos::getValue<tt>(pe);
950 int countChar (std::string inStr,
char cntChar){
952 for (
unsigned int i = 0; i < inStr.size(); ++i){
953 if (inStr[i] == cntChar) {
960 template <
typename tt>
962 std::stringstream ss (std::stringstream::in |std::stringstream::out);
964 std::string tmp =
"";
969 template <
typename tt>
970 tt fromString(std::string obj){
971 std::stringstream ss (std::stringstream::in | std::stringstream::out);
976 throw "Cannot convert string " + obj;
982 void splitString(std::string inStr,
char splitChar, std::string *outSplittedStr){
983 std::stringstream ss (std::stringstream::in | std::stringstream::out);
987 std::string tmp =
"";
988 std::getline(ss, tmp ,splitChar);
989 outSplittedStr[cnt++] = tmp;
1026 void getCoordinateDistributions(std::string coordinate_distributions){
1027 if(coordinate_distributions ==
""){
1028 throw "At least one distribution function must be provided to coordinate generator.";
1031 int argCnt = this->countChar(coordinate_distributions,
',') + 1;
1032 std::string *splittedStr =
new std::string[argCnt];
1033 splitString(coordinate_distributions,
',', splittedStr);
1035 for(
int i = 0; i < argCnt; ){
1038 std::string distName = splittedStr[i++];
1040 if(distName ==
"NORMAL"){
1042 if (this->coordinate_dimension == 3){
1045 if(i + reqArg > argCnt) {
1046 std::string tmp = toString<int>(reqArg);
1049 np_ = fromString<gno_t>(splittedStr[i++]);
1052 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1053 pp.y = fromString<scalar_t>(splittedStr[i++]);
1055 if(this->coordinate_dimension == 3){
1056 pp.z = fromString<scalar_t>(splittedStr[i++]);
1059 scalar_t sd_x = fromString<scalar_t>(splittedStr[i++]);
1060 scalar_t sd_y = fromString<scalar_t>(splittedStr[i++]);
1062 if(this->coordinate_dimension == 3){
1063 sd_z = fromString<scalar_t>(splittedStr[i++]);
1067 }
else if(distName ==
"UNIFORM" ){
1069 if (this->coordinate_dimension == 3){
1072 if(i + reqArg > argCnt) {
1073 std::string tmp = toString<int>(reqArg);
1076 np_ = fromString<gno_t>(splittedStr[i++]);
1077 scalar_t l_x = fromString<scalar_t>(splittedStr[i++]);
1078 scalar_t r_x = fromString<scalar_t>(splittedStr[i++]);
1079 scalar_t l_y = fromString<scalar_t>(splittedStr[i++]);
1080 scalar_t r_y = fromString<scalar_t>(splittedStr[i++]);
1082 scalar_t l_z = 0, r_z = 0;
1084 if(this->coordinate_dimension == 3){
1085 l_z = fromString<scalar_t>(splittedStr[i++]);
1086 r_z = fromString<scalar_t>(splittedStr[i++]);
1089 this->coordinateDistributions[this->distributionCount++] =
new CoordinateUniformDistribution<scalar_t, lno_t,gno_t>( np_, this->coordinate_dimension, l_x, r_x, l_y, r_y, l_z, r_z, this->worldSize );
1090 }
else if (distName ==
"GRID"){
1092 if(this->coordinate_dimension == 3){
1095 if(i + reqArg > argCnt) {
1096 std::string tmp = toString<int>(reqArg);
1100 gno_t np_x = fromString<gno_t>(splittedStr[i++]);
1101 gno_t np_y = fromString<gno_t>(splittedStr[i++]);
1105 if(this->coordinate_dimension == 3){
1106 np_z = fromString<gno_t>(splittedStr[i++]);
1109 np_ = np_x * np_y * np_z;
1110 scalar_t l_x = fromString<scalar_t>(splittedStr[i++]);
1111 scalar_t r_x = fromString<scalar_t>(splittedStr[i++]);
1112 scalar_t l_y = fromString<scalar_t>(splittedStr[i++]);
1113 scalar_t r_y = fromString<scalar_t>(splittedStr[i++]);
1115 scalar_t l_z = 0, r_z = 0;
1117 if(this->coordinate_dimension == 3){
1118 l_z = fromString<scalar_t>(splittedStr[i++]);
1119 r_z = fromString<scalar_t>(splittedStr[i++]);
1122 if(np_x < 1 || np_z < 1 || np_y < 1 ){
1123 throw "Provide at least 1 point along each dimension for grid test.";
1127 (np_x, np_y,np_z, this->coordinate_dimension, l_x, r_x,l_y, r_y, l_z, r_z , this->myRank, this->worldSize);
1131 std::string tmp = toString<int>(this->coordinate_dimension);
1134 this->numGlobalCoords += (gno_t) np_;
1136 delete []splittedStr;
1139 void getProcLoadDistributions(std::string proc_load_distributions){
1141 this->loadDistributions =
new float[this->worldSize];
1142 if(proc_load_distributions ==
""){
1143 float r = 1.0 / this->worldSize;
1144 for(
int i = 0; i < this->worldSize; ++i){
1145 this->loadDistributions[i] = r;
1150 int argCnt = this->countChar(proc_load_distributions,
',') + 1;
1151 if(argCnt != this->worldSize) {
1152 throw "Invalid parameter count load distributions. Given " + toString<int>(argCnt) +
" processor size is " + toString<int>(this->worldSize);
1154 std::string *splittedStr =
new std::string[argCnt];
1155 splitString(proc_load_distributions,
',', splittedStr);
1156 for(
int i = 0; i < argCnt; ++i){
1157 this->loadDistributions[i] = fromString<float>(splittedStr[i]);
1159 delete []splittedStr;
1163 for(
int i = 0; i < this->worldSize; ++i){
1164 sum += this->loadDistributions[i];
1167 throw "Processor load ratios do not sum to 1.0.";
1173 void getHoles(std::string holeDescription){
1174 if(holeDescription ==
""){
1178 int argCnt = this->countChar(holeDescription,
',') + 1;
1179 std::string *splittedStr =
new std::string[argCnt];
1180 splitString(holeDescription,
',', splittedStr);
1182 for(
int i = 0; i < argCnt; ){
1185 std::string shapeName = splittedStr[i++];
1186 if(shapeName ==
"SQUARE" && this->coordinate_dimension == 2){
1187 if(i + 3 > argCnt) {
1191 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1192 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1193 scalar_t edge = fromString<scalar_t>(splittedStr[i++]);
1195 }
else if(shapeName ==
"RECTANGLE" && this->coordinate_dimension == 2){
1196 if(i + 4 > argCnt) {
1200 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1201 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1202 scalar_t edgex = fromString<scalar_t>(splittedStr[i++]);
1203 scalar_t edgey = fromString<scalar_t>(splittedStr[i++]);
1206 }
else if(shapeName ==
"CIRCLE" && this->coordinate_dimension == 2){
1207 if(i + 3 > argCnt) {
1211 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1212 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1213 scalar_t r = fromString<scalar_t>(splittedStr[i++]);
1215 }
else if(shapeName ==
"CUBE" && this->coordinate_dimension == 3){
1216 if(i + 4 > argCnt) {
1220 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1221 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1222 pp.
z = fromString<scalar_t>(splittedStr[i++]);
1223 scalar_t edge = fromString<scalar_t>(splittedStr[i++]);
1225 }
else if(shapeName ==
"RECTANGULAR_PRISM" && this->coordinate_dimension == 3){
1226 if(i + 6 > argCnt) {
1230 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1231 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1232 pp.
z = fromString<scalar_t>(splittedStr[i++]);
1233 scalar_t edgex = fromString<scalar_t>(splittedStr[i++]);
1234 scalar_t edgey = fromString<scalar_t>(splittedStr[i++]);
1235 scalar_t edgez = fromString<scalar_t>(splittedStr[i++]);
1238 }
else if(shapeName ==
"SPHERE" && this->coordinate_dimension == 3){
1239 if(i + 4 > argCnt) {
1243 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1244 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1245 pp.
z = fromString<scalar_t>(splittedStr[i++]);
1246 scalar_t r = fromString<scalar_t>(splittedStr[i++]);
1249 std::string tmp = toString<int>(this->coordinate_dimension);
1253 delete [] splittedStr;
1256 void getWeightDistribution(std::string *weight_distribution_arr,
int wdimension){
1261 std::string weight_distribution = weight_distribution_arr[ii];
1262 if(weight_distribution ==
"")
continue;
1263 if(wcount == wdimension) {
1264 throw "Weight Dimension is provided as " + toString<int>(wdimension) +
". More weight distribution is provided.";
1267 int count = this->countChar(weight_distribution,
' ');
1268 std::string *splittedStr =
new string[count + 1];
1269 this->splitString(weight_distribution,
' ', splittedStr);
1272 scalar_t a1=0,a2=0,a3=0;
1273 scalar_t x1=0,y1=0,z1=0;
1274 scalar_t b1=1,b2=1,b3=1;
1275 scalar_t *steps = NULL;
1276 scalar_t *values= NULL;
1280 for (
int i = 1; i < count + 1; ++i){
1281 int assignmentCount = this->countChar(splittedStr[i],
'=');
1282 if (assignmentCount != 1){
1283 throw "Error at the argument " + splittedStr[i];
1285 std::string parameter_vs_val[2];
1286 this->splitString(splittedStr[i],
'=', parameter_vs_val);
1287 std::string parameter = parameter_vs_val[0];
1288 std::string value = parameter_vs_val[1];
1291 if (parameter ==
"a1"){
1292 a1 = this->fromString<scalar_t>(value);
1294 else if (parameter ==
"a2"){
1295 if(this->coordinate_dimension > 1){
1296 a2 = this->fromString<scalar_t>(value);
1299 throw parameter+
" argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1303 else if (parameter ==
"a3"){
1304 if(this->coordinate_dimension > 2){
1305 a3 = this->fromString<scalar_t>(value);
1308 throw parameter+
" argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1311 else if (parameter ==
"b1"){
1312 b1 = this->fromString<scalar_t>(value);
1314 else if (parameter ==
"b2"){
1315 if(this->coordinate_dimension > 1){
1316 b2 = this->fromString<scalar_t>(value);
1319 throw parameter+
" argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1322 else if (parameter ==
"b3"){
1324 if(this->coordinate_dimension > 2){
1325 b3 = this->fromString<scalar_t>(value);
1328 throw parameter+
" argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1331 else if (parameter ==
"c"){
1332 c = this->fromString<scalar_t>(value);
1334 else if (parameter ==
"x1"){
1335 x1 = this->fromString<scalar_t>(value);
1337 else if (parameter ==
"y1"){
1338 if(this->coordinate_dimension > 1){
1339 y1 = this->fromString<scalar_t>(value);
1342 throw parameter+
" argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1345 else if (parameter ==
"z1"){
1346 if(this->coordinate_dimension > 2){
1347 z1 = this->fromString<scalar_t>(value);
1350 throw parameter+
" argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1353 else if (parameter ==
"steps"){
1354 stepCount = this->countChar(value,
',') + 1;
1355 std::string *stepstr =
new std::string[stepCount];
1356 this->splitString(value,
',', stepstr);
1357 steps =
new scalar_t[stepCount];
1358 for (
int j = 0; j < stepCount; ++j){
1359 steps[j] = this->fromString<scalar_t>(stepstr[j]);
1363 else if (parameter ==
"values"){
1364 valueCount = this->countChar(value,
',') + 1;
1365 std::string *stepstr =
new std::string[valueCount];
1366 this->splitString(value,
',', stepstr);
1367 values =
new scalar_t[valueCount];
1368 for (
int j = 0; j < valueCount; ++j){
1369 values[j] = this->fromString<scalar_t>(stepstr[j]);
1374 throw "Invalid parameter name at " + splittedStr[i];
1378 delete []splittedStr;
1379 if(stepCount + 1!= valueCount){
1380 throw "Step count: " + this->toString<int>(stepCount) +
" must be 1 less than value count: " + this->toString<int>(valueCount);
1384 this->wd[wcount++] =
new SteppedEquation<scalar_t,scalar_t>(a1, a2, a3, b1, b2, b3, c, x1, y1, z1, steps, values, stepCount);
1392 if(wcount != this->numWeightsPerCoord){
1393 throw "Weight Dimension is provided as " + toString<int>(wdimension) +
". But " + toString<int>(wcount)+
" weight distributions are provided.";
1397 void parseParams(
const Teuchos::ParameterList & params){
1399 std::string holeDescription =
"";
1400 std::string proc_load_distributions =
"";
1401 std::string distinctDescription =
"";
1402 std::string coordinate_distributions =
"";
1405 numWeightsPerCoord_parameters[i] =
"";
1409 for (Teuchos::ParameterList::ConstIterator pit = params.begin(); pit != params.end(); ++pit ){
1410 const std::string &
paramName = params.name(pit);
1412 const Teuchos::ParameterEntry &pe = params.entry(pit);
1414 if(paramName.find(
"hole-") == 0){
1415 if(holeDescription ==
""){
1416 holeDescription = getParamVal<std::string>(pe,
paramName);
1419 holeDescription +=
","+getParamVal<std::string>(pe,
paramName);
1422 else if(paramName.find(
"distribution-") == 0){
1423 if(coordinate_distributions ==
"")
1424 coordinate_distributions = getParamVal<std::string>(pe,
paramName);
1426 coordinate_distributions +=
","+getParamVal<std::string>(pe,
paramName);
1431 else if (paramName.find(weight_distribution_string) == 0){
1432 std::string weight_dist_param = paramName.substr(weight_distribution_string.size());
1433 int dash_pos = weight_dist_param.find(
"-");
1434 std::string distribution_index_string = weight_dist_param.substr(0, dash_pos);
1435 int distribution_index = fromString<int>(distribution_index_string);
1437 if(distribution_index >= MAX_WEIGHT_DIM){
1438 throw "Given distribution index:" + distribution_index_string +
" larger than maximum allowed number of weights:" + toString<int>(
MAX_WEIGHT_DIM);
1440 numWeightsPerCoord_parameters[distribution_index] +=
" " + weight_dist_param.substr(dash_pos + 1)+
"="+ getParamVal<std::string>(pe,
paramName);
1442 else if(paramName ==
"dim"){
1443 int dim = fromString<int>(getParamVal<std::string>(pe,
paramName));
1444 if(dim < 2 && dim > 3){
1447 this->coordinate_dimension = dim;
1450 else if(paramName ==
"wdim"){
1451 int dim = fromString<int>(getParamVal<std::string>(pe,
paramName));
1452 if(dim < 1 && dim > MAX_WEIGHT_DIM){
1455 this->numWeightsPerCoord = dim;
1458 else if(paramName ==
"predistribution"){
1459 int pre_distribution = fromString<int>(getParamVal<std::string>(pe,
paramName));
1460 if(pre_distribution < 0 && pre_distribution > 3){
1463 this->predistribution = pre_distribution;
1466 else if(paramName ==
"perturbation_ratio"){
1467 float perturbation = fromString<float>(getParamVal<std::string>(pe,
paramName));
1468 if(perturbation < 0 && perturbation > 1){
1471 this->perturbation_ratio = perturbation;
1476 else if(paramName ==
"proc_load_distributions"){
1477 proc_load_distributions = getParamVal<std::string>(pe,
paramName);
1478 this->loadDistSet =
true;
1481 else if(paramName ==
"distinct_coordinates"){
1482 distinctDescription = getParamVal<std::string>(pe,
paramName);
1483 if(distinctDescription ==
"T"){
1484 this->distinctCoordSet =
true;
1485 }
else if(distinctDescription ==
"F"){
1486 this->distinctCoordSet =
true;
1488 throw "Invalid parameter for distinct_coordinates: " + distinctDescription +
". Candidates: T or F." ;
1492 else if(paramName ==
"out_file"){
1493 this->outfile = getParamVal<std::string>(pe,
paramName);
1502 if(this->coordinate_dimension == 0){
1503 throw "Dimension must be provided to coordinate generator.";
1522 if (this->loadDistSet && this->distinctCoordSet){
1523 throw "distinct_coordinates and proc_load_distributions parameters cannot be satisfied together.";
1525 this->getHoles(holeDescription);
1527 this->getProcLoadDistributions(proc_load_distributions);
1528 this->getCoordinateDistributions(coordinate_distributions);
1529 this->getWeightDistribution(numWeightsPerCoord_parameters, this->numWeightsPerCoord);
1536 catch(std::string s){
1544 catch(
char const* s){
1553 for (
int i = 0; i < this->holeCount; ++i){
1554 delete this->holes[i];
1560 delete []loadDistributions;
1562 if(coordinateDistributions){
1564 for (
int i = 0; i < this->distributionCount; ++i){
1565 delete this->coordinateDistributions[i];
1567 free (this->coordinateDistributions);
1570 for (
int i = 0; i < this->numWeightsPerCoord; ++i){
1576 if(this->numWeightsPerCoord){
1577 for(
int i = 0; i < this->numWeightsPerCoord; ++i)
1578 delete [] this->wghts[i];
1579 delete []this->wghts;
1581 if(this->coordinate_dimension){
1582 for(
int i = 0; i < this->coordinate_dimension; ++i)
1583 delete [] this->coords[i];
1584 delete [] this->coords;
1590 std::cout <<
"\nGeometric Generator Parameter File Format:" << std::endl;
1591 std::cout <<
"- dim=Coordinate Dimension: 2 or 3" << std::endl;
1592 std::cout <<
"- Available distributions:" << std::endl;
1593 std::cout <<
"\tUNIFORM: -> distribution-1=UNIFORM,NUMCOORDINATES,XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << std::endl;
1594 std::cout <<
"\tGRID: -> distribution-2=GRID,XLENGTH,YLENGTH{,ZLENGTH},XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << std::endl;
1595 std::cout <<
"\tNORMAL: -> distribution-3=NORMAL,XCENTER,YCENTER{,ZCENTER},XSD,YSD,{,ZSD}" << std::endl;
1596 std::cout <<
"- wdim=numWeightsPerCoord: There should be as many weight function as number of weights per coord." << std::endl;
1597 std::cout <<
"- Weight Equation: w = (a1 * (x - x1)^b1) + (a2 * (y - y1)^b2) + (a3 * (z - z1)^b3) + c" << std::endl;
1598 std::cout <<
"Parameter settings:" << std::endl;
1599 std::cout <<
"\tWeightDistribution-1-a1=a1 " << std::endl;
1600 std::cout <<
"\tWeightDistribution-1-a2=a2 " << std::endl;
1601 std::cout <<
"\tWeightDistribution-1-a3=a3 " << std::endl;
1602 std::cout <<
"\tWeightDistribution-1-b1=b1 " << std::endl;
1603 std::cout <<
"\tWeightDistribution-1-b2=b2 " << std::endl;
1604 std::cout <<
"\tWeightDistribution-1-b3=b3 " << std::endl;
1605 std::cout <<
"\tWeightDistribution-1-x1=x1 " << std::endl;
1606 std::cout <<
"\tWeightDistribution-1-y1=y1 " << std::endl;
1607 std::cout <<
"\tWeightDistribution-1-z1=z1 " << std::endl;
1608 std::cout <<
"\tWeightDistribution-1-c=c " << std::endl;
1609 std::cout <<
"\tIt is possible to set step function to the result of weight equation." << std::endl;
1610 std::cout <<
"\tWeightDistribution-1-steps=step1,step2,step3:increasing order" << std::endl;
1611 std::cout <<
"\tWeightDistribution-1-values=val1,val2,val3,val4." << std::endl;
1612 std::cout <<
"\t\tIf w < step1 -> w = val1" << std::endl;
1613 std::cout <<
"\t\tElse if w < step2 -> w = val2" << std::endl;
1614 std::cout <<
"\t\tElse if w < step3 -> w = val3" << std::endl;
1615 std::cout <<
"\t\tElse -> w = val4" << std::endl;
1616 std::cout <<
"- Holes:" << std::endl;
1617 std::cout <<
"\thole-1:SPHERE,XCENTER,YCENTER,ZCENTER,RADIUS (only for dim=3)" << std::endl;
1618 std::cout <<
"\thole-2:CUBE,XCENTER,YCENTER,ZCENTER,EDGE (only for dim=3)" << std::endl;
1619 std::cout <<
"\thole-3:RECTANGULAR_PRISM,XCENTER,YCENTER,ZCENTER,XEDGE,YEDGE,ZEDGE (only for dim=3)" << std::endl;
1620 std::cout <<
"\thole-4:SQUARE,XCENTER,YCENTER,EDGE (only for dim=2)" << std::endl;
1621 std::cout <<
"\thole-5:RECTANGLE,XCENTER,YCENTER,XEDGE,YEDGE (only for dim=2)" << std::endl;
1622 std::cout <<
"\thole-6:CIRCLE,XCENTER,YCENTER,RADIUS (only for dim=2)" << std::endl;
1623 std::cout <<
"- out_file:out_file_path : if provided output will be written to files." << std::endl;
1624 std::cout <<
"- proc_load_distributions:ratio_0, ratio_1, ratio_2....ratio_n. Loads of each processor, should be as many as MPI ranks and should sum up to 1." << std::endl;
1625 std::cout <<
"- predistribution:distribution_option. Predistribution of the coordinates to the processors. 0 for NONE, 1 RCB, 2 MJ, 3 BLOCK." << std::endl;
1626 std::cout <<
"- perturbation_ratio:the percent of the local data which will be perturbed in order to simulate the changes in the dynamic partitioning. Float value between 0 and 1." << std::endl;
1633 this->coordinate_dimension = 0;
1634 this->numWeightsPerCoord = 0;
1635 this->worldSize = comm_->getSize();
1636 this->numGlobalCoords = 0;
1637 this->loadDistributions = NULL;
1639 this->coordinateDistributions = NULL;
1640 this->holeCount = 0;
1641 this->distributionCount = 0;
1643 this->predistribution = 0;
1644 this->perturbation_ratio = 0;
1655 this->loadDistSet =
false;
1656 this->distinctCoordSet =
false;
1657 this->myRank = comm_->getRank();
1660 this->parseParams(params);
1662 catch(std::string s){
1664 print_description();
1672 lno_t myPointCount = 0;
1673 this->numGlobalCoords = 0;
1675 gno_t prefixSum = 0;
1676 for(
int i = 0; i < this->distributionCount; ++i){
1677 for(
int ii = 0; ii < this->worldSize; ++ii){
1678 lno_t increment = lno_t (this->coordinateDistributions[i]->numPoints * this->loadDistributions[ii]);
1679 if (gno_t(ii) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1682 this->numGlobalCoords += increment;
1684 prefixSum += increment;
1687 myPointCount += lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
1688 if (gno_t(myRank) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1693 this->coords =
new scalar_t *[this->coordinate_dimension];
1694 for(
int i = 0; i < this->coordinate_dimension; ++i){
1695 this->coords[i] =
new scalar_t[myPointCount];
1698 for (
int ii = 0; ii < this->coordinate_dimension; ++ii){
1699 #ifdef HAVE_ZOLTAN2_OMP 1700 #pragma omp parallel for 1702 for(lno_t i = 0; i < myPointCount; ++i){
1703 this->coords[ii][i] = 0;
1707 this->numLocalCoords = 0;
1708 srand ((myRank + 1) * this->numLocalCoords);
1709 for (
int i = 0; i < distributionCount; ++i){
1711 lno_t requestedPointCount = lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
1712 if (gno_t(myRank) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1713 requestedPointCount += 1;
1717 this->coordinateDistributions[i]->
GetPoints(requestedPointCount,this->coords, this->numLocalCoords, this->holes, this->holeCount, this->loadDistributions, myRank);
1718 this->numLocalCoords += requestedPointCount;
1740 if (this->predistribution){
1747 if (this->perturbation_ratio > 0){
1748 this->perturb_data(this->perturbation_ratio, scale);
1769 if (this->distinctCoordSet){
1773 if(this->outfile !=
""){
1775 std::ofstream myfile;
1776 myfile.open ((this->outfile + toString<int>(myRank)).c_str());
1777 for(lno_t i = 0; i < this->numLocalCoords; ++i){
1779 myfile << this->coords[0][i];
1780 if(this->coordinate_dimension > 1){
1781 myfile <<
" " << this->coords[1][i];
1783 if(this->coordinate_dimension > 2){
1784 myfile <<
" " << this->coords[2][i];
1786 myfile << std::endl;
1790 if (this->myRank == 0){
1791 std::ofstream gnuplotfile(
"gnu.gnuplot");
1792 for(
int i = 0; i < this->worldSize; ++i){
1794 if (this->coordinate_dimension == 2){
1800 gnuplotfile << s <<
" \"" << (this->outfile + toString<int>(i)) <<
"\"" << std::endl;
1802 gnuplotfile <<
"pause -1" << std::endl;
1803 gnuplotfile.close();
1816 if (this->numWeightsPerCoord > 0){
1817 this->wghts =
new scalar_t *[this->numWeightsPerCoord];
1818 for(
int i = 0; i < this->numWeightsPerCoord; ++i){
1819 this->wghts[i] =
new scalar_t[this->numLocalCoords];
1823 for(
int ii = 0; ii < this->numWeightsPerCoord; ++ii){
1824 switch(this->coordinate_dimension){
1826 #ifdef HAVE_ZOLTAN2_OMP 1827 #pragma omp parallel for 1829 for (lno_t i = 0; i < this->numLocalCoords; ++i){
1830 this->wghts[ii][i] = this->wd[ii]->
get1DWeight(this->coords[0][i]);
1834 #ifdef HAVE_ZOLTAN2_OMP 1835 #pragma omp parallel for 1837 for (lno_t i = 0; i < this->numLocalCoords; ++i){
1838 this->wghts[ii][i] = this->wd[ii]->
get2DWeight(this->coords[0][i], this->coords[1][i]);
1842 #ifdef HAVE_ZOLTAN2_OMP 1843 #pragma omp parallel for 1845 for (lno_t i = 0; i < this->numLocalCoords; ++i){
1846 this->wghts[ii][i] = this->wd[ii]->
get3DWeight(this->coords[0][i], this->coords[1][i], this->coords[2][i]);
1857 scalar_t *maxCoords=
new scalar_t [this->coordinate_dimension];
1858 scalar_t *minCoords=
new scalar_t [this->coordinate_dimension];
1859 for (
int dim = 0; dim < this->coordinate_dimension; ++dim){
1860 minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
1861 for (lno_t i = 1; i < this->numLocalCoords; ++i){
1862 if (minCoords[dim] > this->coords[dim][i]){
1863 minCoords[dim] = this->coords[dim][i];
1866 if (maxCoords[dim] < this->coords[dim][i]){
1867 maxCoords[dim] = this->coords[dim][i];
1874 scalar_t center = (maxCoords[dim] + minCoords[dim]) / 2;
1876 minCoords[dim] = center - (center - minCoords[dim]) * scale;
1877 maxCoords[dim] = (maxCoords[dim] - center) * scale + center;
1881 gno_t numLocalToPerturb = gno_t(this->numLocalCoords*used_perturbation_ratio);
1883 for (
int dim = 0; dim < this->coordinate_dimension; ++dim){
1884 scalar_t range = maxCoords[dim] - minCoords[dim];
1885 for (gno_t i = 0; i < numLocalToPerturb; ++i){
1886 this->coords[dim][i] = (rand() / double (RAND_MAX)) * (range) + minCoords[dim];
1899 void getBestSurface (
int remaining,
int *dimProcs,
int dim,
int currentDim,
int &bestSurface,
int *bestDimProcs){
1901 if (currentDim < dim - 1){
1903 while(ipx <= remaining) {
1904 if(remaining % ipx == 0) {
1905 int nremain = remaining / ipx;
1906 dimProcs[currentDim] = ipx;
1907 getBestSurface (nremain, dimProcs, dim, currentDim + 1, bestSurface, bestDimProcs);
1913 dimProcs[currentDim] = remaining;
1915 for (
int i = 0; i < dim; ++i) surface += dimProcs[i];
1916 if (surface < bestSurface){
1917 bestSurface = surface;
1918 for (
int i = 0; i < dim; ++i) bestDimProcs[i] = dimProcs[i];
1926 scalar_t *maxCoords=
new scalar_t [this->coordinate_dimension];
1927 scalar_t *minCoords=
new scalar_t [this->coordinate_dimension];
1928 for (
int dim = 0; dim < this->coordinate_dimension; ++dim){
1929 minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
1930 for (lno_t i = 1; i < this->numLocalCoords; ++i){
1931 if (minCoords[dim] > this->coords[dim][i]){
1932 minCoords[dim] = this->coords[dim][i];
1935 if (maxCoords[dim] < this->coords[dim][i]){
1936 maxCoords[dim] = this->coords[dim][i];
1941 reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MAX,
1942 this->coordinate_dimension,
1947 reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MIN,
1948 this->coordinate_dimension,
1962 typedef SSIZE_T ssize_t;
1969 scalar_t *maxCoords=
new scalar_t [this->coordinate_dimension];
1970 scalar_t *minCoords=
new scalar_t [this->coordinate_dimension];
1972 this->getMinMaxCoords(maxCoords, minCoords);
1978 int remaining = this->worldSize;
1979 int coord_dim = this->coordinate_dimension;
1980 int *dimProcs =
new int[coord_dim];
1981 int *bestDimProcs =
new int[coord_dim];
1984 int bestSurface = 0;
1985 dimProcs[0] = remaining;
1986 for (
int i = 1; i < coord_dim; ++i) dimProcs[i] = 1;
1987 for (
int i = 0; i < coord_dim; ++i) bestSurface += dimProcs[i];
1988 for (
int i = 0; i < coord_dim; ++i) bestDimProcs[i] = dimProcs[i];
1990 getBestSurface ( remaining, dimProcs, coord_dim, currentDim, bestSurface, bestDimProcs);
1998 int *shiftProcCount =
new int[coord_dim];
2002 int remainingProc = this->worldSize;
2003 for (
int dim = 0; dim < coord_dim; ++dim){
2004 remainingProc = remainingProc / bestDimProcs[dim];
2005 shiftProcCount[dim] = remainingProc;
2008 scalar_t *dim_slices =
new scalar_t[coord_dim];
2009 for (
int dim = 0; dim < coord_dim; ++dim){
2010 scalar_t dim_range = maxCoords[dim] - minCoords[dim];
2011 scalar_t slice = dim_range / bestDimProcs[dim];
2012 dim_slices[dim] = slice;
2019 gno_t *numPointsInParts =
new gno_t[this->worldSize];
2020 gno_t *numGlobalPointsInParts =
new gno_t[this->worldSize];
2021 gno_t *numPointsInPartsInclusiveUptoMyIndex =
new gno_t[this->worldSize];
2023 gno_t *partBegins =
new gno_t [this->worldSize];
2024 gno_t *partNext =
new gno_t[this->numLocalCoords];
2027 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2030 for (
int i = 0; i < this->worldSize; ++i) {
2034 for (
int i = 0; i < this->worldSize; ++i)
2035 numPointsInParts[i] = 0;
2037 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2039 for (
int dim = 0; dim < coord_dim; ++dim){
2040 int shift = int ((this->coords[dim][i] - minCoords[dim]) / dim_slices[dim]);
2041 if (shift >= bestDimProcs[dim]){
2042 shift = bestDimProcs[dim] - 1;
2044 shift = shift * shiftProcCount[dim];
2047 numPointsInParts[partIndex] += 1;
2048 coordinate_grid_parts[i] = partIndex;
2050 partNext[i] = partBegins[partIndex];
2051 partBegins[partIndex] = i;
2058 reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2061 numGlobalPointsInParts);
2064 Teuchos::scan<int,gno_t>(
2065 *(this->comm), Teuchos::REDUCE_SUM,
2068 numPointsInPartsInclusiveUptoMyIndex
2087 gno_t optimal_num = gno_t(this->numGlobalCoords/
double(this->worldSize)+0.5);
2089 if (this->myRank == 0){
2090 gno_t totalSize = 0;
2091 for (
int i = 0; i < this->worldSize; ++i){
2092 std::cout <<
"me:" << this->myRank <<
" NumPoints in part:" << i <<
" is: " << numGlobalPointsInParts[i] << std::endl;
2093 totalSize += numGlobalPointsInParts[i];
2095 std::cout <<
"Total:" << totalSize <<
" ng:" << this->numGlobalCoords << std::endl;
2096 std::cout <<
"optimal_num:" << optimal_num << std::endl;
2099 ssize_t *extraInPart =
new ssize_t [this->worldSize];
2101 ssize_t extraExchange = 0;
2102 for (
int i = 0; i < this->worldSize; ++i){
2103 extraInPart[i] = numGlobalPointsInParts[i] - optimal_num;
2104 extraExchange += extraInPart[i];
2106 if (extraExchange != 0){
2108 if (extraExchange < 0) addition = 1;
2109 for (ssize_t i = 0; i < extraExchange; ++i){
2110 extraInPart[i % this->worldSize] += addition;
2118 int overloadedPartCount = 0;
2119 int *overloadedPartIndices =
new int [this->worldSize];
2122 int underloadedPartCount = 0;
2123 int *underloadedPartIndices =
new int [this->worldSize];
2125 for (
int i = 0; i < this->worldSize; ++i){
2126 if(extraInPart[i] > 0){
2127 overloadedPartIndices[overloadedPartCount++] = i;
2129 else if(extraInPart[i] < 0){
2130 underloadedPartIndices[underloadedPartCount++] = i;
2134 int underloadpartindex = underloadedPartCount - 1;
2137 int numPartsISendFrom = 0;
2138 int *mySendFromParts =
new int[this->worldSize * 2];
2139 gno_t *mySendFromPartsCounts =
new gno_t[this->worldSize * 2];
2141 int numPartsISendTo = 0;
2142 int *mySendParts =
new int[this->worldSize * 2];
2143 gno_t *mySendCountsToParts =
new gno_t[this->worldSize * 2];
2152 for (
int i = overloadedPartCount - 1; i >= 0; --i){
2155 int overloadPart = overloadedPartIndices[i];
2156 gno_t overload = extraInPart[overloadPart];
2157 gno_t myload = numPointsInParts[overloadPart];
2161 gno_t inclusiveLoadUpToMe = numPointsInPartsInclusiveUptoMyIndex[overloadPart];
2164 gno_t exclusiveLoadUptoMe = inclusiveLoadUpToMe - myload;
2167 if (exclusiveLoadUptoMe >= overload){
2171 overloadedPartIndices[i] = -1;
2172 extraInPart[overloadPart] = 0;
2174 while (overload > 0){
2175 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2176 ssize_t underload = extraInPart[nextUnderLoadedPart];
2177 ssize_t left = overload + underload;
2180 extraInPart[nextUnderLoadedPart] = 0;
2181 underloadedPartIndices[underloadpartindex--] = -1;
2184 extraInPart[nextUnderLoadedPart] = left;
2189 else if (exclusiveLoadUptoMe < overload){
2192 gno_t mySendCount = overload - exclusiveLoadUptoMe;
2194 gno_t sendAfterMe = 0;
2197 if (mySendCount > myload){
2198 sendAfterMe = mySendCount - myload;
2199 mySendCount = myload;
2205 mySendFromParts[numPartsISendFrom] = overloadPart;
2206 mySendFromPartsCounts[numPartsISendFrom++] = mySendCount;
2209 while (exclusiveLoadUptoMe > 0){
2210 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2211 ssize_t underload = extraInPart[nextUnderLoadedPart];
2212 ssize_t left = exclusiveLoadUptoMe + underload;
2215 extraInPart[nextUnderLoadedPart] = 0;
2216 underloadedPartIndices[underloadpartindex--] = -1;
2219 extraInPart[nextUnderLoadedPart] = left;
2221 exclusiveLoadUptoMe = left;
2225 while (mySendCount > 0){
2226 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2227 ssize_t underload = extraInPart[nextUnderLoadedPart];
2228 ssize_t left = mySendCount + underload;
2231 mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2232 mySendCountsToParts[numPartsISendTo++] = -underload;
2234 extraInPart[nextUnderLoadedPart] = 0;
2235 underloadedPartIndices[underloadpartindex--] = -1;
2238 extraInPart[nextUnderLoadedPart] = left;
2240 mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2241 mySendCountsToParts[numPartsISendTo++] = mySendCount;
2247 while (sendAfterMe > 0){
2248 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2249 ssize_t underload = extraInPart[nextUnderLoadedPart];
2250 ssize_t left = sendAfterMe + underload;
2253 extraInPart[nextUnderLoadedPart] = 0;
2254 underloadedPartIndices[underloadpartindex--] = -1;
2257 extraInPart[nextUnderLoadedPart] = left;
2268 for (
int i = 0 ; i < numPartsISendFrom; ++i){
2271 int sendFromPart = mySendFromParts[i];
2272 ssize_t sendCount = mySendFromPartsCounts[i];
2273 while(sendCount > 0){
2274 int partToSendIndex = numPartsISendTo - 1;
2275 int partToSend = mySendParts[partToSendIndex];
2277 int sendCountToThisPart = mySendCountsToParts[partToSendIndex];
2281 if (sendCountToThisPart <= sendCount){
2282 mySendParts[partToSendIndex] = 0;
2283 mySendCountsToParts[partToSendIndex] = 0;
2285 sendCount -= sendCountToThisPart;
2288 mySendCountsToParts[partToSendIndex] = sendCountToThisPart - sendCount;
2289 sendCountToThisPart = sendCount;
2294 gno_t toChange = partBegins[sendFromPart];
2295 gno_t previous_begin = partBegins[partToSend];
2298 for (
int k = 0; k < sendCountToThisPart - 1; ++k){
2299 coordinate_grid_parts[toChange] = partToSend;
2300 toChange = partNext[toChange];
2302 coordinate_grid_parts[toChange] = partToSend;
2304 gno_t newBegin = partNext[toChange];
2305 partNext[toChange] = previous_begin;
2306 partBegins[partToSend] = partBegins[sendFromPart];
2307 partBegins[sendFromPart] = newBegin;
2316 for (
int i = 0; i < this->worldSize; ++i) numPointsInParts[i] = 0;
2318 for (
int i = 0; i < this->numLocalCoords; ++i){
2319 numPointsInParts[coordinate_grid_parts[i]] += 1;
2322 reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2325 numGlobalPointsInParts);
2326 if (this->myRank == 0){
2327 std::cout <<
"reassigning" << std::endl;
2328 gno_t totalSize = 0;
2329 for (
int i = 0; i < this->worldSize; ++i){
2330 std::cout <<
"me:" << this->myRank <<
" NumPoints in part:" << i <<
" is: " << numGlobalPointsInParts[i] << std::endl;
2331 totalSize += numGlobalPointsInParts[i];
2333 std::cout <<
"Total:" << totalSize <<
" ng:" << this->numGlobalCoords << std::endl;
2336 delete []mySendCountsToParts;
2337 delete []mySendParts;
2338 delete []mySendFromPartsCounts;
2339 delete []mySendFromParts;
2340 delete []underloadedPartIndices;
2341 delete []overloadedPartIndices;
2342 delete []extraInPart;
2344 delete []partBegins;
2345 delete []numPointsInPartsInclusiveUptoMyIndex;
2346 delete []numPointsInParts;
2347 delete []numGlobalPointsInParts;
2349 delete []shiftProcCount;
2350 delete []bestDimProcs;
2351 delete []dim_slices;
2360 Tpetra::Distributor distributor(comm);
2361 ArrayView<const int> pIds( coordinate_grid_parts, this->numLocalCoords);
2367 gno_t numMyNewGnos = distributor.createFromSends(pIds);
2371 this->numLocalCoords = numMyNewGnos;
2374 ArrayRCP<scalar_t> recvBuf2(distributor.getTotalReceiveLength());
2376 for (
int i = 0; i < this->coordinate_dimension; ++i){
2377 ArrayView<scalar_t> s(this->coords[i], this->numLocalCoords);
2378 distributor.doPostsAndWaits<scalar_t>(s, 1, recvBuf2());
2379 delete [] this->coords[i];
2380 this->coords[i] =
new scalar_t[this->numLocalCoords];
2381 for (lno_t j = 0; j < this->numLocalCoords; ++j){
2382 this->coords[i][j] = recvBuf2[j];
2390 int coord_dim = this->coordinate_dimension;
2392 lno_t numLocalPoints = this->numLocalCoords;
2393 gno_t numGlobalPoints = this->numGlobalCoords;
2398 RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
2399 new Tpetra::Map<lno_t, gno_t, node_t> (numGlobalPoints, numLocalPoints, 0, comm));
2401 Teuchos::Array<Teuchos::ArrayView<const scalar_t> > coordView(coord_dim);
2405 for (
int i=0; i < coord_dim; i++){
2406 if(numLocalPoints > 0){
2407 Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalPoints);
2410 Teuchos::ArrayView<const scalar_t> a;
2415 RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmVector = RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >(
2416 new Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t>( mp, coordView.view(0, coord_dim), coord_dim));
2419 RCP<const tMVector_t> coordsConst = Teuchos::rcp_const_cast<
const tMVector_t>(tmVector);
2420 std::vector<const scalar_t *>
weights;
2421 std::vector <int> stride;
2425 inputAdapter_t ia(coordsConst,weights, stride);
2427 Teuchos::RCP <Teuchos::ParameterList> params ;
2428 params =RCP <Teuchos::ParameterList> (
new Teuchos::ParameterList,
true);
2431 params->set(
"algorithm",
"multijagged");
2432 params->set(
"num_global_parts", this->worldSize);
2444 #ifdef HAVE_ZOLTAN2_MPI 2458 const typename inputAdapter_t::part_t *partIds = problem->
getSolution().getPartListView();
2460 for (lno_t i = 0; i < this->numLocalCoords;++i){
2461 coordinate_grid_parts[i] = partIds[i];
2470 int rank = this->myRank;
2471 int nprocs = this->worldSize;
2474 MEMORY_CHECK(rank==0 || rank==nprocs-1,
"After initializing MPI");
2479 string memoryOn(
"memoryOn");
2480 string memoryOff(
"memoryOff");
2481 bool doMemory=
false;
2482 int numGlobalParts = nprocs;
2486 string balanceCount(
"balance_object_count");
2487 string balanceWeight(
"balance_object_weight");
2488 string mcnorm1(
"multicriteria_minimize_total_weight");
2489 string mcnorm2(
"multicriteria_balance_total_maximum");
2490 string mcnorm3(
"multicriteria_minimize_maximum_weight");
2491 string objective(balanceWeight);
2494 CommandLineProcessor commandLine(
false,
true);
2497 int input_option = 0;
2498 commandLine.setOption(
"input_option", &input_option,
2499 "whether to use mesh creation, geometric generator, or file input");
2500 string inputFile =
"";
2502 commandLine.setOption(
"input_file", &inputFile,
2503 "the input file for geometric generator or file input");
2506 commandLine.setOption(
"size", &numGlobalCoords,
2507 "Approximate number of global coordinates.");
2508 commandLine.setOption(
"numParts", &numGlobalParts,
2509 "Number of parts (default is one per proc).");
2510 commandLine.setOption(
"nWeights", &nWeights,
2511 "Number of weights per coordinate, zero implies uniform weights.");
2512 commandLine.setOption(
"debug", &debugLevel,
"Zoltan1 debug level");
2513 commandLine.setOption(
"remap",
"no-remap", &remap,
2514 "Zoltan1 REMAP parameter; disabled by default for scalability testing");
2515 commandLine.setOption(
"timers", &dummyTimer,
"ignored");
2516 commandLine.setOption(memoryOn.c_str(), memoryOff.c_str(), &doMemory,
2517 "do memory profiling");
2519 string doc(balanceCount);
2520 doc.append(
": ignore weights\n");
2521 doc.append(balanceWeight);
2522 doc.append(
": balance on first weight\n");
2523 doc.append(mcnorm1);
2524 doc.append(
": given multiple weights, balance their total.\n");
2525 doc.append(mcnorm3);
2526 doc.append(
": given multiple weights, " 2527 "balance the maximum for each coordinate.\n");
2528 doc.append(mcnorm2);
2529 doc.append(
": given multiple weights, balance the L2 norm of the weights.\n");
2530 commandLine.setOption(
"objective", &objective, doc.c_str());
2532 CommandLineProcessor::EParseCommandLineReturn rc =
2533 commandLine.parse(0, NULL);
2537 if (rc != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
2538 if (rc == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) {
2539 if (rank==0) std::cout <<
"PASS" << std::endl;
2543 if (rank==0) std::cout <<
"FAIL" << std::endl;
2552 int coord_dim = this->coordinate_dimension;
2555 RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
2556 new Tpetra::Map<lno_t, gno_t, node_t> (this->numGlobalCoords, this->numLocalCoords, 0, this->comm));
2558 Teuchos::Array<Teuchos::ArrayView<const scalar_t> > coordView(coord_dim);
2559 for (
int i=0; i < coord_dim; i++){
2560 if(numLocalCoords > 0){
2561 Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalCoords);
2564 Teuchos::ArrayView<const scalar_t> a;
2569 tMVector_t *tmVector =
new tMVector_t( mp, coordView.view(0, coord_dim), coord_dim);
2572 dots_.
weights.resize(nWeights);
2575 MEMORY_CHECK(doMemory && rank==0,
"After creating input");
2580 int aok = Zoltan_Initialize(0,NULL, &ver);
2583 printf(
"Zoltan_Initialize failed\n");
2587 struct Zoltan_Struct *zz;
2588 zz = Zoltan_Create(MPI_COMM_WORLD);
2590 Zoltan_Set_Param(zz,
"LB_METHOD",
"RCB");
2591 Zoltan_Set_Param(zz,
"LB_APPROACH",
"PARTITION");
2592 Zoltan_Set_Param(zz,
"CHECK_GEOM",
"0");
2593 Zoltan_Set_Param(zz,
"NUM_GID_ENTRIES",
"1");
2594 Zoltan_Set_Param(zz,
"NUM_LID_ENTRIES",
"0");
2595 Zoltan_Set_Param(zz,
"RETURN_LISTS",
"PART");
2596 std::ostringstream oss;
2597 oss << numGlobalParts;
2598 Zoltan_Set_Param(zz,
"NUM_GLOBAL_PARTS", oss.str().c_str());
2601 Zoltan_Set_Param(zz,
"DEBUG_LEVEL", oss.str().c_str());
2604 Zoltan_Set_Param(zz,
"REMAP",
"1");
2606 Zoltan_Set_Param(zz,
"REMAP",
"0");
2608 if (objective != balanceCount){
2611 Zoltan_Set_Param(zz,
"OBJ_WEIGHT_DIM", oss.str().c_str());
2613 if (objective == mcnorm1)
2614 Zoltan_Set_Param(zz,
"RCB_MULTICRITERIA_NORM",
"1");
2615 else if (objective == mcnorm2)
2616 Zoltan_Set_Param(zz,
"RCB_MULTICRITERIA_NORM",
"2");
2617 else if (objective == mcnorm3)
2618 Zoltan_Set_Param(zz,
"RCB_MULTICRITERIA_NORM",
"3");
2621 Zoltan_Set_Param(zz,
"OBJ_WEIGHT_DIM",
"0");
2624 Zoltan_Set_Num_Obj_Fn(zz, getNumObj<tMVector_t>, &dots_);
2625 Zoltan_Set_Obj_List_Fn(zz, getObjList<tMVector_t>, &dots_);
2626 Zoltan_Set_Num_Geom_Fn(zz, getDim<tMVector_t>, &dots_);
2627 Zoltan_Set_Geom_Multi_Fn(zz, getCoords<tMVector_t>, &dots_);
2629 int changes, numGidEntries, numLidEntries, numImport, numExport;
2630 ZOLTAN_ID_PTR importGlobalGids, importLocalGids;
2631 ZOLTAN_ID_PTR exportGlobalGids, exportLocalGids;
2632 int *importProcs, *importToPart, *exportProcs, *exportToPart;
2634 MEMORY_CHECK(doMemory && rank==0,
"Before Zoltan_LB_Partition");
2637 aok = Zoltan_LB_Partition(zz, &changes, &numGidEntries, &numLidEntries,
2638 &numImport, &importGlobalGids, &importLocalGids,
2639 &importProcs, &importToPart,
2640 &numExport, &exportGlobalGids, &exportLocalGids,
2641 &exportProcs, &exportToPart);
2644 MEMORY_CHECK(doMemory && rank==0,
"After Zoltan_LB_Partition");
2646 for (lno_t i = 0; i < numLocalCoords; i++)
2647 coordinate_grid_parts[i] = exportToPart[i];
2648 Zoltan_Destroy(&zz);
2649 MEMORY_CHECK(doMemory && rank==0,
"After Zoltan_Destroy");
2655 int *coordinate_grid_parts =
new int[this->numLocalCoords];
2656 switch (this->predistribution){
2658 this->predistributeRCB(coordinate_grid_parts);
2662 this->predistributeMJ(coordinate_grid_parts);
2666 blockPartition(coordinate_grid_parts);
2669 this->distribute_points(coordinate_grid_parts);
2671 delete []coordinate_grid_parts;
2682 return this->numWeightsPerCoord;
2685 return this->coordinate_dimension;
2688 return this->numLocalCoords;
2691 return this->numGlobalCoords;
2695 return this->coords;
2703 for(
int ii = 0; ii < this->coordinate_dimension; ++ii){
2704 #ifdef HAVE_ZOLTAN2_OMP 2705 #pragma omp parallel for 2707 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2708 c[ii][i] = this->coords[ii][i];
2714 for(
int ii = 0; ii < this->numWeightsPerCoord; ++ii){
2715 #ifdef HAVE_ZOLTAN2_OMP 2716 #pragma omp parallel for 2718 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2719 w[ii][i] = this->wghts[ii][i];
RectangleHole(CoordinatePoint< T > center_, T edge_x_, T edge_y_)
virtual weighttype get1DWeight(T x)
virtual bool isInArea(CoordinatePoint< T > dot)
void getCoords(void *data, int numGid, int numLid, int numObj, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int dim, double *coords_, int *ierr)
CoordinateGridDistribution(gno_t alongX, gno_t alongY, gno_t alongZ, int dim, T l_x, T r_x, T l_y, T r_y, T l_z, T r_z, int myRank_, int wSize)
virtual weighttype get1DWeight(T x)=0
std::vector< std::vector< float > > weights
void print(T x, T y, T z)
virtual ~CoordinateDistribution()
void getLocalCoordinatesCopy(scalar_t **c)
void GetPoints(lno_t requestedPointcount, CoordinatePoint< T > *points, Hole< T > **holes, lno_t holeCount, float *sharedRatios_, int myRank)
virtual bool isInArea(CoordinatePoint< T > dot)
virtual bool isInArea(CoordinatePoint< T > dot)
int predistributeRCB(int *coordinate_grid_parts)
virtual CoordinatePoint< T > getPoint(gno_t pindex, unsigned int &state)
virtual weighttype get2DWeight(T x, T y)
void GetPoints(lno_t requestedPointcount, T **coords, lno_t tindex, Hole< T > **holes, lno_t holeCount, float *sharedRatios_, int myRank)
Hole(CoordinatePoint< T > center_, T edge1_, T edge2_, T edge3_)
scalar_t ** getLocalWeightsView()
SphereHole(CoordinatePoint< T > center_, T edge_)
CircleHole(CoordinatePoint< T > center_, T edge_)
int getCoordinateDimension()
virtual ~RectangularPrismHole()
Defines the PartitioningSolution class.
void perturb_data(float used_perturbation_ratio, int scale)
Start perturbation function##########################//
virtual weighttype get3DWeight(T x, T y, T z)
virtual weighttype get3DWeight(T x, T y, T z)=0
Defines the XpetraMultiVectorAdapter.
void blockPartition(int *coordinate_grid_parts)
RectangularPrismHole(CoordinatePoint< T > center_, T edge_x_, T edge_y_, T edge_z_)
CoordinateNormalDistribution(gno_t np_, int dim, CoordinatePoint< T > center_, T sd_x, T sd_y, T sd_z, int wSize)
virtual bool isInArea(CoordinatePoint< T > dot)
virtual bool isInArea(CoordinatePoint< T > dot)
SteppedEquation(T a1_, T a2_, T a3_, T b1_, T b2_, T b3_, T c_, T x1_, T y1_, T z1_, T *steps_, T *values_, int stepCount_)
int getNumObj(void *data, int *ierr)
const std::string shapes[]
scalar_t ** getLocalCoordinatesView()
virtual weighttype getWeight(CoordinatePoint< T > p)
int getDim(void *data, int *ierr)
virtual ~SteppedEquation()
void getObjList(void *data, int numGid, int numLid, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int num_wgts, float *obj_wgts, int *ierr)
const std::string weight_distribution_string
virtual CoordinatePoint< T > getPoint(gno_t pindex, unsigned int &state)
virtual weighttype get2DWeight(T x, T y)=0
An adapter for Xpetra::MultiVector.
void getBestSurface(int remaining, int *dimProcs, int dim, int currentDim, int &bestSurface, int *bestDimProcs)
Start Predistribution functions######################//
Expression is a generic following method.
std::string toString(tt obj)
Converts the given object to string.
CoordinatePoint< T > center
void distribute_points(int *coordinate_grid_parts)
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
virtual ~WeightDistribution()
PartitioningProblem sets up partitioning problems for the user.
GeometricGenerator(Teuchos::ParameterList ¶ms, const RCP< const Teuchos::Comm< int > > &comm_)
int getNumWeights()
##END Predistribution functions######################//
CoordinatePoint< T > center
#define INVALID_SHAPE_ARG(SHAPE, REQUIRED)
lno_t getNumLocalCoords()
virtual ~CoordinateNormalDistribution()
#define INVALIDSHAPE(STR, DIM)
#define CATCH_EXCEPTIONS(pp)
Defines the PartitioningProblem class.
CoordinateDistribution(gno_t np_, int dim, int wSize)
int predistributeMJ(int *coordinate_grid_parts)
virtual ~CoordinateGridDistribution()
gno_t getNumGlobalCoords()
SquareHole(CoordinatePoint< T > center_, T edge_)
void solve(bool updateInputData=true)
Direct the problem to create a solution.
CubeHole(CoordinatePoint< T > center_, T edge_)
void getLocalWeightsCopy(scalar_t **w)
virtual bool isInArea(CoordinatePoint< T > dot)
void getMinMaxCoords(scalar_t *globalMaxCoords, scalar_t *globalMinCoords)
#define MEMORY_CHECK(iPrint, msg)