Zoltan2
GeometricGenerator.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 #ifndef GEOMETRICGENERATOR
46 #define GEOMETRICGENERATOR
47 
48 #include <Teuchos_Comm.hpp>
49 #include <Teuchos_ParameterList.hpp>
50 #include <Teuchos_FilteredIterator.hpp>
51 #include <Teuchos_ParameterEntry.hpp>
52 #include <iostream>
53 #include <ctime>
54 #include <limits>
55 #include <climits>
56 #include <string>
57 #include <cstdlib>
58 #include <sstream>
59 #include <fstream>
60 #include <Tpetra_MultiVector_decl.hpp>
63 #include <Teuchos_ArrayViewDecl.hpp>
64 #include <Teuchos_RCP.hpp>
65 #include <Tpetra_Distributor.hpp>
67 
68 
69 #include <zoltan.h>
70 
71 #ifdef _MSC_VER
72 #define NOMINMAX
73 #include <windows.h>
74 #endif
75 
76 using Teuchos::CommandLineProcessor;
77 using Teuchos::ArrayView;
78 using std::string;
79 using Teuchos::ArrayRCP;
80 
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; \
86  return -1; \
87  } \
88  catch (std::logic_error &e) { \
89  std::cout << "Logic exception returned from " << pp << ": " \
90  << e.what() << " FAIL" << std::endl; \
91  return -1; \
92  } \
93  catch (std::bad_alloc &e) { \
94  std::cout << "Bad_alloc exception returned from " << pp << ": " \
95  << e.what() << " FAIL" << std::endl; \
96  return -1; \
97  } \
98  catch (std::exception &e) { \
99  std::cout << "Unknown exception returned from " << pp << ": " \
100  << e.what() << " FAIL" << std::endl; \
101  return -1; \
102  }
103 
104 
105 
106 
107 template <typename tMVector_t>
108 class DOTS{
109 public:
110  std::vector<std::vector<float> > weights;
112 };
113 
114 template <typename tMVector_t>
115 int getNumObj(void *data, int *ierr)
116 {
117  *ierr = 0;
118  DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
119  return dots_->coordinates->getLocalLength();
120 }
122 template <typename tMVector_t>
123 void getCoords(void *data, int numGid, int numLid,
124  int numObj, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids,
125  int dim, double *coords_, int *ierr)
126 {
127  typedef typename tMVector_t::scalar_type scalar_t;
128 
129  // I know that Zoltan asks for coordinates in gid order.
130  if (dim == 3){
131  *ierr = 0;
132  DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
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]);
141  }
142  }
143  else {
144  *ierr = 0;
145  DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
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]);
152  }
153  }
154 }
155 
156 template <typename tMVector_t>
157 int getDim(void *data, int *ierr)
158 {
159  *ierr = 0;
160  DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
161  int dim = dots_->coordinates->getNumVectors();
162 
163  return dim;
164 }
165 
167 template <typename tMVector_t>
168 void getObjList(void *data, int numGid, int numLid,
169  ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids,
170  int num_wgts, float *obj_wgts, int *ierr)
171 {
172  typedef typename tMVector_t::global_ordinal_type gno_t;
173 
174  *ierr = 0;
175  DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
176 
177  size_t localLen = dots_->coordinates->getLocalLength();
178  const gno_t *ids =
179  dots_->coordinates->getMap()->getNodeElementList().getRawPtr();
180 
181  if (sizeof(ZOLTAN_ID_TYPE) == sizeof(gno_t))
182  memcpy(gids, ids, sizeof(ZOLTAN_ID_TYPE) * localLen);
183  else
184  for (size_t i=0; i < localLen; i++)
185  gids[i] = static_cast<ZOLTAN_ID_TYPE>(ids[i]);
186 
187  if (num_wgts > 0){
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];
192  }
193 }
194 
195 
197 const std::string shapes[] = {"SQUARE", "RECTANGLE", "CIRCLE", "CUBE", "RECTANGULAR_PRISM", "SPHERE"};
198 #define SHAPE_COUNT 6
199 
201 const std::string distribution[] = {"distribution", "uniform"};
202 #define DISTRIBUTION_COUNT 2
203 
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"
208 
209 #define INVALID_SHAPE_ARG(SHAPE, REQUIRED) "Invalid argument count for shape " + SHAPE + ". Requires " + REQUIRED + " argument(s)."
210 #define MAX_ITER_ALLOWED 500
211 
212 const std::string weight_distribution_string = "WeightDistribution-";
213 
214 template <typename T>
216  T x;
217  T y;
218  T z;
219  /*
220  bool isInArea(int dim, T *minCoords, T *maxCoords){
221  dim = min(3, dim);
222  for (int i = 0; i < dim; ++i){
223  bool maybe = true;
224  switch(i){
225  case 0:
226  maybe = x < maxCoords[i] && x > maxCoords[i];
227  break;
228  case 1:
229  maybe = y < maxCoords[i] && y > maxCoords[i];
230  break;
231  case 2:
232  maybe = z < maxCoords[i] && z > maxCoords[i];
233  break;
234  }
235  if (!maybe){
236  return false;
237  }
238  }
239  return true;
240  }
241  */
243  x = 0;y=0;z=0;
244  }
245 };
246 
247 template <typename T>
248 class Hole{
249 public:
251  T edge1, edge2, edge3;
252  Hole(CoordinatePoint<T> center_, T edge1_, T edge2_, T edge3_){
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_;
259  }
260  virtual bool isInArea(CoordinatePoint<T> dot) = 0;
261  virtual ~Hole(){}
262 };
263 
264 template <typename T>
265 class SquareHole:public Hole<T>{
266 public:
267  SquareHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
268  virtual ~SquareHole(){}
269  virtual bool isInArea(CoordinatePoint<T> dot){
270  return fabs(dot.x - this->center.x) < this->edge1 / 2 && fabs(dot.y - this->center.y) < this->edge1 / 2;
271  }
272 };
273 
274 template <typename T>
275 class RectangleHole:public Hole<T>{
276 public:
277  RectangleHole(CoordinatePoint<T> center_ , T edge_x_, T edge_y_): Hole<T>(center_, edge_x_, edge_y_, 0){}
278  virtual bool isInArea(CoordinatePoint<T> dot){
279  return fabs(dot.x - this->center.x) < this->edge1 / 2 && fabs(dot.y - this->center.y) < this->edge2 / 2;
280  }
281  virtual ~RectangleHole(){}
282 };
283 
284 template <typename T>
285 class CircleHole:public Hole<T>{
286 public:
287  CircleHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
288  virtual ~CircleHole(){}
289  virtual bool isInArea(CoordinatePoint<T> dot){
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;
291  }
292 };
293 
294 template <typename T>
295 class CubeHole:public Hole<T>{
296 public:
297  CubeHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
298  virtual ~CubeHole(){}
299  virtual bool isInArea(CoordinatePoint<T> dot){
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;
301  }
302 };
303 
304 template <typename T>
305 class RectangularPrismHole:public Hole<T>{
306 public:
307  RectangularPrismHole(CoordinatePoint<T> center_ , T edge_x_, T edge_y_, T edge_z_): Hole<T>(center_, edge_x_, edge_y_, edge_z_){}
309  virtual bool isInArea(CoordinatePoint<T> dot){
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;
311  }
312 };
313 
314 template <typename T>
315 class SphereHole:public Hole<T>{
316 public:
317  SphereHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
318  virtual ~SphereHole(){}
319  virtual bool isInArea(CoordinatePoint<T> dot){
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))
323  <
324  this->edge1 * this->edge1 * this->edge1;
325  }
326 };
327 
328 template <typename T, typename weighttype>
330 public:
331  virtual weighttype getWeight(CoordinatePoint<T> P)=0;
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;
336  virtual ~WeightDistribution(){};
337 };
338 
357 template <typename T, typename weighttype>
358 class SteppedEquation:public WeightDistribution<T, weighttype>{
359  T a1,a2,a3;
360  T b1,b2,b3;
361  T c;
362  T x1,y1,z1;
363 
364  T *steps;
365  weighttype *values;
366  int stepCount;
367 public:
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>(){
369  this->a1 = a1_;
370  this->a2 = a2_;
371  this->a3 = a3_;
372  this->b1 = b1_;
373  this->b2 = b2_;
374  this->b3 = b3_;
375  this->c = c_;
376  this->x1 = x1_;
377  this->y1 = y1_;
378  this->z1 = z1_;
379 
380 
381  this->stepCount = stepCount_;
382  if(this->stepCount > 0){
383  this->steps = new T[this->stepCount];
384  this->values = new T[this->stepCount + 1];
385 
386  for (int i = 0; i < this->stepCount; ++i){
387  this->steps[i] = steps_[i];
388  this->values[i] = values_[i];
389  }
390  this->values[this->stepCount] = values_[this->stepCount];
391  }
392 
393  }
394 
395  virtual ~SteppedEquation(){
396  if(this->stepCount > 0){
397  delete [] this->steps;
398  delete [] this->values;
399  }
400  }
401 
402 
403  virtual weighttype get1DWeight(T x){
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];
408  }
409  return this->values[this->stepCount];
410  }
411  else {
412  return weighttype(expressionRes);
413  }
414  }
415 
416  virtual weighttype get2DWeight(T x, T y){
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];
421  }
422  return this->values[this->stepCount];
423  }
424  else {
425  return weighttype(expressionRes);
426  }
427  }
428 
429  void print (T x, T y, T z){
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;
433 
434  }
435 
436  virtual weighttype get3DWeight(T x, T y, T z){
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;
438 
439  //this->print(x,y,z);
440  if(this->stepCount > 0){
441  for (int i = 0; i < this->stepCount; ++i){
442  if (expressionRes < this->steps[i]) {
443  //std::cout << "0exp:" << expressionRes << " step:" << steps[i] << " value:" << values[i] << std::endl;
444  return this->values[i];
445  }
446  }
447 
448  //std::cout << "1exp:" << expressionRes << " step:" << steps[stepCount] << " value:" << values[stepCount] << std::endl;
449  return this->values[this->stepCount];
450  }
451  else {
452  return weighttype(expressionRes);
453  }
454  }
455 
456  virtual weighttype getWeight(CoordinatePoint<T> p){
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];
461  }
462  return this->values[this->stepCount];
463  }
464  else {
465  return weighttype(expressionRes);
466  }
467  }
468 };
469 
470 template <typename T, typename lno_t, typename gno_t>
472 public:
473  gno_t numPoints;
475  lno_t requested;
479 
480  CoordinateDistribution(gno_t np_, int dim, int wSize) :
481  numPoints(np_), dimension(dim), requested(0), assignedPrevious(0),
482  worldSize(wSize){}
483 
484  virtual CoordinatePoint<T> getPoint(gno_t point_index, unsigned int &state) = 0;
485  virtual T getXCenter() = 0;
486  virtual T getXRadius() =0;
487 
488  void GetPoints(lno_t requestedPointcount, CoordinatePoint<T> *points /*preallocated sized numPoints*/,
489  Hole<T> **holes, lno_t holeCount,
490  float *sharedRatios_, int myRank){
491 
492  for (int i = 0; i < myRank; ++i){
493  //std::cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< std::endl;
494  this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
495  if (i < this->numPoints % this->worldSize){
496  this->assignedPrevious += 1;
497  }
498  }
499 
500  this->requested = requestedPointcount;
501 
502  unsigned int slice = UINT_MAX/(this->worldSize);
503  unsigned int stateBegin = myRank * slice;
504 
505 //#ifdef HAVE_ZOLTAN2_OMP
506 //#pragma omp parallel for
507 //#endif
508 #ifdef HAVE_ZOLTAN2_OMP
509 #pragma omp parallel
510 #endif
511  {
512  int me = 0;
513  int tsize = 1;
514 #ifdef HAVE_ZOLTAN2_OMP
515  me = omp_get_thread_num();
516  tsize = omp_get_num_threads();
517 #endif
518  unsigned int state = stateBegin + me * slice/(tsize);
519 
520 #ifdef HAVE_ZOLTAN2_OMP
521 #pragma omp for
522 #endif
523  for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
524  lno_t iteration = 0;
525  while(1){
526  if(++iteration > MAX_ITER_ALLOWED) {
527  throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
528  }
529  CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, state);
530 
531  bool isInHole = false;
532  for(lno_t i = 0; i < holeCount; ++i){
533  if(holes[i][0].isInArea(p)){
534  isInHole = true;
535  break;
536  }
537  }
538  if(isInHole) continue;
539  points[cnt].x = p.x;
540 
541  points[cnt].y = p.y;
542  points[cnt].z = p.z;
543  break;
544  }
545  }
546  }
547 //#pragma omp parallel
548  /*
549  {
550 
551  lno_t cnt = 0;
552  lno_t iteration = 0;
553  while (cnt < requestedPointcount){
554  if(++iteration > MAX_ITER_ALLOWED) {
555  throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
556  }
557  CoordinatePoint <T> p = this->getPoint();
558 
559  bool isInHole = false;
560  for(lno_t i = 0; i < holeCount; ++i){
561  if(holes[i][0].isInArea(p)){
562  isInHole = true;
563  break;
564  }
565  }
566  if(isInHole) continue;
567  iteration = 0;
568  points[cnt].x = p.x;
569  points[cnt].y = p.y;
570  points[cnt].z = p.z;
571  ++cnt;
572  }
573  }
574  */
575  }
576 
577  void GetPoints(lno_t requestedPointcount, T **coords/*preallocated sized numPoints*/, lno_t tindex,
578  Hole<T> **holes, lno_t holeCount,
579  float *sharedRatios_, int myRank){
580 
581  for (int i = 0; i < myRank; ++i){
582  //std::cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< std::endl;
583  this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
584  if (gno_t(i) < this->numPoints % this->worldSize){
585  this->assignedPrevious += 1;
586  }
587  }
588 
589  this->requested = requestedPointcount;
590 
591  unsigned int slice = UINT_MAX/(this->worldSize);
592  unsigned int stateBegin = myRank * slice;
593 #ifdef HAVE_ZOLTAN2_OMP
594 #pragma omp parallel
595 #endif
596  {
597  int me = 0;
598  int tsize = 1;
599 #ifdef HAVE_ZOLTAN2_OMP
600  me = omp_get_thread_num();
601  tsize = omp_get_num_threads();
602 #endif
603  unsigned int state = stateBegin + me * (slice/(tsize));
604  /*
605 #pragma omp critical
606  {
607 
608  std::cout << "myRank:" << me << " stateBeg:" << stateBegin << " tsize:" << tsize << " state:" << state << " slice: " << slice / tsize << std::endl;
609  }
610  */
611 #ifdef HAVE_ZOLTAN2_OMP
612 #pragma omp for
613 #endif
614  for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
615  lno_t iteration = 0;
616  while(1){
617  if(++iteration > MAX_ITER_ALLOWED) {
618  throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
619  }
620  CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, state);
621 
622  bool isInHole = false;
623  for(lno_t i = 0; i < holeCount; ++i){
624  if(holes[i][0].isInArea(p)){
625  isInHole = true;
626  break;
627  }
628  }
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;
635  }
636  }
637  break;
638  }
639  }
640  }
641  }
642 };
643 
644 template <typename T, typename lno_t, typename gno_t>
646 public:
651 
652 
653  virtual T getXCenter(){
654  return center.x;
655  }
656  virtual T getXRadius(){
657  return standartDevx;
658  }
659 
660  CoordinateNormalDistribution(gno_t np_, int dim, CoordinatePoint<T> center_ ,
661  T sd_x, T sd_y, T sd_z, int wSize) :
662  CoordinateDistribution<T,lno_t,gno_t>(np_,dim,wSize),
663  standartDevx(sd_x), standartDevy(sd_y), standartDevz(sd_z)
664  {
665  this->center.x = center_.x;
666  this->center.y = center_.y;
667  this->center.z = center_.z;
668  }
669 
670  virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
671 
672  pindex = 0; // not used in normal distribution.
674 
675  for(int i = 0; i < this->dimension; ++i){
676  switch(i){
677  case 0:
678  p.x = normalDist(this->center.x, this->standartDevx, state);
679  break;
680  case 1:
681  p.y = normalDist(this->center.y, this->standartDevy, state);
682  break;
683  case 2:
684  p.z = normalDist(this->center.z, this->standartDevz, state);
685  break;
686  default:
687  throw "unsupported dimension";
688  }
689  }
690  return p;
691  }
692 
694 private:
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;
699  if (!derived) {
700  do {
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);
705 
706  polarsqrt=sqrt(-2.0*log(normalsquared)/normalsquared);
707  storedDerivation=normal1*polarsqrt;
708  derived=true;
709  return normal2*polarsqrt*sd + center_;
710  }
711  else {
712  derived=false;
713  return storedDerivation*sd + center_;
714  }
715  }
716 };
717 
718 template <typename T, typename lno_t, typename gno_t>
720 public:
727 
728 
729  virtual T getXCenter(){
730  return (rightMostx - leftMostx)/2 + leftMostx;
731  }
732  virtual T getXRadius(){
733  return (rightMostx - leftMostx)/2;
734  }
735 
736 
737  CoordinateUniformDistribution(gno_t np_, int dim, T l_x, T r_x, T l_y, T r_y,
738  T l_z, T r_z, int wSize ) :
739  CoordinateDistribution<T,lno_t,gno_t>(np_,dim,wSize),
740  leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y),
741  leftMostz(l_z), rightMostz(r_z){}
742 
744  virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
745 
746 
747  pindex = 0; //not used in uniform dist.
749  for(int i = 0; i < this->dimension; ++i){
750  switch(i){
751  case 0:
752  p.x = uniformDist(this->leftMostx, this->rightMostx, state);
753  break;
754  case 1:
755  p.y = uniformDist(this->leftMosty, this->rightMosty, state);
756  break;
757  case 2:
758  p.z = uniformDist(this->leftMostz, this->rightMostz, state);
759  break;
760  default:
761  throw "unsupported dimension";
762  }
763  }
764  return p;
765  }
766 
767 private:
768 
769  T uniformDist(T a, T b, unsigned int &state)
770  {
771  return (b-a)*(rand_r(&state) / double(RAND_MAX)) + a;
772  }
773 };
774 
775 template <typename T, typename lno_t, typename gno_t>
777 public:
784  gno_t along_X, along_Y, along_Z;
785  //T currentX, currentY, currentZ;
787  int myRank;
788  T xstep, ystep, zstep;
789  gno_t xshift, yshift, zshift;
790 
791  virtual T getXCenter(){
792  return (rightMostx - leftMostx)/2 + leftMostx;
793  }
794  virtual T getXRadius(){
795  return (rightMostx - leftMostx)/2;
796  }
797 
798 
799  CoordinateGridDistribution(gno_t alongX, gno_t alongY, gno_t alongZ, int dim,
800  T l_x, T r_x, T l_y, T r_y, T l_z, T r_z ,
801  int myRank_, int wSize) :
802  CoordinateDistribution<T,lno_t,gno_t>(alongX * alongY * alongZ,dim,wSize),
803  leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y), leftMostz(l_z), rightMostz(r_z), myRank(myRank_){
804  //currentX = leftMostx, currentY = leftMosty, currentZ = leftMostz;
805  this->processCnt = 0;
806  this->along_X = alongX; this->along_Y = alongY; this->along_Z = alongZ;
807 
808  if(this->along_X > 1)
809  this->xstep = (rightMostx - leftMostx) / (alongX - 1);
810  else
811  this->xstep = 1;
812  if(this->along_Y > 1)
813  this->ystep = (rightMosty - leftMosty) / (alongY - 1);
814  else
815  this->ystep = 1;
816  if(this->along_Z > 1)
817  this->zstep = (rightMostz - leftMostz) / (alongZ - 1);
818  else
819  this->zstep = 1;
820  xshift = 0; yshift = 0; zshift = 0;
821  }
822 
824  virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
825  //lno_t before = processCnt + this->assignedPrevious;
826  //std::cout << "before:" << processCnt << " " << this->assignedPrevious << std::endl;
827  //lno_t xshift = 0, yshift = 0, zshift = 0;
828 
829  //lno_t tmp = before % (this->along_X * this->along_Y);
830  //xshift = tmp % this->along_X;
831  //yshift = tmp / this->along_X;
832  //zshift = before / (this->along_X * this->along_Y);
833 
834  state = 0; //not used here
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;
838 
839  //this->xshift = pindex / (along_Z * along_Y);
840  // this->zshift = (pindex % (along_Z * along_Y)) / along_Y;
841  //this->yshift = (pindex % (along_Z * along_Y)) % along_Y;
842 
844  p.x = xshift * this->xstep + leftMostx;
845  p.y = yshift * this->ystep + leftMosty;
846  p.z = zshift * this->zstep + leftMostz;
847 /*
848  ++xshift;
849  if(xshift == this->along_X){
850  ++yshift;
851  xshift = 0;
852  if(yshift == this->along_Y){
853  ++zshift;
854  yshift = 0;
855  }
856  }
857 */
858  /*
859  if(this->processCnt == 0){
860  this->xshift = this->assignedPrevious / (along_Z * along_Y);
861  //this->yshift = (this->assignedPrevious % (along_X * along_Y)) / along_X;
862  this->zshift = (this->assignedPrevious % (along_Z * along_Y)) / along_Y;
863  //this->xshift = (this->assignedPrevious % (along_X * along_Y)) % along_X;
864  this->yshift = (this->assignedPrevious % (along_Z * along_Y)) % along_Y;
865  ++this->processCnt;
866  }
867 
868  CoordinatePoint <T> p;
869  p.x = xshift * this->xstep + leftMostx;
870  p.y = yshift * this->ystep + leftMosty;
871  p.z = zshift * this->zstep + leftMostz;
872 
873  ++yshift;
874  if(yshift == this->along_Y){
875  ++zshift;
876  yshift = 0;
877  if(zshift == this->along_Z){
878  ++xshift;
879  zshift = 0;
880  }
881 
882  }
883  */
884  /*
885  if(this->requested - 1 > this->processCnt){
886  this->processCnt++;
887  } else {
888  std::cout << "req:" << this->requested << " pr:" <<this->processCnt << std::endl;
889  std::cout << "p:" << p.x << " " << p.y << " " << p.z << std::endl ;
890  std::cout << "s:" << xshift << " " << yshift << " " << zshift << std::endl ;
891  std::cout << "st:" << this->xstep << " " << this->ystep << " " << this->zstep << std::endl ;
892  }
893  */
894  return p;
895  }
896 
897 private:
898 
899 };
900 
901 template <typename scalar_t, typename lno_t, typename gno_t, typename node_t>
903 private:
904  Hole<scalar_t> **holes; //to represent if there is any hole in the input
905  int holeCount;
906  int coordinate_dimension; //dimension of the geometry
907  gno_t numGlobalCoords; //global number of coordinates requested to be created.
908  lno_t numLocalCoords;
909  float *loadDistributions; //sized as the number of processors, the load of each processor.
910  bool loadDistSet;
911  bool distinctCoordSet;
912  CoordinateDistribution<scalar_t, lno_t,gno_t> **coordinateDistributions;
913  int distributionCount;
914  //CoordinatePoint<scalar_t> *points;
915  scalar_t **coords;
916  scalar_t **wghts;
918  int numWeightsPerCoord;
919  int predistribution;
920  RCP<const Teuchos::Comm<int> > comm;
921  //RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmVector;
922  //RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmwVector;
923  int worldSize;
924  int myRank;
925  scalar_t minx;
926  scalar_t maxx;
927  scalar_t miny;
928  scalar_t maxy;
929  scalar_t minz;
930  scalar_t maxz;
931  std::string outfile;
932  float perturbation_ratio;
933 
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;
936 
937 
938  template <typename tt>
939  tt getParamVal( const Teuchos::ParameterEntry& pe, const std::string &paramname){
940  tt returnVal;
941  try {
942  returnVal = Teuchos::getValue<tt>(pe);
943  }
944  catch (...){
945  throw INVALID(paramname);
946  }
947  return returnVal;
948  }
949 
950  int countChar (std::string inStr, char cntChar){
951  int cnt = 0;
952  for (unsigned int i = 0; i < inStr.size(); ++i){
953  if (inStr[i] == cntChar) {
954  cnt++;
955  }
956  }
957  return cnt;
958  }
959 
960  template <typename tt>
961  std::string toString(tt obj){
962  std::stringstream ss (std::stringstream::in |std::stringstream::out);
963  ss << obj;
964  std::string tmp = "";
965  ss >> tmp;
966  return tmp;
967  }
968 
969  template <typename tt>
970  tt fromString(std::string obj){
971  std::stringstream ss (std::stringstream::in | std::stringstream::out);
972  ss << obj;
973  tt tmp;
974  ss >> tmp;
975  if (ss.fail()){
976  throw "Cannot convert string " + obj;
977  }
978  return tmp;
979  }
980 
981 
982  void splitString(std::string inStr, char splitChar, std::string *outSplittedStr){
983  std::stringstream ss (std::stringstream::in | std::stringstream::out);
984  ss << inStr;
985  int cnt = 0;
986  while (!ss.eof()){
987  std::string tmp = "";
988  std::getline(ss, tmp ,splitChar);
989  outSplittedStr[cnt++] = tmp;
990  }
991  }
992 
993  /*
994  void getDistinctCoordinateDescription(std::string distinctDescription){
995 
996  this->distinctCoordinates = new bool[this->dimension];
997  if(distinctDescription != ""){
998  int argCnt = this->countChar(distinctDescription, ',') + 1;
999  if(argCnt != this->dimension) {
1000  throw "Invalid parameter count for distinct_coordinates. Size should be equal to dimension.";
1001  }
1002 
1003  std::string *splittedStr = new std::string[argCnt];
1004  splitString(holeDescription, ',', splittedStr);
1005 
1006  for(int i = 0; i < argCnt; ++i){
1007  if(splittedStr[i] == "T"){
1008  distinctCoordinates[i] = true;
1009  } else if(splittedStr[i] == "F"){
1010  distinctCoordinates[i] = false;
1011  } else {
1012  throw "Invalid parameter " + splittedStr[i] + " for distinct_coordinates.";
1013  }
1014  }
1015  delete []splittedStr;
1016  }
1017  else {
1018  for(int i = 0; i < this->dimension; ++i){
1019  distinctCoordinates[i] = false;
1020  }
1021  }
1022  }
1023  */
1024 
1025 
1026  void getCoordinateDistributions(std::string coordinate_distributions){
1027  if(coordinate_distributions == ""){
1028  throw "At least one distribution function must be provided to coordinate generator.";
1029  }
1030 
1031  int argCnt = this->countChar(coordinate_distributions, ',') + 1;
1032  std::string *splittedStr = new std::string[argCnt];
1033  splitString(coordinate_distributions, ',', splittedStr);
1034  coordinateDistributions = (CoordinateDistribution<scalar_t, lno_t,gno_t> **) malloc(sizeof (CoordinateDistribution<scalar_t, lno_t,gno_t> *) * 1);
1035  for(int i = 0; i < argCnt; ){
1036  coordinateDistributions = (CoordinateDistribution<scalar_t, lno_t,gno_t> **)realloc((void *)coordinateDistributions, (this->distributionCount + 1)* sizeof(CoordinateDistribution<scalar_t, lno_t,gno_t> *));
1037 
1038  std::string distName = splittedStr[i++];
1039  gno_t np_ = 0;
1040  if(distName == "NORMAL"){
1041  int reqArg = 5;
1042  if (this->coordinate_dimension == 3){
1043  reqArg = 7;
1044  }
1045  if(i + reqArg > argCnt) {
1046  std::string tmp = toString<int>(reqArg);
1047  throw INVALID_SHAPE_ARG(distName, tmp);
1048  }
1049  np_ = fromString<gno_t>(splittedStr[i++]);
1051 
1052  pp.x = fromString<scalar_t>(splittedStr[i++]);
1053  pp.y = fromString<scalar_t>(splittedStr[i++]);
1054  pp.z = 0;
1055  if(this->coordinate_dimension == 3){
1056  pp.z = fromString<scalar_t>(splittedStr[i++]);
1057  }
1058 
1059  scalar_t sd_x = fromString<scalar_t>(splittedStr[i++]);
1060  scalar_t sd_y = fromString<scalar_t>(splittedStr[i++]);
1061  scalar_t sd_z = 0;
1062  if(this->coordinate_dimension == 3){
1063  sd_z = fromString<scalar_t>(splittedStr[i++]);
1064  }
1065  this->coordinateDistributions[this->distributionCount++] = new CoordinateNormalDistribution<scalar_t, lno_t,gno_t>(np_, this->coordinate_dimension, pp , sd_x, sd_y, sd_z, this->worldSize );
1066 
1067  } else if(distName == "UNIFORM" ){
1068  int reqArg = 5;
1069  if (this->coordinate_dimension == 3){
1070  reqArg = 7;
1071  }
1072  if(i + reqArg > argCnt) {
1073  std::string tmp = toString<int>(reqArg);
1074  throw INVALID_SHAPE_ARG(distName, tmp);
1075  }
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++]);
1081 
1082  scalar_t l_z = 0, r_z = 0;
1083 
1084  if(this->coordinate_dimension == 3){
1085  l_z = fromString<scalar_t>(splittedStr[i++]);
1086  r_z = fromString<scalar_t>(splittedStr[i++]);
1087  }
1088 
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"){
1091  int reqArg = 6;
1092  if(this->coordinate_dimension == 3){
1093  reqArg = 9;
1094  }
1095  if(i + reqArg > argCnt) {
1096  std::string tmp = toString<int>(reqArg);
1097  throw INVALID_SHAPE_ARG(distName, tmp);
1098  }
1099 
1100  gno_t np_x = fromString<gno_t>(splittedStr[i++]);
1101  gno_t np_y = fromString<gno_t>(splittedStr[i++]);
1102  gno_t np_z = 1;
1103 
1104 
1105  if(this->coordinate_dimension == 3){
1106  np_z = fromString<gno_t>(splittedStr[i++]);
1107  }
1108 
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++]);
1114 
1115  scalar_t l_z = 0, r_z = 0;
1116 
1117  if(this->coordinate_dimension == 3){
1118  l_z = fromString<scalar_t>(splittedStr[i++]);
1119  r_z = fromString<scalar_t>(splittedStr[i++]);
1120  }
1121 
1122  if(np_x < 1 || np_z < 1 || np_y < 1 ){
1123  throw "Provide at least 1 point along each dimension for grid test.";
1124  }
1125  //std::cout << "ly:" << l_y << " ry:" << r_y << std::endl;
1126  this->coordinateDistributions[this->distributionCount++] = new CoordinateGridDistribution<scalar_t, lno_t,gno_t>
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);
1128 
1129  }
1130  else {
1131  std::string tmp = toString<int>(this->coordinate_dimension);
1132  throw INVALIDSHAPE(distName, tmp);
1133  }
1134  this->numGlobalCoords += (gno_t) np_;
1135  }
1136  delete []splittedStr;
1137  }
1138 
1139  void getProcLoadDistributions(std::string proc_load_distributions){
1140 
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;
1146  }
1147  } else{
1148 
1149 
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);
1153  }
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]);
1158  }
1159  delete []splittedStr;
1160 
1161 
1162  float sum = 0;
1163  for(int i = 0; i < this->worldSize; ++i){
1164  sum += this->loadDistributions[i];
1165  }
1166  if (fabs(sum - 1.0) > 10*std::numeric_limits<float>::epsilon()){
1167  throw "Processor load ratios do not sum to 1.0.";
1168  }
1169  }
1170 
1171  }
1172 
1173  void getHoles(std::string holeDescription){
1174  if(holeDescription == ""){
1175  return;
1176  }
1177  this->holes = (Hole<scalar_t> **) malloc(sizeof (Hole <scalar_t>*));
1178  int argCnt = this->countChar(holeDescription, ',') + 1;
1179  std::string *splittedStr = new std::string[argCnt];
1180  splitString(holeDescription, ',', splittedStr);
1181 
1182  for(int i = 0; i < argCnt; ){
1183  holes = (Hole<scalar_t> **)realloc((void *)holes, (this->holeCount + 1)* sizeof(Hole<scalar_t> *));
1184 
1185  std::string shapeName = splittedStr[i++];
1186  if(shapeName == "SQUARE" && this->coordinate_dimension == 2){
1187  if(i + 3 > argCnt) {
1188  throw INVALID_SHAPE_ARG(shapeName, "3");
1189  }
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++]);
1194  this->holes[this->holeCount++] = new SquareHole<scalar_t>(pp, edge);
1195  } else if(shapeName == "RECTANGLE" && this->coordinate_dimension == 2){
1196  if(i + 4 > argCnt) {
1197  throw INVALID_SHAPE_ARG(shapeName, "4");
1198  }
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++]);
1204 
1205  this->holes[this->holeCount++] = new RectangleHole<scalar_t>(pp, edgex,edgey);
1206  } else if(shapeName == "CIRCLE" && this->coordinate_dimension == 2){
1207  if(i + 3 > argCnt) {
1208  throw INVALID_SHAPE_ARG(shapeName, "3");
1209  }
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++]);
1214  this->holes[this->holeCount++] = new CircleHole<scalar_t>(pp, r);
1215  } else if(shapeName == "CUBE" && this->coordinate_dimension == 3){
1216  if(i + 4 > argCnt) {
1217  throw INVALID_SHAPE_ARG(shapeName, "4");
1218  }
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++]);
1224  this->holes[this->holeCount++] = new CubeHole<scalar_t>(pp, edge);
1225  } else if(shapeName == "RECTANGULAR_PRISM" && this->coordinate_dimension == 3){
1226  if(i + 6 > argCnt) {
1227  throw INVALID_SHAPE_ARG(shapeName, "6");
1228  }
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++]);
1236  this->holes[this->holeCount++] = new RectangularPrismHole<scalar_t>(pp, edgex, edgey, edgez);
1237 
1238  } else if(shapeName == "SPHERE" && this->coordinate_dimension == 3){
1239  if(i + 4 > argCnt) {
1240  throw INVALID_SHAPE_ARG(shapeName, "4");
1241  }
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++]);
1247  this->holes[this->holeCount++] = new SphereHole<scalar_t>(pp, r);
1248  } else {
1249  std::string tmp = toString<int>(this->coordinate_dimension);
1250  throw INVALIDSHAPE(shapeName, tmp);
1251  }
1252  }
1253  delete [] splittedStr;
1254  }
1255 
1256  void getWeightDistribution(std::string *weight_distribution_arr, int wdimension){
1257  int wcount = 0;
1258 
1259  this->wd = new WeightDistribution<scalar_t,scalar_t> *[wdimension];
1260  for(int ii = 0; ii < MAX_WEIGHT_DIM; ++ii){
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.";
1265  }
1266 
1267  int count = this->countChar(weight_distribution, ' ');
1268  std::string *splittedStr = new string[count + 1];
1269  this->splitString(weight_distribution, ' ', splittedStr);
1270  //std::cout << count << std::endl;
1271  scalar_t c=1;
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;
1277  int stepCount = 0;
1278  int valueCount = 1;
1279 
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];
1284  }
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];
1289  //std::cout << "parameter:" << parameter << " value:" << value << std::endl;
1290 
1291  if (parameter == "a1"){
1292  a1 = this->fromString<scalar_t>(value);
1293  }
1294  else if (parameter == "a2"){
1295  if(this->coordinate_dimension > 1){
1296  a2 = this->fromString<scalar_t>(value);
1297  }
1298  else {
1299  throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1300  }
1301 
1302  }
1303  else if (parameter == "a3"){
1304  if(this->coordinate_dimension > 2){
1305  a3 = this->fromString<scalar_t>(value);
1306  }
1307  else {
1308  throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1309  }
1310  }
1311  else if (parameter == "b1"){
1312  b1 = this->fromString<scalar_t>(value);
1313  }
1314  else if (parameter == "b2"){
1315  if(this->coordinate_dimension > 1){
1316  b2 = this->fromString<scalar_t>(value);
1317  }
1318  else {
1319  throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1320  }
1321  }
1322  else if (parameter == "b3"){
1323 
1324  if(this->coordinate_dimension > 2){
1325  b3 = this->fromString<scalar_t>(value);
1326  }
1327  else {
1328  throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1329  }
1330  }
1331  else if (parameter == "c"){
1332  c = this->fromString<scalar_t>(value);
1333  }
1334  else if (parameter == "x1"){
1335  x1 = this->fromString<scalar_t>(value);
1336  }
1337  else if (parameter == "y1"){
1338  if(this->coordinate_dimension > 1){
1339  y1 = this->fromString<scalar_t>(value);
1340  }
1341  else {
1342  throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1343  }
1344  }
1345  else if (parameter == "z1"){
1346  if(this->coordinate_dimension > 2){
1347  z1 = this->fromString<scalar_t>(value);
1348  }
1349  else {
1350  throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
1351  }
1352  }
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]);
1360  }
1361  delete [] stepstr;
1362  }
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]);
1370  }
1371  delete [] stepstr;
1372  }
1373  else {
1374  throw "Invalid parameter name at " + splittedStr[i];
1375  }
1376  }
1377 
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);
1381  }
1382 
1383 
1384  this->wd[wcount++] = new SteppedEquation<scalar_t,scalar_t>(a1, a2, a3, b1, b2, b3, c, x1, y1, z1, steps, values, stepCount);
1385 
1386  if(stepCount > 0){
1387  delete [] steps;
1388  delete [] values;
1389 
1390  }
1391  }
1392  if(wcount != this->numWeightsPerCoord){
1393  throw "Weight Dimension is provided as " + toString<int>(wdimension) + ". But " + toString<int>(wcount)+" weight distributions are provided.";
1394  }
1395  }
1396 
1397  void parseParams(const Teuchos::ParameterList & params){
1398  try {
1399  std::string holeDescription = "";
1400  std::string proc_load_distributions = "";
1401  std::string distinctDescription = "";
1402  std::string coordinate_distributions = "";
1403  std::string numWeightsPerCoord_parameters[MAX_WEIGHT_DIM];
1404  for (int i = 0; i < MAX_WEIGHT_DIM; ++i){
1405  numWeightsPerCoord_parameters[i] = "";
1406  }
1407 
1408 
1409  for (Teuchos::ParameterList::ConstIterator pit = params.begin(); pit != params.end(); ++pit ){
1410  const std::string &paramName = params.name(pit);
1411 
1412  const Teuchos::ParameterEntry &pe = params.entry(pit);
1413 
1414  if(paramName.find("hole-") == 0){
1415  if(holeDescription == ""){
1416  holeDescription = getParamVal<std::string>(pe, paramName);
1417  }
1418  else {
1419  holeDescription +=","+getParamVal<std::string>(pe, paramName);
1420  }
1421  }
1422  else if(paramName.find("distribution-") == 0){
1423  if(coordinate_distributions == "")
1424  coordinate_distributions = getParamVal<std::string>(pe, paramName);
1425  else
1426  coordinate_distributions += ","+getParamVal<std::string>(pe, paramName);
1427  //std::cout << coordinate_distributions << std::endl;
1428  //TODO coordinate distribution description
1429  }
1430 
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);
1436 
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);
1439  }
1440  numWeightsPerCoord_parameters[distribution_index] += " " + weight_dist_param.substr(dash_pos + 1)+ "="+ getParamVal<std::string>(pe, paramName);
1441  }
1442  else if(paramName == "dim"){
1443  int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
1444  if(dim < 2 && dim > 3){
1445  throw INVALID(paramName);
1446  } else {
1447  this->coordinate_dimension = dim;
1448  }
1449  }
1450  else if(paramName == "wdim"){
1451  int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
1452  if(dim < 1 && dim > MAX_WEIGHT_DIM){
1453  throw INVALID(paramName);
1454  } else {
1455  this->numWeightsPerCoord = dim;
1456  }
1457  }
1458  else if(paramName == "predistribution"){
1459  int pre_distribution = fromString<int>(getParamVal<std::string>(pe, paramName));
1460  if(pre_distribution < 0 && pre_distribution > 3){
1461  throw INVALID(paramName);
1462  } else {
1463  this->predistribution = pre_distribution;
1464  }
1465  }
1466  else if(paramName == "perturbation_ratio"){
1467  float perturbation = fromString<float>(getParamVal<std::string>(pe, paramName));
1468  if(perturbation < 0 && perturbation > 1){
1469  throw INVALID(paramName);
1470  } else {
1471  this->perturbation_ratio = perturbation;
1472  }
1473  }
1474 
1475 
1476  else if(paramName == "proc_load_distributions"){
1477  proc_load_distributions = getParamVal<std::string>(pe, paramName);
1478  this->loadDistSet = true;
1479  }
1480 
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;
1487  } else {
1488  throw "Invalid parameter for distinct_coordinates: " + distinctDescription + ". Candidates: T or F." ;
1489  }
1490  }
1491 
1492  else if(paramName == "out_file"){
1493  this->outfile = getParamVal<std::string>(pe, paramName);
1494  }
1495 
1496  else {
1497  throw INVALID(paramName);
1498  }
1499  }
1500 
1501 
1502  if(this->coordinate_dimension == 0){
1503  throw "Dimension must be provided to coordinate generator.";
1504  }
1505  /*
1506  if(this->globalNumCoords == 0){
1507  throw "Number of coordinates must be provided to coordinate generator.";
1508  }
1509  */
1510  /*
1511  if(maxx <= minx ){
1512  throw "Error: maxx= "+ toString<scalar_t>(maxx)+ " and minx=" + toString<scalar_t>(minx);
1513  }
1514  if(maxy <= miny ){
1515  throw "Error: maxy= "+ toString<scalar_t>(maxy)+ " and miny=" + toString<scalar_t>(miny);
1516 
1517  }
1518  if(this->dimension == 3 && maxz <= minz ){
1519  throw "Error: maxz= "+ toString<scalar_t>(maxz)+ " and minz=" + toString<scalar_t>(minz);
1520  }
1521  */
1522  if (this->loadDistSet && this->distinctCoordSet){
1523  throw "distinct_coordinates and proc_load_distributions parameters cannot be satisfied together.";
1524  }
1525  this->getHoles(holeDescription);
1526  //this->getDistinctCoordinateDescription(distinctDescription);
1527  this->getProcLoadDistributions(proc_load_distributions);
1528  this->getCoordinateDistributions(coordinate_distributions);
1529  this->getWeightDistribution(numWeightsPerCoord_parameters, this->numWeightsPerCoord);
1530  /*
1531  if(this->numGlobalCoords <= 0){
1532  throw "Must have at least 1 point";
1533  }
1534  */
1535  }
1536  catch(std::string s){
1537  throw s;
1538  }
1539 
1540  catch(char * s){
1541  throw s;
1542  }
1543 
1544  catch(char const* s){
1545  throw s;
1546  }
1547 
1548  }
1549 public:
1550 
1552  if(this->holes){
1553  for (int i = 0; i < this->holeCount; ++i){
1554  delete this->holes[i];
1555  }
1556  free (this->holes);
1557  }
1558 
1559 
1560  delete []loadDistributions; //sized as the number of processors, the load of each processor.
1561  //delete []distinctCoordinates; // if processors have different or same range for coordinates to be created.
1562  if(coordinateDistributions){
1563 
1564  for (int i = 0; i < this->distributionCount; ++i){
1565  delete this->coordinateDistributions[i];
1566  }
1567  free (this->coordinateDistributions);
1568  }
1569  if (this->wd){
1570  for (int i = 0; i < this->numWeightsPerCoord; ++i){
1571  delete this->wd[i];
1572  }
1573  delete []this->wd;
1574  }
1575 
1576  if(this->numWeightsPerCoord){
1577  for(int i = 0; i < this->numWeightsPerCoord; ++i)
1578  delete [] this->wghts[i];
1579  delete []this->wghts;
1580  }
1581  if(this->coordinate_dimension){
1582  for(int i = 0; i < this->coordinate_dimension; ++i)
1583  delete [] this->coords[i];
1584  delete [] this->coords;
1585  }
1586  //delete []this->points;
1587  }
1588 
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;
1627  }
1628 
1629  GeometricGenerator(Teuchos::ParameterList &params, const RCP<const Teuchos::Comm<int> > & comm_){
1630  this->wd = NULL;
1631  this->comm = comm_;
1632  this->holes = NULL; //to represent if there is any hole in the input
1633  this->coordinate_dimension = 0; //dimension of the geometry
1634  this->numWeightsPerCoord = 0;
1635  this->worldSize = comm_->getSize(); //comminication world object.
1636  this->numGlobalCoords = 0; //global number of coordinates requested to be created.
1637  this->loadDistributions = NULL; //sized as the number of processors, the load of each processor.
1638  //this->distinctCoordinates = NULL; // if processors have different or same range for coordinates to be created.
1639  this->coordinateDistributions = NULL;
1640  this->holeCount = 0;
1641  this->distributionCount = 0;
1642  this->outfile = "";
1643  this->predistribution = 0;
1644  this->perturbation_ratio = 0;
1645  //this->points = NULL;
1646 
1647  /*
1648  this->minx = 0;
1649  this->maxx = 0;
1650  this->miny = 0;
1651  this->maxy = 0;
1652  this->minz = 0;
1653  this->maxz = 0;
1654  */
1655  this->loadDistSet = false;
1656  this->distinctCoordSet = false;
1657  this->myRank = comm_->getRank();
1658 
1659  try {
1660  this->parseParams(params);
1661  }
1662  catch(std::string s){
1663  if(myRank == 0){
1664  print_description();
1665  }
1666  throw s;
1667  }
1668 
1669 
1670 
1671 
1672  lno_t myPointCount = 0;
1673  this->numGlobalCoords = 0;
1674 
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){
1680  increment += 1;
1681  }
1682  this->numGlobalCoords += increment;
1683  if(ii < myRank){
1684  prefixSum += increment;
1685  }
1686  }
1687  myPointCount += lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
1688  if (gno_t(myRank) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1689  myPointCount += 1;
1690  }
1691  }
1692 
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];
1696  }
1697 
1698  for (int ii = 0; ii < this->coordinate_dimension; ++ii){
1699 #ifdef HAVE_ZOLTAN2_OMP
1700 #pragma omp parallel for
1701 #endif
1702  for(lno_t i = 0; i < myPointCount; ++i){
1703  this->coords[ii][i] = 0;
1704  }
1705  }
1706 
1707  this->numLocalCoords = 0;
1708  srand ((myRank + 1) * this->numLocalCoords);
1709  for (int i = 0; i < distributionCount; ++i){
1710 
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;
1714  }
1715  //std::cout << "req:" << requestedPointCount << std::endl;
1716  //this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->points + this->numLocalCoords, this->holes, this->holeCount, this->loadDistributions, myRank);
1717  this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->coords, this->numLocalCoords, this->holes, this->holeCount, this->loadDistributions, myRank);
1718  this->numLocalCoords += requestedPointCount;
1719  }
1720 
1721  /*
1722 
1723  if (this->myRank >= 0){
1724  for(lno_t i = 0; i < this->numLocalCoords; ++i){
1725 
1726  std::cout <<"me:" << this->myRank << " "<< this->coords[0][i];
1727  if(this->coordinate_dimension > 1){
1728  std::cout << " " << this->coords[1][i];
1729  }
1730  if(this->coordinate_dimension > 2){
1731  std::cout << " " << this->coords[2][i];
1732  }
1733  std::cout << std::endl;
1734  }
1735  }
1736  */
1737 
1738 
1739 
1740  if (this->predistribution){
1741  redistribute();
1742  }
1743 
1744 
1745 
1746  int scale = 3;
1747  if (this->perturbation_ratio > 0){
1748  this->perturb_data(this->perturbation_ratio, scale);
1749  }
1750  /*
1751  if (this->myRank >= 0){
1752  std::cout << "after distribution" << std::endl;
1753  for(lno_t i = 0; i < this->numLocalCoords; ++i){
1754 
1755  std::cout <<"me:" << this->myRank << " " << this->coords[0][i];
1756  if(this->coordinate_dimension > 1){
1757  std::cout << " " << this->coords[1][i];
1758  }
1759  if(this->coordinate_dimension > 2){
1760  std::cout << " " << this->coords[2][i];
1761  }
1762  std::cout << std::endl;
1763  }
1764  }
1765 
1766  */
1767 
1768 
1769  if (this->distinctCoordSet){
1770  //TODO: Partition and migration.
1771  }
1772 
1773  if(this->outfile != ""){
1774 
1775  std::ofstream myfile;
1776  myfile.open ((this->outfile + toString<int>(myRank)).c_str());
1777  for(lno_t i = 0; i < this->numLocalCoords; ++i){
1778 
1779  myfile << this->coords[0][i];
1780  if(this->coordinate_dimension > 1){
1781  myfile << " " << this->coords[1][i];
1782  }
1783  if(this->coordinate_dimension > 2){
1784  myfile << " " << this->coords[2][i];
1785  }
1786  myfile << std::endl;
1787  }
1788  myfile.close();
1789 
1790  if (this->myRank == 0){
1791  std::ofstream gnuplotfile("gnu.gnuplot");
1792  for(int i = 0; i < this->worldSize; ++i){
1793  string s = "splot";
1794  if (this->coordinate_dimension == 2){
1795  s = "plot";
1796  }
1797  if (i > 0){
1798  s = "replot";
1799  }
1800  gnuplotfile << s << " \"" << (this->outfile + toString<int>(i)) << "\"" << std::endl;
1801  }
1802  gnuplotfile << "pause -1" << std::endl;
1803  gnuplotfile.close();
1804  }
1805  }
1806 
1807 
1808 
1809  /*
1810  Zoltan2::XpetraMultiVectorAdapter < Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> > xmv (RCP < Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> > (tmVector));
1811 
1812  RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmVector2;
1813  Zoltan2::PartitioningSolution< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> > solution;
1814  xmv.applyPartitioningSolution<Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >(this->tmVector, &tmVector2, solution);
1815  */
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];
1820  }
1821  }
1822 
1823  for(int ii = 0; ii < this->numWeightsPerCoord; ++ii){
1824  switch(this->coordinate_dimension){
1825  case 1:
1826 #ifdef HAVE_ZOLTAN2_OMP
1827 #pragma omp parallel for
1828 #endif
1829  for (lno_t i = 0; i < this->numLocalCoords; ++i){
1830  this->wghts[ii][i] = this->wd[ii]->get1DWeight(this->coords[0][i]);
1831  }
1832  break;
1833  case 2:
1834 #ifdef HAVE_ZOLTAN2_OMP
1835 #pragma omp parallel for
1836 #endif
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]);
1839  }
1840  break;
1841  case 3:
1842 #ifdef HAVE_ZOLTAN2_OMP
1843 #pragma omp parallel for
1844 #endif
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]);
1847  }
1848  break;
1849  }
1850  }
1851  }
1852 
1853  //############################################################//
1855  //############################################################//
1856  void perturb_data(float used_perturbation_ratio, int scale){
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];
1864  }
1865 
1866  if (maxCoords[dim] < this->coords[dim][i]){
1867  maxCoords[dim] = this->coords[dim][i];
1868  }
1869  }
1870 
1871 
1872 
1873 
1874  scalar_t center = (maxCoords[dim] + minCoords[dim]) / 2;
1875 
1876  minCoords[dim] = center - (center - minCoords[dim]) * scale;
1877  maxCoords[dim] = (maxCoords[dim] - center) * scale + center;
1878 
1879  }
1880 
1881  gno_t numLocalToPerturb = gno_t(this->numLocalCoords*used_perturbation_ratio);
1882  //std::cout << "numLocalToPerturb :" << numLocalToPerturb << std::endl;
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];
1887 
1888  }
1889  }
1890  delete []maxCoords;
1891  delete []minCoords;
1892  }
1893 
1894  //############################################################//
1896  //############################################################//
1897 
1898  //Returns the partitioning dimension as even as possible
1899  void getBestSurface (int remaining, int *dimProcs, int dim, int currentDim, int &bestSurface, int *bestDimProcs){
1900 
1901  if (currentDim < dim - 1){
1902  int ipx = 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);
1908  }
1909  ipx++;
1910  }
1911  }
1912  else {
1913  dimProcs[currentDim] = remaining;
1914  int surface = 0;
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];
1919  }
1920  }
1921 
1922  }
1923 
1924  //returns min and max coordinates along each dimension
1925  void getMinMaxCoords(scalar_t *globalMaxCoords, scalar_t *globalMinCoords){
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];
1933  }
1934 
1935  if (maxCoords[dim] < this->coords[dim][i]){
1936  maxCoords[dim] = this->coords[dim][i];
1937  }
1938  }
1939  }
1940 
1941  reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MAX,
1942  this->coordinate_dimension,
1943  maxCoords,
1944  globalMaxCoords);
1945 
1946 
1947  reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MIN,
1948  this->coordinate_dimension,
1949  minCoords,
1950  globalMinCoords);
1951 
1952  delete []maxCoords;
1953  delete []minCoords;
1954  }
1955 
1956 
1957  //performs a block partitioning.
1958  //then distributes the points of the overloaded parts to underloaded parts.
1959  void blockPartition(int *coordinate_grid_parts){
1960 
1961 #ifdef _MSC_VER
1962  typedef SSIZE_T ssize_t;
1963 #endif
1964 
1965  //############################################################//
1967  //############################################################//
1968 
1969  scalar_t *maxCoords= new scalar_t [this->coordinate_dimension];
1970  scalar_t *minCoords= new scalar_t [this->coordinate_dimension];
1971  //global min and max coordinates in each dimension.
1972  this->getMinMaxCoords(maxCoords, minCoords);
1973 
1974 
1975  //############################################################//
1977  //############################################################//
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];
1982  int currentDim = 0;
1983 
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];
1989  //divides the parts into dimensions as even as possible.
1990  getBestSurface ( remaining, dimProcs, coord_dim, currentDim, bestSurface, bestDimProcs);
1991 
1992 
1993  delete []dimProcs;
1994 
1995  //############################################################//
1997  //############################################################//
1998  int *shiftProcCount = new int[coord_dim];
1999  //how the consecutive parts along a dimension
2000  //differs in the index.
2001 
2002  int remainingProc = this->worldSize;
2003  for (int dim = 0; dim < coord_dim; ++dim){
2004  remainingProc = remainingProc / bestDimProcs[dim];
2005  shiftProcCount[dim] = remainingProc;
2006  }
2007 
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;
2013  }
2014 
2015  //############################################################//
2017  //############################################################//
2018 
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];
2022 
2023  gno_t *partBegins = new gno_t [this->worldSize];
2024  gno_t *partNext = new gno_t[this->numLocalCoords];
2025 
2026 
2027  for (lno_t i = 0; i < this->numLocalCoords; ++i){
2028  partNext[i] = -1;
2029  }
2030  for (int i = 0; i < this->worldSize; ++i) {
2031  partBegins[i] = 1;
2032  }
2033 
2034  for (int i = 0; i < this->worldSize; ++i)
2035  numPointsInParts[i] = 0;
2036 
2037  for (lno_t i = 0; i < this->numLocalCoords; ++i){
2038  int partIndex = 0;
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;
2043  }
2044  shift = shift * shiftProcCount[dim];
2045  partIndex += shift;
2046  }
2047  numPointsInParts[partIndex] += 1;
2048  coordinate_grid_parts[i] = partIndex;
2049 
2050  partNext[i] = partBegins[partIndex];
2051  partBegins[partIndex] = i;
2052 
2053  }
2054 
2055  //############################################################//
2057  //############################################################//
2058  reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2059  this->worldSize,
2060  numPointsInParts,
2061  numGlobalPointsInParts);
2062 
2063 
2064  Teuchos::scan<int,gno_t>(
2065  *(this->comm), Teuchos::REDUCE_SUM,
2066  this->worldSize,
2067  numPointsInParts,
2068  numPointsInPartsInclusiveUptoMyIndex
2069  );
2070 
2071 
2072 
2073 
2074  /*
2075  gno_t totalSize = 0;
2076  for (int i = 0; i < this->worldSize; ++i){
2077  totalSize += numPointsInParts[i];
2078  }
2079  if (totalSize != this->numLocalCoords){
2080  std::cout << "me:" << this->myRank << " ts:" << totalSize << " nl:" << this->numLocalCoords << std::endl;
2081  }
2082  */
2083 
2084 
2085  //std::cout << "me:" << this->myRank << " ilk print" << std::endl;
2086 
2087  gno_t optimal_num = gno_t(this->numGlobalCoords/double(this->worldSize)+0.5);
2088 #ifdef printparts
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];
2094  }
2095  std::cout << "Total:" << totalSize << " ng:" << this->numGlobalCoords << std::endl;
2096  std::cout << "optimal_num:" << optimal_num << std::endl;
2097  }
2098 #endif
2099  ssize_t *extraInPart = new ssize_t [this->worldSize];
2100 
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];
2105  }
2106  if (extraExchange != 0){
2107  int addition = -1;
2108  if (extraExchange < 0) addition = 1;
2109  for (ssize_t i = 0; i < extraExchange; ++i){
2110  extraInPart[i % this->worldSize] += addition;
2111  }
2112  }
2113 
2114  //############################################################//
2116  //############################################################//
2117 
2118  int overloadedPartCount = 0;
2119  int *overloadedPartIndices = new int [this->worldSize];
2120 
2121 
2122  int underloadedPartCount = 0;
2123  int *underloadedPartIndices = new int [this->worldSize];
2124 
2125  for (int i = 0; i < this->worldSize; ++i){
2126  if(extraInPart[i] > 0){
2127  overloadedPartIndices[overloadedPartCount++] = i;
2128  }
2129  else if(extraInPart[i] < 0){
2130  underloadedPartIndices[underloadedPartCount++] = i;
2131  }
2132  }
2133 
2134  int underloadpartindex = underloadedPartCount - 1;
2135 
2136 
2137  int numPartsISendFrom = 0;
2138  int *mySendFromParts = new int[this->worldSize * 2];
2139  gno_t *mySendFromPartsCounts = new gno_t[this->worldSize * 2];
2140 
2141  int numPartsISendTo = 0;
2142  int *mySendParts = new int[this->worldSize * 2];
2143  gno_t *mySendCountsToParts = new gno_t[this->worldSize * 2];
2144 
2145 
2146  //############################################################//
2151  //############################################################//
2152  for (int i = overloadedPartCount - 1; i >= 0; --i){
2153  //get the overloaded part
2154  //the overload
2155  int overloadPart = overloadedPartIndices[i];
2156  gno_t overload = extraInPart[overloadPart];
2157  gno_t myload = numPointsInParts[overloadPart];
2158 
2159 
2160  //the inclusive load of the processors up to me
2161  gno_t inclusiveLoadUpToMe = numPointsInPartsInclusiveUptoMyIndex[overloadPart];
2162 
2163  //the exclusive load of the processors up to me
2164  gno_t exclusiveLoadUptoMe = inclusiveLoadUpToMe - myload;
2165 
2166 
2167  if (exclusiveLoadUptoMe >= overload){
2168  //this processor does not have to convert anything.
2169  //for this overloaded part.
2170  //set the extra for this processor to zero.
2171  overloadedPartIndices[i] = -1;
2172  extraInPart[overloadPart] = 0;
2173  //then consume underloaded parts.
2174  while (overload > 0){
2175  int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2176  ssize_t underload = extraInPart[nextUnderLoadedPart];
2177  ssize_t left = overload + underload;
2178 
2179  if(left >= 0){
2180  extraInPart[nextUnderLoadedPart] = 0;
2181  underloadedPartIndices[underloadpartindex--] = -1;
2182  }
2183  else {
2184  extraInPart[nextUnderLoadedPart] = left;
2185  }
2186  overload = left;
2187  }
2188  }
2189  else if (exclusiveLoadUptoMe < overload){
2190  //if the previous processors load is not enough
2191  //then this processor should convert some of its elements.
2192  gno_t mySendCount = overload - exclusiveLoadUptoMe;
2193  //how much more needed.
2194  gno_t sendAfterMe = 0;
2195  //if my load is not enough
2196  //then the next processor will continue to convert.
2197  if (mySendCount > myload){
2198  sendAfterMe = mySendCount - myload;
2199  mySendCount = myload;
2200  }
2201 
2202 
2203  //this processors will convert from overloaded part
2204  //as many as mySendCount items.
2205  mySendFromParts[numPartsISendFrom] = overloadPart;
2206  mySendFromPartsCounts[numPartsISendFrom++] = mySendCount;
2207 
2208  //first consume underloaded parts for the previous processors.
2209  while (exclusiveLoadUptoMe > 0){
2210  int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2211  ssize_t underload = extraInPart[nextUnderLoadedPart];
2212  ssize_t left = exclusiveLoadUptoMe + underload;
2213 
2214  if(left >= 0){
2215  extraInPart[nextUnderLoadedPart] = 0;
2216  underloadedPartIndices[underloadpartindex--] = -1;
2217  }
2218  else {
2219  extraInPart[nextUnderLoadedPart] = left;
2220  }
2221  exclusiveLoadUptoMe = left;
2222  }
2223 
2224  //consume underloaded parts for my load.
2225  while (mySendCount > 0){
2226  int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2227  ssize_t underload = extraInPart[nextUnderLoadedPart];
2228  ssize_t left = mySendCount + underload;
2229 
2230  if(left >= 0){
2231  mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2232  mySendCountsToParts[numPartsISendTo++] = -underload;
2233 
2234  extraInPart[nextUnderLoadedPart] = 0;
2235  underloadedPartIndices[underloadpartindex--] = -1;
2236  }
2237  else {
2238  extraInPart[nextUnderLoadedPart] = left;
2239 
2240  mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2241  mySendCountsToParts[numPartsISendTo++] = mySendCount;
2242 
2243  }
2244  mySendCount = left;
2245  }
2246  //consume underloaded parts for the load of the processors after my index.
2247  while (sendAfterMe > 0){
2248  int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2249  ssize_t underload = extraInPart[nextUnderLoadedPart];
2250  ssize_t left = sendAfterMe + underload;
2251 
2252  if(left >= 0){
2253  extraInPart[nextUnderLoadedPart] = 0;
2254  underloadedPartIndices[underloadpartindex--] = -1;
2255  }
2256  else {
2257  extraInPart[nextUnderLoadedPart] = left;
2258  }
2259  sendAfterMe = left;
2260  }
2261  }
2262  }
2263 
2264 
2265  //############################################################//
2267  //############################################################//
2268  for (int i = 0 ; i < numPartsISendFrom; ++i){
2269 
2270  //get the part from which the elements will be converted.
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];
2276 
2277  int sendCountToThisPart = mySendCountsToParts[partToSendIndex];
2278 
2279  //determine which part i should convert to
2280  //and how many to this part.
2281  if (sendCountToThisPart <= sendCount){
2282  mySendParts[partToSendIndex] = 0;
2283  mySendCountsToParts[partToSendIndex] = 0;
2284  --numPartsISendTo;
2285  sendCount -= sendCountToThisPart;
2286  }
2287  else {
2288  mySendCountsToParts[partToSendIndex] = sendCountToThisPart - sendCount;
2289  sendCountToThisPart = sendCount;
2290  sendCount = 0;
2291  }
2292 
2293 
2294  gno_t toChange = partBegins[sendFromPart];
2295  gno_t previous_begin = partBegins[partToSend];
2296 
2297  //do the conversion.
2298  for (int k = 0; k < sendCountToThisPart - 1; ++k){
2299  coordinate_grid_parts[toChange] = partToSend;
2300  toChange = partNext[toChange];
2301  }
2302  coordinate_grid_parts[toChange] = partToSend;
2303 
2304  gno_t newBegin = partNext[toChange];
2305  partNext[toChange] = previous_begin;
2306  partBegins[partToSend] = partBegins[sendFromPart];
2307  partBegins[sendFromPart] = newBegin;
2308  }
2309  }
2310 
2311  //if (this->myRank == 0) std::cout << "4" << std::endl;
2312 
2313 #ifdef printparts
2314 
2315 
2316  for (int i = 0; i < this->worldSize; ++i) numPointsInParts[i] = 0;
2317 
2318  for (int i = 0; i < this->numLocalCoords; ++i){
2319  numPointsInParts[coordinate_grid_parts[i]] += 1;
2320  }
2321 
2322  reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2323  this->worldSize,
2324  numPointsInParts,
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];
2332  }
2333  std::cout << "Total:" << totalSize << " ng:" << this->numGlobalCoords << std::endl;
2334  }
2335 #endif
2336  delete []mySendCountsToParts;
2337  delete []mySendParts;
2338  delete []mySendFromPartsCounts;
2339  delete []mySendFromParts;
2340  delete []underloadedPartIndices;
2341  delete []overloadedPartIndices;
2342  delete []extraInPart;
2343  delete []partNext;
2344  delete []partBegins;
2345  delete []numPointsInPartsInclusiveUptoMyIndex;
2346  delete []numPointsInParts;
2347  delete []numGlobalPointsInParts;
2348 
2349  delete []shiftProcCount;
2350  delete []bestDimProcs;
2351  delete []dim_slices;
2352  delete []minCoords;
2353  delete []maxCoords;
2354  }
2355 
2356  //given the part numbers for each local coordinate,
2357  //distributes the coordinates to the corresponding processors.
2358  void distribute_points(int *coordinate_grid_parts){
2359 
2360  Tpetra::Distributor distributor(comm);
2361  ArrayView<const int> pIds( coordinate_grid_parts, this->numLocalCoords);
2362  /*
2363  for (int i = 0 ; i < this->numLocalCoords; ++i){
2364  std::cout << "me:" << this->myRank << " to part:" << coordinate_grid_parts[i] << std::endl;
2365  }
2366  */
2367  gno_t numMyNewGnos = distributor.createFromSends(pIds);
2368 
2369  //std::cout << "distribution step 1 me:" << this->myRank << " numLocal:" <<numMyNewGnos << " old:" << numLocalCoords << std::endl;
2370 
2371  this->numLocalCoords = numMyNewGnos;
2372 
2373 
2374  ArrayRCP<scalar_t> recvBuf2(distributor.getTotalReceiveLength());
2375 
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];
2383  }
2384 
2385  }
2386  }
2387 
2388  //calls MJ for p = numProcs
2389  int predistributeMJ(int *coordinate_grid_parts){
2390  int coord_dim = this->coordinate_dimension;
2391 
2392  lno_t numLocalPoints = this->numLocalCoords;
2393  gno_t numGlobalPoints = this->numGlobalCoords;
2394 
2395 
2396  //T **weight = NULL;
2397  //typedef Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> tMVector_t;
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));
2400 
2401  Teuchos::Array<Teuchos::ArrayView<const scalar_t> > coordView(coord_dim);
2402 
2403 
2404 
2405  for (int i=0; i < coord_dim; i++){
2406  if(numLocalPoints > 0){
2407  Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalPoints);
2408  coordView[i] = a;
2409  } else{
2410  Teuchos::ArrayView<const scalar_t> a;
2411  coordView[i] = a;
2412  }
2413  }
2414 
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));
2417 
2418 
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;
2422 
2423  typedef Zoltan2::XpetraMultiVectorAdapter<tMVector_t> inputAdapter_t;
2424  //inputAdapter_t ia(coordsConst);
2425  inputAdapter_t ia(coordsConst,weights, stride);
2426 
2427  Teuchos::RCP <Teuchos::ParameterList> params ;
2428  params =RCP <Teuchos::ParameterList> (new Teuchos::ParameterList, true);
2429 
2430 
2431  params->set("algorithm", "multijagged");
2432  params->set("num_global_parts", this->worldSize);
2433 
2434  //TODO we need to fix the setting parts.
2435  //Although MJ sets the parts with
2436  //currently the part setting is not correct when migration is done.
2437  //params->set("migration_check_option", 2);
2438 
2439 
2441 
2442 
2443  try {
2444 #ifdef HAVE_ZOLTAN2_MPI
2445  problem = new Zoltan2::PartitioningProblem<inputAdapter_t>(&ia, params.getRawPtr(),
2446  MPI_COMM_WORLD);
2447 #else
2448  problem = new Zoltan2::PartitioningProblem<inputAdapter_t>(&ia, params.getRawPtr());
2449 #endif
2450  }
2451  CATCH_EXCEPTIONS("PartitioningProblem()")
2452 
2453  try {
2454  problem->solve();
2455  }
2456  CATCH_EXCEPTIONS("solve()")
2457 
2458  const typename inputAdapter_t::part_t *partIds = problem->getSolution().getPartListView();
2459 
2460  for (lno_t i = 0; i < this->numLocalCoords;++i){
2461  coordinate_grid_parts[i] = partIds[i];
2462  //std::cout << "me:" << this->myRank << " i:" << i << " goes to:" << partIds[i] << std::endl;
2463  }
2464  delete problem;
2465  return 0;
2466  }
2467 
2468  //calls RCP for p = numProcs
2469  int predistributeRCB(int *coordinate_grid_parts){
2470  int rank = this->myRank;
2471  int nprocs = this->worldSize;
2472  DOTS<tMVector_t> dots_;
2473 
2474  MEMORY_CHECK(rank==0 || rank==nprocs-1, "After initializing MPI");
2475 
2476 
2477  int nWeights = 0;
2478  int debugLevel=0;
2479  string memoryOn("memoryOn");
2480  string memoryOff("memoryOff");
2481  bool doMemory=false;
2482  int numGlobalParts = nprocs;
2483  int dummyTimer=0;
2484  bool remap=0;
2485 
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); // default
2492 
2493  // Process command line input
2494  CommandLineProcessor commandLine(false, true);
2495  //commandLine.setOption("size", &numGlobalCoords,
2496  // "Approximate number of global coordinates.");
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 = "";
2501 
2502  commandLine.setOption("input_file", &inputFile,
2503  "the input file for geometric generator or file input");
2504 
2505 
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");
2518 
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());
2531 
2532  CommandLineProcessor::EParseCommandLineReturn rc =
2533  commandLine.parse(0, NULL);
2534 
2535 
2536 
2537  if (rc != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
2538  if (rc == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) {
2539  if (rank==0) std::cout << "PASS" << std::endl;
2540  return 1;
2541  }
2542  else {
2543  if (rank==0) std::cout << "FAIL" << std::endl;
2544  return 0;
2545  }
2546  }
2547 
2548  //MEMORY_CHECK(doMemory && rank==0, "After processing parameters");
2549 
2550  // Create the data structure
2551 
2552  int coord_dim = this->coordinate_dimension;
2553 
2554 
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));
2557 
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);
2562  coordView[i] = a;
2563  } else{
2564  Teuchos::ArrayView<const scalar_t> a;
2565  coordView[i] = a;
2566  }
2567  }
2568 
2569  tMVector_t *tmVector = new tMVector_t( mp, coordView.view(0, coord_dim), coord_dim);
2570 
2571  dots_.coordinates = tmVector;
2572  dots_.weights.resize(nWeights);
2573 
2574 
2575  MEMORY_CHECK(doMemory && rank==0, "After creating input");
2576 
2577  // Now call Zoltan to partition the problem.
2578 
2579  float ver;
2580  int aok = Zoltan_Initialize(0,NULL, &ver);
2581 
2582  if (aok != 0){
2583  printf("Zoltan_Initialize failed\n");
2584  exit(0);
2585  }
2586 
2587  struct Zoltan_Struct *zz;
2588  zz = Zoltan_Create(MPI_COMM_WORLD);
2589 
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());
2599  oss.str("");
2600  oss << debugLevel;
2601  Zoltan_Set_Param(zz, "DEBUG_LEVEL", oss.str().c_str());
2602 
2603  if (remap)
2604  Zoltan_Set_Param(zz, "REMAP", "1");
2605  else
2606  Zoltan_Set_Param(zz, "REMAP", "0");
2607 
2608  if (objective != balanceCount){
2609  oss.str("");
2610  oss << nWeights;
2611  Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", oss.str().c_str());
2612 
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");
2619  }
2620  else{
2621  Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", "0");
2622  }
2623 
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_);
2628 
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;
2633 
2634  MEMORY_CHECK(doMemory && rank==0, "Before Zoltan_LB_Partition");
2635 
2636 
2637  aok = Zoltan_LB_Partition(zz, &changes, &numGidEntries, &numLidEntries,
2638  &numImport, &importGlobalGids, &importLocalGids,
2639  &importProcs, &importToPart,
2640  &numExport, &exportGlobalGids, &exportLocalGids,
2641  &exportProcs, &exportToPart);
2642 
2643 
2644  MEMORY_CHECK(doMemory && rank==0, "After Zoltan_LB_Partition");
2645 
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");
2650 
2651  delete dots_.coordinates;
2652  return 0;
2653 }
2655  int *coordinate_grid_parts = new int[this->numLocalCoords];
2656  switch (this->predistribution){
2657  case 1:
2658  this->predistributeRCB(coordinate_grid_parts);
2659  break;
2660  case 2:
2661 
2662  this->predistributeMJ(coordinate_grid_parts);
2663  break;
2664  case 3:
2665  //block
2666  blockPartition(coordinate_grid_parts);
2667  break;
2668  }
2669  this->distribute_points(coordinate_grid_parts);
2670 
2671  delete []coordinate_grid_parts;
2672 
2673 
2674  }
2675 
2676  //############################################################//
2678  //############################################################//
2679 
2680 
2682  return this->numWeightsPerCoord;
2683  }
2685  return this->coordinate_dimension;
2686  }
2688  return this->numLocalCoords;
2689  }
2691  return this->numGlobalCoords;
2692  }
2693 
2695  return this->coords;
2696  }
2697 
2698  scalar_t **getLocalWeightsView(){
2699  return this->wghts;
2700  }
2701 
2702  void getLocalCoordinatesCopy( scalar_t ** c){
2703  for(int ii = 0; ii < this->coordinate_dimension; ++ii){
2704 #ifdef HAVE_ZOLTAN2_OMP
2705 #pragma omp parallel for
2706 #endif
2707  for (lno_t i = 0; i < this->numLocalCoords; ++i){
2708  c[ii][i] = this->coords[ii][i];
2709  }
2710  }
2711  }
2712 
2713  void getLocalWeightsCopy(scalar_t **w){
2714  for(int ii = 0; ii < this->numWeightsPerCoord; ++ii){
2715 #ifdef HAVE_ZOLTAN2_OMP
2716 #pragma omp parallel for
2717 #endif
2718  for (lno_t i = 0; i < this->numLocalCoords; ++i){
2719  w[ii][i] = this->wghts[ii][i];
2720  }
2721  }
2722  }
2723 };
2724 }
2725 
2726 #endif /* GEOMETRICGENERATOR */
RectangleHole(CoordinatePoint< T > center_, T edge_x_, T edge_y_)
virtual weighttype get1DWeight(T x)
virtual bool isInArea(CoordinatePoint< T > dot)
#define INVALID(STR)
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
paramName
Definition: xml2dox.py:207
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)
Tpetra::Map< zlno_t, zgno_t, znode_t > tMap_t
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_)
CoordinateUniformDistribution(gno_t np_, int dim, T l_x, T r_x, T l_y, T r_y, T l_z, T r_z, int wSize)
SphereHole(CoordinatePoint< T > center_, T edge_)
CircleHole(CoordinatePoint< T > center_, T edge_)
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[]
virtual weighttype getWeight(CoordinatePoint< T > p)
int getDim(void *data, int *ierr)
void getObjList(void *data, int numGid, int numLid, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int num_wgts, float *obj_wgts, int *ierr)
outfile
Definition: xml2dox.py:11
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.
void distribute_points(int *coordinate_grid_parts)
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
PartitioningProblem sets up partitioning problems for the user.
#define MAX_ITER_ALLOWED
GeometricGenerator(Teuchos::ParameterList &params, const RCP< const Teuchos::Comm< int > > &comm_)
Tpetra::MultiVector< double, int, int > tMVector_t
int getNumWeights()
##END Predistribution functions######################//
CoordinatePoint< T > center
#define INVALID_SHAPE_ARG(SHAPE, REQUIRED)
#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)
#define epsilon
SquareHole(CoordinatePoint< T > center_, T edge_)
void solve(bool updateInputData=true)
Direct the problem to create a solution.
virtual CoordinatePoint< T > getPoint(gno_t pindex, unsigned int &state)
CubeHole(CoordinatePoint< T > center_, T edge_)
#define MAX_WEIGHT_DIM
virtual bool isInArea(CoordinatePoint< T > dot)
void getMinMaxCoords(scalar_t *globalMaxCoords, scalar_t *globalMinCoords)
#define MEMORY_CHECK(iPrint, msg)