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;
81 namespace GeometricGen{
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()->getLocalElementList().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>
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))
328 template <
typename T,
typename weighttype>
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>
490 float *sharedRatios_,
int myRank){
492 for (
int i = 0; i < myRank; ++i){
495 if (i < this->numPoints % this->
worldSize){
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;
579 float *sharedRatios_,
int myRank){
581 for (
int i = 0; i < myRank; ++i){
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;
632 coords[1][cnt + tindex] = p.
y;
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) :
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 T polarsqrt, normalsquared, normal1, normal2;
698 normal1=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
699 normal2=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
700 normalsquared=normal1*normal1+normal2*normal2;
701 }
while ( normalsquared>=1.0 || normalsquared == 0.0);
703 polarsqrt=sqrt(-2.0*log(normalsquared)/normalsquared);
704 return normal2*polarsqrt*sd + center_;
708 template <
typename T,
typename lno_t,
typename gno_t>
728 T l_z, T r_z,
int wSize ) :
739 for(
int i = 0; i < this->
dimension; ++i){
751 throw "unsupported dimension";
759 T uniformDist(T a, T b,
unsigned int &state)
761 return (b-a)*(rand_r(&state) / double(RAND_MAX)) + a;
765 template <
typename T,
typename lno_t,
typename gno_t>
790 T l_x, T r_x, T l_y, T r_y, T l_z, T r_z ,
791 int myRank_,
int wSize) :
891 template <
typename scalar_t,
typename lno_t,
typename gno_t,
typename node_t>
896 int coordinate_dimension;
897 gno_t numGlobalCoords;
898 lno_t numLocalCoords;
899 float *loadDistributions;
901 bool distinctCoordSet;
903 int distributionCount;
908 int numWeightsPerCoord;
910 RCP<const Teuchos::Comm<int> > comm;
922 float perturbation_ratio;
924 typedef Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> tMVector_t;
925 typedef Tpetra::Map<lno_t, gno_t, node_t> tMap_t;
928 template <
typename tt>
929 tt getParamVal(
const Teuchos::ParameterEntry& pe,
const std::string ¶mname){
932 returnVal = Teuchos::getValue<tt>(pe);
940 int countChar (std::string inStr,
char cntChar){
942 for (
unsigned int i = 0; i < inStr.size(); ++i){
943 if (inStr[i] == cntChar) {
950 template <
typename tt>
951 tt fromString(std::string obj){
952 std::stringstream ss (std::stringstream::in | std::stringstream::out);
957 throw "Cannot convert string " + obj;
963 void splitString(std::string inStr,
char splitChar, std::string *outSplittedStr){
964 std::stringstream ss (std::stringstream::in | std::stringstream::out);
968 std::string tmp =
"";
969 std::getline(ss, tmp ,splitChar);
970 outSplittedStr[cnt++] = tmp;
1007 void getCoordinateDistributions(std::string coordinate_distributions){
1008 if(coordinate_distributions ==
""){
1009 throw "At least one distribution function must be provided to coordinate generator.";
1012 int argCnt = this->countChar(coordinate_distributions,
',') + 1;
1013 std::string *splittedStr =
new std::string[argCnt];
1014 splitString(coordinate_distributions,
',', splittedStr);
1016 for(
int i = 0; i < argCnt; ){
1019 std::string distName = splittedStr[i++];
1021 if(distName ==
"NORMAL"){
1023 if (this->coordinate_dimension == 3){
1026 if(i + reqArg > argCnt) {
1027 std::string tmp = Teuchos::toString<int>(reqArg);
1030 np_ = fromString<gno_t>(splittedStr[i++]);
1033 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1034 pp.y = fromString<scalar_t>(splittedStr[i++]);
1036 if(this->coordinate_dimension == 3){
1037 pp.z = fromString<scalar_t>(splittedStr[i++]);
1040 scalar_t sd_x = fromString<scalar_t>(splittedStr[i++]);
1041 scalar_t sd_y = fromString<scalar_t>(splittedStr[i++]);
1043 if(this->coordinate_dimension == 3){
1044 sd_z = fromString<scalar_t>(splittedStr[i++]);
1048 }
else if(distName ==
"UNIFORM" ){
1050 if (this->coordinate_dimension == 3){
1053 if(i + reqArg > argCnt) {
1054 std::string tmp = Teuchos::toString<int>(reqArg);
1057 np_ = fromString<gno_t>(splittedStr[i++]);
1058 scalar_t l_x = fromString<scalar_t>(splittedStr[i++]);
1059 scalar_t r_x = fromString<scalar_t>(splittedStr[i++]);
1060 scalar_t l_y = fromString<scalar_t>(splittedStr[i++]);
1061 scalar_t r_y = fromString<scalar_t>(splittedStr[i++]);
1063 scalar_t l_z = 0, r_z = 0;
1065 if(this->coordinate_dimension == 3){
1066 l_z = fromString<scalar_t>(splittedStr[i++]);
1067 r_z = fromString<scalar_t>(splittedStr[i++]);
1070 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 );
1071 }
else if (distName ==
"GRID"){
1073 if(this->coordinate_dimension == 3){
1076 if(i + reqArg > argCnt) {
1077 std::string tmp = Teuchos::toString<int>(reqArg);
1081 gno_t np_x = fromString<gno_t>(splittedStr[i++]);
1082 gno_t np_y = fromString<gno_t>(splittedStr[i++]);
1086 if(this->coordinate_dimension == 3){
1087 np_z = fromString<gno_t>(splittedStr[i++]);
1090 np_ = np_x * np_y * np_z;
1091 scalar_t l_x = fromString<scalar_t>(splittedStr[i++]);
1092 scalar_t r_x = fromString<scalar_t>(splittedStr[i++]);
1093 scalar_t l_y = fromString<scalar_t>(splittedStr[i++]);
1094 scalar_t r_y = fromString<scalar_t>(splittedStr[i++]);
1096 scalar_t l_z = 0, r_z = 0;
1098 if(this->coordinate_dimension == 3){
1099 l_z = fromString<scalar_t>(splittedStr[i++]);
1100 r_z = fromString<scalar_t>(splittedStr[i++]);
1103 if(np_x < 1 || np_z < 1 || np_y < 1 ){
1104 throw "Provide at least 1 point along each dimension for grid test.";
1108 (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);
1112 std::string tmp = Teuchos::toString<int>(this->coordinate_dimension);
1115 this->numGlobalCoords += (
gno_t) np_;
1117 delete []splittedStr;
1120 void getProcLoadDistributions(std::string proc_load_distributions){
1122 this->loadDistributions =
new float[this->worldSize];
1123 if(proc_load_distributions ==
""){
1124 float r = 1.0 / this->worldSize;
1125 for(
int i = 0; i < this->worldSize; ++i){
1126 this->loadDistributions[i] = r;
1131 int argCnt = this->countChar(proc_load_distributions,
',') + 1;
1132 if(argCnt != this->worldSize) {
1133 throw "Invalid parameter count load distributions. Given " + Teuchos::toString<int>(argCnt) +
" processor size is " + Teuchos::toString<int>(this->worldSize);
1135 std::string *splittedStr =
new std::string[argCnt];
1136 splitString(proc_load_distributions,
',', splittedStr);
1137 for(
int i = 0; i < argCnt; ++i){
1138 this->loadDistributions[i] = fromString<float>(splittedStr[i]);
1140 delete []splittedStr;
1144 for(
int i = 0; i < this->worldSize; ++i){
1145 sum += this->loadDistributions[i];
1148 throw "Processor load ratios do not sum to 1.0.";
1154 void getHoles(std::string holeDescription){
1155 if(holeDescription ==
""){
1159 int argCnt = this->countChar(holeDescription,
',') + 1;
1160 std::string *splittedStr =
new std::string[argCnt];
1161 splitString(holeDescription,
',', splittedStr);
1163 for(
int i = 0; i < argCnt; ){
1166 std::string shapeName = splittedStr[i++];
1167 if(shapeName ==
"SQUARE" && this->coordinate_dimension == 2){
1168 if(i + 3 > argCnt) {
1172 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1173 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1174 scalar_t edge = fromString<scalar_t>(splittedStr[i++]);
1176 }
else if(shapeName ==
"RECTANGLE" && this->coordinate_dimension == 2){
1177 if(i + 4 > argCnt) {
1181 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1182 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1183 scalar_t edgex = fromString<scalar_t>(splittedStr[i++]);
1184 scalar_t edgey = fromString<scalar_t>(splittedStr[i++]);
1187 }
else if(shapeName ==
"CIRCLE" && this->coordinate_dimension == 2){
1188 if(i + 3 > argCnt) {
1192 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1193 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1194 scalar_t r = fromString<scalar_t>(splittedStr[i++]);
1196 }
else if(shapeName ==
"CUBE" && this->coordinate_dimension == 3){
1197 if(i + 4 > argCnt) {
1201 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1202 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1203 pp.
z = fromString<scalar_t>(splittedStr[i++]);
1204 scalar_t edge = fromString<scalar_t>(splittedStr[i++]);
1206 }
else if(shapeName ==
"RECTANGULAR_PRISM" && this->coordinate_dimension == 3){
1207 if(i + 6 > argCnt) {
1211 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1212 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1213 pp.
z = fromString<scalar_t>(splittedStr[i++]);
1214 scalar_t edgex = fromString<scalar_t>(splittedStr[i++]);
1215 scalar_t edgey = fromString<scalar_t>(splittedStr[i++]);
1216 scalar_t edgez = fromString<scalar_t>(splittedStr[i++]);
1219 }
else if(shapeName ==
"SPHERE" && this->coordinate_dimension == 3){
1220 if(i + 4 > argCnt) {
1224 pp.
x = fromString<scalar_t>(splittedStr[i++]);
1225 pp.
y = fromString<scalar_t>(splittedStr[i++]);
1226 pp.
z = fromString<scalar_t>(splittedStr[i++]);
1227 scalar_t r = fromString<scalar_t>(splittedStr[i++]);
1230 std::string tmp = Teuchos::toString<int>(this->coordinate_dimension);
1234 delete [] splittedStr;
1237 void getWeightDistribution(std::string *weight_distribution_arr,
int wdimension){
1242 std::string weight_distribution = weight_distribution_arr[ii];
1243 if(weight_distribution ==
"")
continue;
1244 if(wcount == wdimension) {
1245 throw "Weight Dimension is provided as " + Teuchos::toString<int>(wdimension) +
". More weight distribution is provided.";
1248 int count = this->countChar(weight_distribution,
' ');
1249 std::string *splittedStr =
new string[count + 1];
1250 this->splitString(weight_distribution,
' ', splittedStr);
1253 scalar_t a1=0,a2=0,a3=0;
1254 scalar_t x1=0,y1=0,z1=0;
1255 scalar_t b1=1,b2=1,b3=1;
1256 scalar_t *steps = NULL;
1257 scalar_t *values= NULL;
1261 for (
int i = 1; i < count + 1; ++i){
1262 int assignmentCount = this->countChar(splittedStr[i],
'=');
1263 if (assignmentCount != 1){
1264 throw "Error at the argument " + splittedStr[i];
1266 std::string parameter_vs_val[2];
1267 this->splitString(splittedStr[i],
'=', parameter_vs_val);
1268 std::string parameter = parameter_vs_val[0];
1269 std::string value = parameter_vs_val[1];
1272 if (parameter ==
"a1"){
1273 a1 = this->fromString<scalar_t>(value);
1275 else if (parameter ==
"a2"){
1276 if(this->coordinate_dimension > 1){
1277 a2 = this->fromString<scalar_t>(value);
1280 throw parameter+
" argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1284 else if (parameter ==
"a3"){
1285 if(this->coordinate_dimension > 2){
1286 a3 = this->fromString<scalar_t>(value);
1289 throw parameter+
" argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1292 else if (parameter ==
"b1"){
1293 b1 = this->fromString<scalar_t>(value);
1295 else if (parameter ==
"b2"){
1296 if(this->coordinate_dimension > 1){
1297 b2 = this->fromString<scalar_t>(value);
1300 throw parameter+
" argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1303 else if (parameter ==
"b3"){
1305 if(this->coordinate_dimension > 2){
1306 b3 = this->fromString<scalar_t>(value);
1309 throw parameter+
" argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1312 else if (parameter ==
"c"){
1313 c = this->fromString<scalar_t>(value);
1315 else if (parameter ==
"x1"){
1316 x1 = this->fromString<scalar_t>(value);
1318 else if (parameter ==
"y1"){
1319 if(this->coordinate_dimension > 1){
1320 y1 = this->fromString<scalar_t>(value);
1323 throw parameter+
" argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1326 else if (parameter ==
"z1"){
1327 if(this->coordinate_dimension > 2){
1328 z1 = this->fromString<scalar_t>(value);
1331 throw parameter+
" argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1334 else if (parameter ==
"steps"){
1335 stepCount = this->countChar(value,
',') + 1;
1336 std::string *stepstr =
new std::string[stepCount];
1337 this->splitString(value,
',', stepstr);
1338 steps =
new scalar_t[stepCount];
1339 for (
int j = 0; j < stepCount; ++j){
1340 steps[j] = this->fromString<scalar_t>(stepstr[j]);
1344 else if (parameter ==
"values"){
1345 valueCount = this->countChar(value,
',') + 1;
1346 std::string *stepstr =
new std::string[valueCount];
1347 this->splitString(value,
',', stepstr);
1348 values =
new scalar_t[valueCount];
1349 for (
int j = 0; j < valueCount; ++j){
1350 values[j] = this->fromString<scalar_t>(stepstr[j]);
1355 throw "Invalid parameter name at " + splittedStr[i];
1359 delete []splittedStr;
1360 if(stepCount + 1!= valueCount){
1361 throw "Step count: " + Teuchos::toString<int>(stepCount) +
" must be 1 less than value count: " + Teuchos::toString<int>(valueCount);
1365 this->wd[wcount++] =
new SteppedEquation<scalar_t,scalar_t>(a1, a2, a3, b1, b2, b3, c, x1, y1, z1, steps, values, stepCount);
1373 if(wcount != this->numWeightsPerCoord){
1374 throw "Weight Dimension is provided as " + Teuchos::toString<int>(wdimension) +
". But " + Teuchos::toString<int>(wcount)+
" weight distributions are provided.";
1378 void parseParams(
const Teuchos::ParameterList & params){
1380 std::string holeDescription =
"";
1381 std::string proc_load_distributions =
"";
1382 std::string distinctDescription =
"";
1383 std::string coordinate_distributions =
"";
1386 numWeightsPerCoord_parameters[i] =
"";
1390 for (Teuchos::ParameterList::ConstIterator pit = params.begin(); pit != params.end(); ++pit ){
1391 const std::string ¶mName = params.name(pit);
1393 const Teuchos::ParameterEntry &pe = params.entry(pit);
1395 if(paramName.find(
"hole-") == 0){
1396 if(holeDescription ==
""){
1397 holeDescription = getParamVal<std::string>(pe, paramName);
1400 holeDescription +=
","+getParamVal<std::string>(pe, paramName);
1403 else if(paramName.find(
"distribution-") == 0){
1404 if(coordinate_distributions ==
"")
1405 coordinate_distributions = getParamVal<std::string>(pe, paramName);
1407 coordinate_distributions +=
","+getParamVal<std::string>(pe, paramName);
1414 int dash_pos = weight_dist_param.find(
"-");
1415 std::string distribution_index_string = weight_dist_param.substr(0, dash_pos);
1416 int distribution_index = fromString<int>(distribution_index_string);
1418 if(distribution_index >= MAX_WEIGHT_DIM){
1419 throw "Given distribution index:" + distribution_index_string +
" larger than maximum allowed number of weights:" + Teuchos::toString<int>(
MAX_WEIGHT_DIM);
1421 numWeightsPerCoord_parameters[distribution_index] +=
" " + weight_dist_param.substr(dash_pos + 1)+
"="+ getParamVal<std::string>(pe, paramName);
1423 else if(paramName ==
"dim"){
1424 int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
1425 if(dim < 2 || dim > 3){
1428 this->coordinate_dimension = dim;
1431 else if(paramName ==
"wdim"){
1432 int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
1433 if(dim < 1 && dim > MAX_WEIGHT_DIM){
1436 this->numWeightsPerCoord = dim;
1439 else if(paramName ==
"predistribution"){
1440 int pre_distribution = fromString<int>(getParamVal<std::string>(pe, paramName));
1441 if(pre_distribution < 0 || pre_distribution > 3){
1444 this->predistribution = pre_distribution;
1447 else if(paramName ==
"perturbation_ratio"){
1448 float perturbation = fromString<float>(getParamVal<std::string>(pe, paramName));
1449 if(perturbation < 0 && perturbation > 1){
1452 this->perturbation_ratio = perturbation;
1457 else if(paramName ==
"proc_load_distributions"){
1458 proc_load_distributions = getParamVal<std::string>(pe, paramName);
1459 this->loadDistSet =
true;
1462 else if(paramName ==
"distinct_coordinates"){
1463 distinctDescription = getParamVal<std::string>(pe, paramName);
1464 if(distinctDescription ==
"T"){
1465 this->distinctCoordSet =
true;
1466 }
else if(distinctDescription ==
"F"){
1467 this->distinctCoordSet =
true;
1469 throw "Invalid parameter for distinct_coordinates: " + distinctDescription +
". Candidates: T or F." ;
1473 else if(paramName ==
"out_file"){
1474 this->outfile = getParamVal<std::string>(pe, paramName);
1483 if(this->coordinate_dimension == 0){
1484 throw "Dimension must be provided to coordinate generator.";
1503 if (this->loadDistSet && this->distinctCoordSet){
1504 throw "distinct_coordinates and proc_load_distributions parameters cannot be satisfied together.";
1506 this->getHoles(holeDescription);
1508 this->getProcLoadDistributions(proc_load_distributions);
1509 this->getCoordinateDistributions(coordinate_distributions);
1510 this->getWeightDistribution(numWeightsPerCoord_parameters, this->numWeightsPerCoord);
1517 catch(std::string &s){
1518 throw std::runtime_error(s);
1521 catch(
const char *s){
1522 throw std::runtime_error(s);
1529 for (
int i = 0; i < this->holeCount; ++i){
1530 delete this->holes[i];
1536 delete []loadDistributions;
1538 if(coordinateDistributions){
1540 for (
int i = 0; i < this->distributionCount; ++i){
1541 delete this->coordinateDistributions[i];
1543 free (this->coordinateDistributions);
1546 for (
int i = 0; i < this->numWeightsPerCoord; ++i){
1552 if(this->numWeightsPerCoord){
1553 for(
int i = 0; i < this->numWeightsPerCoord; ++i)
1554 delete [] this->wghts[i];
1555 delete []this->wghts;
1557 if(this->coordinate_dimension){
1558 for(
int i = 0; i < this->coordinate_dimension; ++i)
1559 delete [] this->coords[i];
1560 delete [] this->coords;
1566 std::cout <<
"\nGeometric Generator Parameter File Format:" << std::endl;
1567 std::cout <<
"- dim=Coordinate Dimension: 2 or 3" << std::endl;
1568 std::cout <<
"- Available distributions:" << std::endl;
1569 std::cout <<
"\tUNIFORM: -> distribution-1=UNIFORM,NUMCOORDINATES,XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << std::endl;
1570 std::cout <<
"\tGRID: -> distribution-2=GRID,XLENGTH,YLENGTH{,ZLENGTH},XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << std::endl;
1571 std::cout <<
"\tNORMAL: -> distribution-3=NORMAL,XCENTER,YCENTER{,ZCENTER},XSD,YSD,{,ZSD}" << std::endl;
1572 std::cout <<
"- wdim=numWeightsPerCoord: There should be as many weight function as number of weights per coord." << std::endl;
1573 std::cout <<
"- Weight Equation: w = (a1 * (x - x1)^b1) + (a2 * (y - y1)^b2) + (a3 * (z - z1)^b3) + c" << std::endl;
1574 std::cout <<
"Parameter settings:" << std::endl;
1575 std::cout <<
"\tWeightDistribution-1-a1=a1 " << std::endl;
1576 std::cout <<
"\tWeightDistribution-1-a2=a2 " << std::endl;
1577 std::cout <<
"\tWeightDistribution-1-a3=a3 " << std::endl;
1578 std::cout <<
"\tWeightDistribution-1-b1=b1 " << std::endl;
1579 std::cout <<
"\tWeightDistribution-1-b2=b2 " << std::endl;
1580 std::cout <<
"\tWeightDistribution-1-b3=b3 " << std::endl;
1581 std::cout <<
"\tWeightDistribution-1-x1=x1 " << std::endl;
1582 std::cout <<
"\tWeightDistribution-1-y1=y1 " << std::endl;
1583 std::cout <<
"\tWeightDistribution-1-z1=z1 " << std::endl;
1584 std::cout <<
"\tWeightDistribution-1-c=c " << std::endl;
1585 std::cout <<
"\tIt is possible to set step function to the result of weight equation." << std::endl;
1586 std::cout <<
"\tWeightDistribution-1-steps=step1,step2,step3:increasing order" << std::endl;
1587 std::cout <<
"\tWeightDistribution-1-values=val1,val2,val3,val4." << std::endl;
1588 std::cout <<
"\t\tIf w < step1 -> w = val1" << std::endl;
1589 std::cout <<
"\t\tElse if w < step2 -> w = val2" << std::endl;
1590 std::cout <<
"\t\tElse if w < step3 -> w = val3" << std::endl;
1591 std::cout <<
"\t\tElse -> w = val4" << std::endl;
1592 std::cout <<
"- Holes:" << std::endl;
1593 std::cout <<
"\thole-1:SPHERE,XCENTER,YCENTER,ZCENTER,RADIUS (only for dim=3)" << std::endl;
1594 std::cout <<
"\thole-2:CUBE,XCENTER,YCENTER,ZCENTER,EDGE (only for dim=3)" << std::endl;
1595 std::cout <<
"\thole-3:RECTANGULAR_PRISM,XCENTER,YCENTER,ZCENTER,XEDGE,YEDGE,ZEDGE (only for dim=3)" << std::endl;
1596 std::cout <<
"\thole-4:SQUARE,XCENTER,YCENTER,EDGE (only for dim=2)" << std::endl;
1597 std::cout <<
"\thole-5:RECTANGLE,XCENTER,YCENTER,XEDGE,YEDGE (only for dim=2)" << std::endl;
1598 std::cout <<
"\thole-6:CIRCLE,XCENTER,YCENTER,RADIUS (only for dim=2)" << std::endl;
1599 std::cout <<
"- out_file:out_file_path : if provided output will be written to files." << std::endl;
1600 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;
1601 std::cout <<
"- predistribution:distribution_option. Predistribution of the coordinates to the processors. 0 for NONE, 1 RCB, 2 MJ, 3 BLOCK." << std::endl;
1602 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;
1609 this->coordinate_dimension = 0;
1610 this->numWeightsPerCoord = 0;
1611 this->worldSize = comm_->getSize();
1612 this->numGlobalCoords = 0;
1613 this->loadDistributions = NULL;
1615 this->coordinateDistributions = NULL;
1616 this->holeCount = 0;
1617 this->distributionCount = 0;
1619 this->predistribution = 0;
1620 this->perturbation_ratio = 0;
1631 this->loadDistSet =
false;
1632 this->distinctCoordSet =
false;
1633 this->myRank = comm_->getRank();
1636 this->parseParams(params);
1638 catch(std::string &s){
1642 throw std::runtime_error(s);
1645 catch(
const char *s){
1649 throw std::runtime_error(s);
1653 lno_t myPointCount = 0;
1654 this->numGlobalCoords = 0;
1656 gno_t prefixSum = 0;
1657 for(
int i = 0; i < this->distributionCount; ++i){
1658 for(
int ii = 0; ii < this->worldSize; ++ii){
1659 lno_t increment =
lno_t (this->coordinateDistributions[i]->numPoints * this->loadDistributions[ii]);
1660 if (
gno_t(ii) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1663 this->numGlobalCoords += increment;
1665 prefixSum += increment;
1668 myPointCount +=
lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
1669 if (
gno_t(myRank) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1674 this->coords =
new scalar_t *[this->coordinate_dimension];
1675 for(
int i = 0; i < this->coordinate_dimension; ++i){
1676 this->coords[i] =
new scalar_t[myPointCount];
1679 for (
int ii = 0; ii < this->coordinate_dimension; ++ii){
1680 #ifdef HAVE_ZOLTAN2_OMP
1681 #pragma omp parallel for
1683 for(
lno_t i = 0; i < myPointCount; ++i){
1684 this->coords[ii][i] = 0;
1688 this->numLocalCoords = 0;
1689 srand ((myRank + 1) * this->numLocalCoords);
1690 for (
int i = 0; i < distributionCount; ++i){
1692 lno_t requestedPointCount =
lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
1693 if (
gno_t(myRank) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1694 requestedPointCount += 1;
1698 this->coordinateDistributions[i]->
GetPoints(requestedPointCount,this->coords, this->numLocalCoords, this->holes, this->holeCount, this->loadDistributions, myRank);
1699 this->numLocalCoords += requestedPointCount;
1721 if (this->predistribution){
1728 if (this->perturbation_ratio > 0){
1750 if (this->distinctCoordSet){
1754 if(this->outfile !=
""){
1756 std::ofstream myfile;
1757 myfile.open ((this->outfile + Teuchos::toString<int>(myRank)).c_str());
1758 for(
lno_t i = 0; i < this->numLocalCoords; ++i){
1760 myfile << this->coords[0][i];
1761 if(this->coordinate_dimension > 1){
1762 myfile <<
" " << this->coords[1][i];
1764 if(this->coordinate_dimension > 2){
1765 myfile <<
" " << this->coords[2][i];
1767 myfile << std::endl;
1771 if (this->myRank == 0){
1772 std::ofstream gnuplotfile(
"gnu.gnuplot");
1773 for(
int i = 0; i < this->worldSize; ++i){
1775 if (this->coordinate_dimension == 2){
1781 gnuplotfile << s <<
" \"" << (this->outfile + Teuchos::toString<int>(i)) <<
"\"" << std::endl;
1783 gnuplotfile <<
"pause -1" << std::endl;
1784 gnuplotfile.close();
1797 if (this->numWeightsPerCoord > 0){
1798 this->wghts =
new scalar_t *[this->numWeightsPerCoord];
1799 for(
int i = 0; i < this->numWeightsPerCoord; ++i){
1800 this->wghts[i] =
new scalar_t[this->numLocalCoords];
1804 for(
int ii = 0; ii < this->numWeightsPerCoord; ++ii){
1805 switch(this->coordinate_dimension){
1808 #ifdef HAVE_ZOLTAN2_OMP
1809 #pragma omp parallel for
1811 for (
lno_t i = 0; i < this->numLocalCoords; ++i){
1812 this->wghts[ii][i] = this->wd[ii]->
get1DWeight(this->coords[0][i]);
1818 #ifdef HAVE_ZOLTAN2_OMP
1819 #pragma omp parallel for
1821 for (
lno_t i = 0; i < this->numLocalCoords; ++i){
1822 this->wghts[ii][i] = this->wd[ii]->
get2DWeight(this->coords[0][i], this->coords[1][i]);
1828 #ifdef HAVE_ZOLTAN2_OMP
1829 #pragma omp parallel for
1831 for (
lno_t i = 0; i < this->numLocalCoords; ++i){
1832 this->wghts[ii][i] = this->wd[ii]->
get3DWeight(this->coords[0][i], this->coords[1][i], this->coords[2][i]);
1844 scalar_t *maxCoords=
new scalar_t [this->coordinate_dimension];
1845 scalar_t *minCoords=
new scalar_t [this->coordinate_dimension];
1846 for (
int dim = 0; dim < this->coordinate_dimension; ++dim){
1847 minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
1848 for (
lno_t i = 1; i < this->numLocalCoords; ++i){
1849 if (minCoords[dim] > this->coords[dim][i]){
1850 minCoords[dim] = this->coords[dim][i];
1853 if (maxCoords[dim] < this->coords[dim][i]){
1854 maxCoords[dim] = this->coords[dim][i];
1861 scalar_t center = (maxCoords[dim] + minCoords[dim]) / 2;
1863 minCoords[dim] = center - (center - minCoords[dim]) * scale;
1864 maxCoords[dim] = (maxCoords[dim] - center) * scale + center;
1868 gno_t numLocalToPerturb =
gno_t(this->numLocalCoords*used_perturbation_ratio);
1870 for (
int dim = 0; dim < this->coordinate_dimension; ++dim){
1871 scalar_t range = maxCoords[dim] - minCoords[dim];
1872 for (
gno_t i = 0; i < numLocalToPerturb; ++i){
1873 this->coords[dim][i] = (rand() / double (RAND_MAX)) * (range) + minCoords[dim];
1886 void getBestSurface (
int remaining,
int *dimProcs,
int dim,
int currentDim,
int &bestSurface,
int *bestDimProcs){
1888 if (currentDim < dim - 1){
1890 while(ipx <= remaining) {
1891 if(remaining % ipx == 0) {
1892 int nremain = remaining / ipx;
1893 dimProcs[currentDim] = ipx;
1894 getBestSurface (nremain, dimProcs, dim, currentDim + 1, bestSurface, bestDimProcs);
1900 dimProcs[currentDim] = remaining;
1902 for (
int i = 0; i < dim; ++i) surface += dimProcs[i];
1903 if (surface < bestSurface){
1904 bestSurface = surface;
1905 for (
int i = 0; i < dim; ++i) bestDimProcs[i] = dimProcs[i];
1913 scalar_t *maxCoords=
new scalar_t [this->coordinate_dimension];
1914 scalar_t *minCoords=
new scalar_t [this->coordinate_dimension];
1915 for (
int dim = 0; dim < this->coordinate_dimension; ++dim){
1916 minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
1917 for (
lno_t i = 1; i < this->numLocalCoords; ++i){
1918 if (minCoords[dim] > this->coords[dim][i]){
1919 minCoords[dim] = this->coords[dim][i];
1922 if (maxCoords[dim] < this->coords[dim][i]){
1923 maxCoords[dim] = this->coords[dim][i];
1928 reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MAX,
1929 this->coordinate_dimension,
1934 reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MIN,
1935 this->coordinate_dimension,
1949 typedef SSIZE_T ssize_t;
1956 scalar_t *maxCoords=
new scalar_t [this->coordinate_dimension];
1957 scalar_t *minCoords=
new scalar_t [this->coordinate_dimension];
1965 int remaining = this->worldSize;
1966 int coord_dim = this->coordinate_dimension;
1967 int *dimProcs =
new int[coord_dim];
1968 int *bestDimProcs =
new int[coord_dim];
1971 int bestSurface = 0;
1972 dimProcs[0] = remaining;
1973 for (
int i = 1; i < coord_dim; ++i) dimProcs[i] = 1;
1974 for (
int i = 0; i < coord_dim; ++i) bestSurface += dimProcs[i];
1975 for (
int i = 0; i < coord_dim; ++i) bestDimProcs[i] = dimProcs[i];
1977 getBestSurface ( remaining, dimProcs, coord_dim, currentDim, bestSurface, bestDimProcs);
1985 int *shiftProcCount =
new int[coord_dim];
1989 int remainingProc = this->worldSize;
1990 for (
int dim = 0; dim < coord_dim; ++dim){
1991 remainingProc = remainingProc / bestDimProcs[dim];
1992 shiftProcCount[dim] = remainingProc;
1995 scalar_t *dim_slices =
new scalar_t[coord_dim];
1996 for (
int dim = 0; dim < coord_dim; ++dim){
1997 scalar_t dim_range = maxCoords[dim] - minCoords[dim];
1998 scalar_t slice = dim_range / bestDimProcs[dim];
1999 dim_slices[dim] = slice;
2006 gno_t *numPointsInParts =
new gno_t[this->worldSize];
2007 gno_t *numGlobalPointsInParts =
new gno_t[this->worldSize];
2008 gno_t *numPointsInPartsInclusiveUptoMyIndex =
new gno_t[this->worldSize];
2010 gno_t *partBegins =
new gno_t [this->worldSize];
2011 gno_t *partNext =
new gno_t[this->numLocalCoords];
2014 for (
lno_t i = 0; i < this->numLocalCoords; ++i){
2017 for (
int i = 0; i < this->worldSize; ++i) {
2021 for (
int i = 0; i < this->worldSize; ++i)
2022 numPointsInParts[i] = 0;
2024 for (
lno_t i = 0; i < this->numLocalCoords; ++i){
2026 for (
int dim = 0; dim < coord_dim; ++dim){
2027 int shift = int ((this->coords[dim][i] - minCoords[dim]) / dim_slices[dim]);
2028 if (shift >= bestDimProcs[dim]){
2029 shift = bestDimProcs[dim] - 1;
2031 shift = shift * shiftProcCount[dim];
2034 numPointsInParts[partIndex] += 1;
2035 coordinate_grid_parts[i] = partIndex;
2037 partNext[i] = partBegins[partIndex];
2038 partBegins[partIndex] = i;
2045 reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2048 numGlobalPointsInParts);
2051 Teuchos::scan<int,gno_t>(
2052 *(this->comm), Teuchos::REDUCE_SUM,
2055 numPointsInPartsInclusiveUptoMyIndex
2074 gno_t optimal_num =
gno_t(this->numGlobalCoords/
double(this->worldSize)+0.5);
2076 if (this->myRank == 0){
2077 gno_t totalSize = 0;
2078 for (
int i = 0; i < this->worldSize; ++i){
2079 std::cout <<
"me:" << this->myRank <<
" NumPoints in part:" << i <<
" is: " << numGlobalPointsInParts[i] << std::endl;
2080 totalSize += numGlobalPointsInParts[i];
2082 std::cout <<
"Total:" << totalSize <<
" ng:" << this->numGlobalCoords << std::endl;
2083 std::cout <<
"optimal_num:" << optimal_num << std::endl;
2086 ssize_t *extraInPart =
new ssize_t [this->worldSize];
2088 ssize_t extraExchange = 0;
2089 for (
int i = 0; i < this->worldSize; ++i){
2090 extraInPart[i] = numGlobalPointsInParts[i] - optimal_num;
2091 extraExchange += extraInPart[i];
2093 if (extraExchange != 0){
2095 if (extraExchange < 0) addition = 1;
2096 for (ssize_t i = 0; i < extraExchange; ++i){
2097 extraInPart[i % this->worldSize] += addition;
2105 int overloadedPartCount = 0;
2106 int *overloadedPartIndices =
new int [this->worldSize];
2109 int underloadedPartCount = 0;
2110 int *underloadedPartIndices =
new int [this->worldSize];
2112 for (
int i = 0; i < this->worldSize; ++i){
2113 if(extraInPart[i] > 0){
2114 overloadedPartIndices[overloadedPartCount++] = i;
2116 else if(extraInPart[i] < 0){
2117 underloadedPartIndices[underloadedPartCount++] = i;
2121 int underloadpartindex = underloadedPartCount - 1;
2124 int numPartsISendFrom = 0;
2125 int *mySendFromParts =
new int[this->worldSize * 2];
2126 gno_t *mySendFromPartsCounts =
new gno_t[this->worldSize * 2];
2128 int numPartsISendTo = 0;
2129 int *mySendParts =
new int[this->worldSize * 2];
2130 gno_t *mySendCountsToParts =
new gno_t[this->worldSize * 2];
2139 for (
int i = overloadedPartCount - 1; i >= 0; --i){
2142 int overloadPart = overloadedPartIndices[i];
2143 gno_t overload = extraInPart[overloadPart];
2144 gno_t myload = numPointsInParts[overloadPart];
2148 gno_t inclusiveLoadUpToMe = numPointsInPartsInclusiveUptoMyIndex[overloadPart];
2151 gno_t exclusiveLoadUptoMe = inclusiveLoadUpToMe - myload;
2154 if (exclusiveLoadUptoMe >= overload){
2158 overloadedPartIndices[i] = -1;
2159 extraInPart[overloadPart] = 0;
2161 while (overload > 0){
2162 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2163 ssize_t underload = extraInPart[nextUnderLoadedPart];
2164 ssize_t left = overload + underload;
2167 extraInPart[nextUnderLoadedPart] = 0;
2168 underloadedPartIndices[underloadpartindex--] = -1;
2171 extraInPart[nextUnderLoadedPart] = left;
2176 else if (exclusiveLoadUptoMe < overload){
2179 gno_t mySendCount = overload - exclusiveLoadUptoMe;
2181 gno_t sendAfterMe = 0;
2184 if (mySendCount > myload){
2185 sendAfterMe = mySendCount - myload;
2186 mySendCount = myload;
2192 mySendFromParts[numPartsISendFrom] = overloadPart;
2193 mySendFromPartsCounts[numPartsISendFrom++] = mySendCount;
2196 while (exclusiveLoadUptoMe > 0){
2197 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2198 ssize_t underload = extraInPart[nextUnderLoadedPart];
2199 ssize_t left = exclusiveLoadUptoMe + underload;
2202 extraInPart[nextUnderLoadedPart] = 0;
2203 underloadedPartIndices[underloadpartindex--] = -1;
2206 extraInPart[nextUnderLoadedPart] = left;
2208 exclusiveLoadUptoMe = left;
2212 while (mySendCount > 0){
2213 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2214 ssize_t underload = extraInPart[nextUnderLoadedPart];
2215 ssize_t left = mySendCount + underload;
2218 mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2219 mySendCountsToParts[numPartsISendTo++] = -underload;
2221 extraInPart[nextUnderLoadedPart] = 0;
2222 underloadedPartIndices[underloadpartindex--] = -1;
2225 extraInPart[nextUnderLoadedPart] = left;
2227 mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2228 mySendCountsToParts[numPartsISendTo++] = mySendCount;
2234 while (sendAfterMe > 0){
2235 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2236 ssize_t underload = extraInPart[nextUnderLoadedPart];
2237 ssize_t left = sendAfterMe + underload;
2240 extraInPart[nextUnderLoadedPart] = 0;
2241 underloadedPartIndices[underloadpartindex--] = -1;
2244 extraInPart[nextUnderLoadedPart] = left;
2255 for (
int i = 0 ; i < numPartsISendFrom; ++i){
2258 int sendFromPart = mySendFromParts[i];
2259 ssize_t sendCount = mySendFromPartsCounts[i];
2260 while(sendCount > 0){
2261 int partToSendIndex = numPartsISendTo - 1;
2262 int partToSend = mySendParts[partToSendIndex];
2264 int sendCountToThisPart = mySendCountsToParts[partToSendIndex];
2268 if (sendCountToThisPart <= sendCount){
2269 mySendParts[partToSendIndex] = 0;
2270 mySendCountsToParts[partToSendIndex] = 0;
2272 sendCount -= sendCountToThisPart;
2275 mySendCountsToParts[partToSendIndex] = sendCountToThisPart - sendCount;
2276 sendCountToThisPart = sendCount;
2281 gno_t toChange = partBegins[sendFromPart];
2282 gno_t previous_begin = partBegins[partToSend];
2285 for (
int k = 0; k < sendCountToThisPart - 1; ++k){
2286 coordinate_grid_parts[toChange] = partToSend;
2287 toChange = partNext[toChange];
2289 coordinate_grid_parts[toChange] = partToSend;
2291 gno_t newBegin = partNext[toChange];
2292 partNext[toChange] = previous_begin;
2293 partBegins[partToSend] = partBegins[sendFromPart];
2294 partBegins[sendFromPart] = newBegin;
2303 for (
int i = 0; i < this->worldSize; ++i) numPointsInParts[i] = 0;
2305 for (
int i = 0; i < this->numLocalCoords; ++i){
2306 numPointsInParts[coordinate_grid_parts[i]] += 1;
2309 reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2312 numGlobalPointsInParts);
2313 if (this->myRank == 0){
2314 std::cout <<
"reassigning" << std::endl;
2315 gno_t totalSize = 0;
2316 for (
int i = 0; i < this->worldSize; ++i){
2317 std::cout <<
"me:" << this->myRank <<
" NumPoints in part:" << i <<
" is: " << numGlobalPointsInParts[i] << std::endl;
2318 totalSize += numGlobalPointsInParts[i];
2320 std::cout <<
"Total:" << totalSize <<
" ng:" << this->numGlobalCoords << std::endl;
2323 delete []mySendCountsToParts;
2324 delete []mySendParts;
2325 delete []mySendFromPartsCounts;
2326 delete []mySendFromParts;
2327 delete []underloadedPartIndices;
2328 delete []overloadedPartIndices;
2329 delete []extraInPart;
2331 delete []partBegins;
2332 delete []numPointsInPartsInclusiveUptoMyIndex;
2333 delete []numPointsInParts;
2334 delete []numGlobalPointsInParts;
2336 delete []shiftProcCount;
2337 delete []bestDimProcs;
2338 delete []dim_slices;
2347 Tpetra::Distributor distributor(comm);
2348 ArrayView<const int> pIds( coordinate_grid_parts, this->numLocalCoords);
2349 gno_t numMyNewGnos = distributor.createFromSends(pIds);
2352 Kokkos::View<scalar_t*, Kokkos::HostSpace> recvBuf2(
2353 Kokkos::ViewAllocateWithoutInitializing(
"recvBuf2"),
2356 for (
int i = 0; i < this->coordinate_dimension; ++i){
2357 Kokkos::View<scalar_t*, Kokkos::HostSpace> s;
2358 if (this->numLocalCoords > 0)
2359 s = Kokkos::View<scalar_t *, Kokkos::HostSpace>(
2360 this->coords[i], this->numLocalCoords);
2362 distributor.doPostsAndWaits(s, 1, recvBuf2);
2364 delete [] this->coords[i];
2365 this->coords[i] =
new scalar_t[numMyNewGnos];
2366 for (
lno_t j = 0; j < numMyNewGnos; ++j){
2367 this->coords[i][j] = recvBuf2[j];
2371 this->numLocalCoords = numMyNewGnos;
2376 int coord_dim = this->coordinate_dimension;
2378 lno_t numLocalPoints = this->numLocalCoords;
2379 gno_t numGlobalPoints = this->numGlobalCoords;
2384 RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
2385 new Tpetra::Map<lno_t, gno_t, node_t> (numGlobalPoints, numLocalPoints, 0, comm));
2387 Teuchos::Array<Teuchos::ArrayView<const scalar_t> > coordView(coord_dim);
2391 for (
int i=0; i < coord_dim; i++){
2392 if(numLocalPoints > 0){
2393 Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalPoints);
2396 Teuchos::ArrayView<const scalar_t> a;
2401 RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmVector = RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >(
2402 new Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t>( mp, coordView.view(0, coord_dim), coord_dim));
2405 RCP<const tMVector_t> coordsConst = Teuchos::rcp_const_cast<
const tMVector_t>(tmVector);
2406 std::vector<const scalar_t *>
weights;
2407 std::vector <int> stride;
2411 inputAdapter_t ia(coordsConst,weights, stride);
2413 Teuchos::RCP <Teuchos::ParameterList> params ;
2414 params =RCP <Teuchos::ParameterList> (
new Teuchos::ParameterList,
true);
2417 params->set(
"algorithm",
"multijagged");
2418 params->set(
"num_global_parts", this->worldSize);
2430 #ifdef HAVE_ZOLTAN2_MPI
2446 for (
lno_t i = 0; i < this->numLocalCoords;++i){
2447 coordinate_grid_parts[i] = partIds[i];
2456 int rank = this->myRank;
2457 int nprocs = this->worldSize;
2460 MEMORY_CHECK(rank==0 || rank==nprocs-1,
"After initializing MPI");
2465 string memoryOn(
"memoryOn");
2466 string memoryOff(
"memoryOff");
2467 bool doMemory=
false;
2468 int numGlobalParts = nprocs;
2472 string balanceCount(
"balance_object_count");
2473 string balanceWeight(
"balance_object_weight");
2474 string mcnorm1(
"multicriteria_minimize_total_weight");
2475 string mcnorm2(
"multicriteria_balance_total_maximum");
2476 string mcnorm3(
"multicriteria_minimize_maximum_weight");
2477 string objective(balanceWeight);
2480 CommandLineProcessor commandLine(
false,
true);
2483 int input_option = 0;
2484 commandLine.setOption(
"input_option", &input_option,
2485 "whether to use mesh creation, geometric generator, or file input");
2486 string inputFile =
"";
2488 commandLine.setOption(
"input_file", &inputFile,
2489 "the input file for geometric generator or file input");
2492 commandLine.setOption(
"size", &numGlobalCoords,
2493 "Approximate number of global coordinates.");
2494 commandLine.setOption(
"numParts", &numGlobalParts,
2495 "Number of parts (default is one per proc).");
2496 commandLine.setOption(
"nWeights", &nWeights,
2497 "Number of weights per coordinate, zero implies uniform weights.");
2498 commandLine.setOption(
"debug", &debugLevel,
"Zoltan1 debug level");
2499 commandLine.setOption(
"remap",
"no-remap", &remap,
2500 "Zoltan1 REMAP parameter; disabled by default for scalability testing");
2501 commandLine.setOption(
"timers", &dummyTimer,
"ignored");
2502 commandLine.setOption(memoryOn.c_str(), memoryOff.c_str(), &doMemory,
2503 "do memory profiling");
2505 string doc(balanceCount);
2506 doc.append(
": ignore weights\n");
2507 doc.append(balanceWeight);
2508 doc.append(
": balance on first weight\n");
2509 doc.append(mcnorm1);
2510 doc.append(
": given multiple weights, balance their total.\n");
2511 doc.append(mcnorm3);
2512 doc.append(
": given multiple weights, "
2513 "balance the maximum for each coordinate.\n");
2514 doc.append(mcnorm2);
2515 doc.append(
": given multiple weights, balance the L2 norm of the weights.\n");
2516 commandLine.setOption(
"objective", &objective, doc.c_str());
2518 CommandLineProcessor::EParseCommandLineReturn rc =
2519 commandLine.parse(0, NULL);
2523 if (rc != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
2524 if (rc == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) {
2525 if (rank==0) std::cout <<
"PASS" << std::endl;
2529 if (rank==0) std::cout <<
"FAIL" << std::endl;
2538 int coord_dim = this->coordinate_dimension;
2541 RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
2542 new Tpetra::Map<lno_t, gno_t, node_t> (this->numGlobalCoords, this->numLocalCoords, 0, this->comm));
2544 Teuchos::Array<Teuchos::ArrayView<const scalar_t> > coordView(coord_dim);
2545 for (
int i=0; i < coord_dim; i++){
2546 if(numLocalCoords > 0){
2547 Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalCoords);
2550 Teuchos::ArrayView<const scalar_t> a;
2555 tMVector_t *tmVector =
new tMVector_t( mp, coordView.view(0, coord_dim), coord_dim);
2558 dots_.
weights.resize(nWeights);
2561 MEMORY_CHECK(doMemory && rank==0,
"After creating input");
2566 int aok = Zoltan_Initialize(0,NULL, &ver);
2569 printf(
"Zoltan_Initialize failed\n");
2573 struct Zoltan_Struct *zz;
2574 zz = Zoltan_Create(MPI_COMM_WORLD);
2576 Zoltan_Set_Param(zz,
"LB_METHOD",
"RCB");
2577 Zoltan_Set_Param(zz,
"LB_APPROACH",
"PARTITION");
2578 Zoltan_Set_Param(zz,
"CHECK_GEOM",
"0");
2579 Zoltan_Set_Param(zz,
"NUM_GID_ENTRIES",
"1");
2580 Zoltan_Set_Param(zz,
"NUM_LID_ENTRIES",
"0");
2581 Zoltan_Set_Param(zz,
"RETURN_LISTS",
"PART");
2582 std::ostringstream oss;
2583 oss << numGlobalParts;
2584 Zoltan_Set_Param(zz,
"NUM_GLOBAL_PARTS", oss.str().c_str());
2587 Zoltan_Set_Param(zz,
"DEBUG_LEVEL", oss.str().c_str());
2590 Zoltan_Set_Param(zz,
"REMAP",
"1");
2592 Zoltan_Set_Param(zz,
"REMAP",
"0");
2594 if (objective != balanceCount){
2597 Zoltan_Set_Param(zz,
"OBJ_WEIGHT_DIM", oss.str().c_str());
2599 if (objective == mcnorm1)
2600 Zoltan_Set_Param(zz,
"RCB_MULTICRITERIA_NORM",
"1");
2601 else if (objective == mcnorm2)
2602 Zoltan_Set_Param(zz,
"RCB_MULTICRITERIA_NORM",
"2");
2603 else if (objective == mcnorm3)
2604 Zoltan_Set_Param(zz,
"RCB_MULTICRITERIA_NORM",
"3");
2607 Zoltan_Set_Param(zz,
"OBJ_WEIGHT_DIM",
"0");
2610 Zoltan_Set_Num_Obj_Fn(zz, getNumObj<tMVector_t>, &dots_);
2611 Zoltan_Set_Obj_List_Fn(zz, getObjList<tMVector_t>, &dots_);
2612 Zoltan_Set_Num_Geom_Fn(zz, getDim<tMVector_t>, &dots_);
2613 Zoltan_Set_Geom_Multi_Fn(zz, getCoords<tMVector_t>, &dots_);
2615 int changes, numGidEntries, numLidEntries, numImport, numExport;
2616 ZOLTAN_ID_PTR importGlobalGids, importLocalGids;
2617 ZOLTAN_ID_PTR exportGlobalGids, exportLocalGids;
2618 int *importProcs, *importToPart, *exportProcs, *exportToPart;
2620 MEMORY_CHECK(doMemory && rank==0,
"Before Zoltan_LB_Partition");
2623 aok = Zoltan_LB_Partition(zz, &changes, &numGidEntries, &numLidEntries,
2624 &numImport, &importGlobalGids, &importLocalGids,
2625 &importProcs, &importToPart,
2626 &numExport, &exportGlobalGids, &exportLocalGids,
2627 &exportProcs, &exportToPart);
2630 MEMORY_CHECK(doMemory && rank==0,
"After Zoltan_LB_Partition");
2632 for (
lno_t i = 0; i < numLocalCoords; i++)
2633 coordinate_grid_parts[i] = exportToPart[i];
2634 Zoltan_Destroy(&zz);
2635 MEMORY_CHECK(doMemory && rank==0,
"After Zoltan_Destroy");
2641 int *coordinate_grid_parts =
new int[this->numLocalCoords];
2642 switch (this->predistribution){
2657 delete []coordinate_grid_parts;
2668 return this->numWeightsPerCoord;
2671 return this->coordinate_dimension;
2674 return this->numLocalCoords;
2677 return this->numGlobalCoords;
2681 return this->coords;
2689 for(
int ii = 0; ii < this->coordinate_dimension; ++ii){
2690 #ifdef HAVE_ZOLTAN2_OMP
2691 #pragma omp parallel for
2693 for (
lno_t i = 0; i < this->numLocalCoords; ++i){
2694 c[ii][i] = this->coords[ii][i];
2700 for(
int ii = 0; ii < this->numWeightsPerCoord; ++ii){
2701 #ifdef HAVE_ZOLTAN2_OMP
2702 #pragma omp parallel for
2704 for (
lno_t i = 0; i < this->numLocalCoords; ++i){
2705 w[ii][i] = this->wghts[ii][i];
virtual bool isInArea(CoordinatePoint< T > dot)=0
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 weighttype getWeight(CoordinatePoint< T > P)=0
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)
map_t::global_ordinal_type gno_t
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_)
SparseMatrixAdapter_t::part_t part_t
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)
map_t::local_ordinal_type lno_t
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.
CoordinatePoint< T > center
virtual CoordinatePoint< T > getPoint(gno_t point_index, unsigned int &state)=0
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_)
Tpetra::MultiVector< zscalar_t, zlno_t, zgno_t, znode_t > tMVector_t
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)