Zoltan2
Loading...
Searching...
No Matches
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
76using Teuchos::CommandLineProcessor;
77using Teuchos::ArrayView;
78using std::string;
79using Teuchos::ArrayRCP;
80
81namespace 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
107template <typename tMVector_t>
108class DOTS{
109public:
110 std::vector<std::vector<float> > weights;
112};
113
114template <typename tMVector_t>
115int getNumObj(void *data, int *ierr)
116{
117 *ierr = 0;
118 DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
119 return dots_->coordinates->getLocalLength();
120}
122template <typename tMVector_t>
123void 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
156template <typename tMVector_t>
157int 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
167template <typename tMVector_t>
168void 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()->getLocalElementList().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
197const std::string shapes[] = {"SQUARE", "RECTANGLE", "CIRCLE", "CUBE", "RECTANGULAR_PRISM", "SPHERE"};
198#define SHAPE_COUNT 6
199
201const 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
212const std::string weight_distribution_string = "WeightDistribution-";
213
214template <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
247template <typename T>
248class Hole{
249public:
250 CoordinatePoint<T> center;
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
264template <typename T>
265class SquareHole:public Hole<T>{
266public:
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
274template <typename T>
275class RectangleHole:public Hole<T>{
276public:
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
284template <typename T>
285class CircleHole:public Hole<T>{
286public:
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
294template <typename T>
295class CubeHole:public Hole<T>{
296public:
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
304template <typename T>
306public:
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
314template <typename T>
315class SphereHole:public Hole<T>{
316public:
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
328template <typename T, typename weighttype>
330public:
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;
337};
338
357template <typename T, typename weighttype>
358class 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;
367public:
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
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
470template <typename T, typename lno_t, typename gno_t>
472public:
479
480 CoordinateDistribution(gno_t np_, int dim, int wSize) :
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
644template <typename T, typename lno_t, typename gno_t>
646public:
647 CoordinatePoint<T> center;
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.
673 CoordinatePoint <T> p;
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
694private:
695 T normalDist(T center_, T sd, unsigned int &state) {
696 T polarsqrt, normalsquared, normal1, normal2;
697 do {
698 normal1=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
699 normal2=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
700 normalsquared=normal1*normal1+normal2*normal2;
701 } while ( normalsquared>=1.0 || normalsquared == 0.0);
702
703 polarsqrt=sqrt(-2.0*log(normalsquared)/normalsquared);
704 return normal2*polarsqrt*sd + center_;
705 }
706};
707
708template <typename T, typename lno_t, typename gno_t>
710public:
717
718
719 virtual T getXCenter(){
720 return (rightMostx - leftMostx)/2 + leftMostx;
721 }
722 virtual T getXRadius(){
723 return (rightMostx - leftMostx)/2;
724 }
725
726
727 CoordinateUniformDistribution(gno_t np_, int dim, T l_x, T r_x, T l_y, T r_y,
728 T l_z, T r_z, int wSize ) :
729 CoordinateDistribution<T,lno_t,gno_t>(np_,dim,wSize),
730 leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y),
731 leftMostz(l_z), rightMostz(r_z){}
732
734 virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
735
736
737 //pindex = 0; //not used in uniform dist.
738 CoordinatePoint <T> p;
739 for(int i = 0; i < this->dimension; ++i){
740 switch(i){
741 case 0:
742 p.x = uniformDist(this->leftMostx, this->rightMostx, state);
743 break;
744 case 1:
745 p.y = uniformDist(this->leftMosty, this->rightMosty, state);
746 break;
747 case 2:
748 p.z = uniformDist(this->leftMostz, this->rightMostz, state);
749 break;
750 default:
751 throw "unsupported dimension";
752 }
753 }
754 return p;
755 }
756
757private:
758
759 T uniformDist(T a, T b, unsigned int &state)
760 {
761 return (b-a)*(rand_r(&state) / double(RAND_MAX)) + a;
762 }
763};
764
765template <typename T, typename lno_t, typename gno_t>
767public:
775 //T currentX, currentY, currentZ;
780
781 virtual T getXCenter(){
782 return (rightMostx - leftMostx)/2 + leftMostx;
783 }
784 virtual T getXRadius(){
785 return (rightMostx - leftMostx)/2;
786 }
787
788
789 CoordinateGridDistribution(gno_t alongX, gno_t alongY, gno_t alongZ, int dim,
790 T l_x, T r_x, T l_y, T r_y, T l_z, T r_z ,
791 int myRank_, int wSize) :
792 CoordinateDistribution<T,lno_t,gno_t>(alongX * alongY * alongZ,dim,wSize),
793 leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y), leftMostz(l_z), rightMostz(r_z), myRank(myRank_){
794 //currentX = leftMostx, currentY = leftMosty, currentZ = leftMostz;
795 this->processCnt = 0;
796 this->along_X = alongX; this->along_Y = alongY; this->along_Z = alongZ;
797
798 if(this->along_X > 1)
799 this->xstep = (rightMostx - leftMostx) / (alongX - 1);
800 else
801 this->xstep = 1;
802 if(this->along_Y > 1)
803 this->ystep = (rightMosty - leftMosty) / (alongY - 1);
804 else
805 this->ystep = 1;
806 if(this->along_Z > 1)
807 this->zstep = (rightMostz - leftMostz) / (alongZ - 1);
808 else
809 this->zstep = 1;
810 xshift = 0; yshift = 0; zshift = 0;
811 }
812
814 virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
815 //lno_t before = processCnt + this->assignedPrevious;
816 //std::cout << "before:" << processCnt << " " << this->assignedPrevious << std::endl;
817 //lno_t xshift = 0, yshift = 0, zshift = 0;
818
819 //lno_t tmp = before % (this->along_X * this->along_Y);
820 //xshift = tmp % this->along_X;
821 //yshift = tmp / this->along_X;
822 //zshift = before / (this->along_X * this->along_Y);
823
824 state = 0; //not used here
825 this->zshift = pindex / (along_X * along_Y);
826 this->yshift = (pindex % (along_X * along_Y)) / along_X;
827 this->xshift = (pindex % (along_X * along_Y)) % along_X;
828
829 //this->xshift = pindex / (along_Z * along_Y);
830 // this->zshift = (pindex % (along_Z * along_Y)) / along_Y;
831 //this->yshift = (pindex % (along_Z * along_Y)) % along_Y;
832
833 CoordinatePoint <T> p;
834 p.x = xshift * this->xstep + leftMostx;
835 p.y = yshift * this->ystep + leftMosty;
836 p.z = zshift * this->zstep + leftMostz;
837/*
838 ++xshift;
839 if(xshift == this->along_X){
840 ++yshift;
841 xshift = 0;
842 if(yshift == this->along_Y){
843 ++zshift;
844 yshift = 0;
845 }
846 }
847*/
848 /*
849 if(this->processCnt == 0){
850 this->xshift = this->assignedPrevious / (along_Z * along_Y);
851 //this->yshift = (this->assignedPrevious % (along_X * along_Y)) / along_X;
852 this->zshift = (this->assignedPrevious % (along_Z * along_Y)) / along_Y;
853 //this->xshift = (this->assignedPrevious % (along_X * along_Y)) % along_X;
854 this->yshift = (this->assignedPrevious % (along_Z * along_Y)) % along_Y;
855 ++this->processCnt;
856 }
857
858 CoordinatePoint <T> p;
859 p.x = xshift * this->xstep + leftMostx;
860 p.y = yshift * this->ystep + leftMosty;
861 p.z = zshift * this->zstep + leftMostz;
862
863 ++yshift;
864 if(yshift == this->along_Y){
865 ++zshift;
866 yshift = 0;
867 if(zshift == this->along_Z){
868 ++xshift;
869 zshift = 0;
870 }
871
872 }
873 */
874 /*
875 if(this->requested - 1 > this->processCnt){
876 this->processCnt++;
877 } else {
878 std::cout << "req:" << this->requested << " pr:" <<this->processCnt << std::endl;
879 std::cout << "p:" << p.x << " " << p.y << " " << p.z << std::endl ;
880 std::cout << "s:" << xshift << " " << yshift << " " << zshift << std::endl ;
881 std::cout << "st:" << this->xstep << " " << this->ystep << " " << this->zstep << std::endl ;
882 }
883 */
884 return p;
885 }
886
887private:
888
889};
890
891template <typename scalar_t, typename lno_t, typename gno_t, typename node_t>
893private:
894 Hole<scalar_t> **holes; //to represent if there is any hole in the input
895 int holeCount;
896 int coordinate_dimension; //dimension of the geometry
897 gno_t numGlobalCoords; //global number of coordinates requested to be created.
898 lno_t numLocalCoords;
899 float *loadDistributions; //sized as the number of processors, the load of each processor.
900 bool loadDistSet;
901 bool distinctCoordSet;
902 CoordinateDistribution<scalar_t, lno_t,gno_t> **coordinateDistributions;
903 int distributionCount;
904 //CoordinatePoint<scalar_t> *points;
905 scalar_t **coords;
906 scalar_t **wghts;
907 WeightDistribution<scalar_t,scalar_t> **wd;
908 int numWeightsPerCoord;
909 int predistribution;
910 RCP<const Teuchos::Comm<int> > comm;
911 //RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmVector;
912 //RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmwVector;
913 int worldSize;
914 int myRank;
915 scalar_t minx;
916 scalar_t maxx;
917 scalar_t miny;
918 scalar_t maxy;
919 scalar_t minz;
920 scalar_t maxz;
921 std::string outfile;
922 float perturbation_ratio;
923
924 typedef Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> tMVector_t;
925 typedef Tpetra::Map<lno_t, gno_t, node_t> tMap_t;
926
927
928 template <typename tt>
929 tt getParamVal( const Teuchos::ParameterEntry& pe, const std::string &paramname){
930 tt returnVal;
931 try {
932 returnVal = Teuchos::getValue<tt>(pe);
933 }
934 catch (...){
935 throw INVALID(paramname);
936 }
937 return returnVal;
938 }
939
940 int countChar (std::string inStr, char cntChar){
941 int cnt = 0;
942 for (unsigned int i = 0; i < inStr.size(); ++i){
943 if (inStr[i] == cntChar) {
944 cnt++;
945 }
946 }
947 return cnt;
948 }
949
950 template <typename tt>
951 tt fromString(std::string obj){
952 std::stringstream ss (std::stringstream::in | std::stringstream::out);
953 ss << obj;
954 tt tmp;
955 ss >> tmp;
956 if (ss.fail()){
957 throw "Cannot convert string " + obj;
958 }
959 return tmp;
960 }
961
962
963 void splitString(std::string inStr, char splitChar, std::string *outSplittedStr){
964 std::stringstream ss (std::stringstream::in | std::stringstream::out);
965 ss << inStr;
966 int cnt = 0;
967 while (!ss.eof()){
968 std::string tmp = "";
969 std::getline(ss, tmp ,splitChar);
970 outSplittedStr[cnt++] = tmp;
971 }
972 }
973
974 /*
975 void getDistinctCoordinateDescription(std::string distinctDescription){
976
977 this->distinctCoordinates = new bool[this->dimension];
978 if(distinctDescription != ""){
979 int argCnt = this->countChar(distinctDescription, ',') + 1;
980 if(argCnt != this->dimension) {
981 throw "Invalid parameter count for distinct_coordinates. Size should be equal to dimension.";
982 }
983
984 std::string *splittedStr = new std::string[argCnt];
985 splitString(holeDescription, ',', splittedStr);
986
987 for(int i = 0; i < argCnt; ++i){
988 if(splittedStr[i] == "T"){
989 distinctCoordinates[i] = true;
990 } else if(splittedStr[i] == "F"){
991 distinctCoordinates[i] = false;
992 } else {
993 throw "Invalid parameter " + splittedStr[i] + " for distinct_coordinates.";
994 }
995 }
996 delete []splittedStr;
997 }
998 else {
999 for(int i = 0; i < this->dimension; ++i){
1000 distinctCoordinates[i] = false;
1001 }
1002 }
1003 }
1004 */
1005
1006
1007 void getCoordinateDistributions(std::string coordinate_distributions){
1008 if(coordinate_distributions == ""){
1009 throw "At least one distribution function must be provided to coordinate generator.";
1010 }
1011
1012 int argCnt = this->countChar(coordinate_distributions, ',') + 1;
1013 std::string *splittedStr = new std::string[argCnt];
1014 splitString(coordinate_distributions, ',', splittedStr);
1015 coordinateDistributions = (CoordinateDistribution<scalar_t, lno_t,gno_t> **) malloc(sizeof (CoordinateDistribution<scalar_t, lno_t,gno_t> *) * 1);
1016 for(int i = 0; i < argCnt; ){
1017 coordinateDistributions = (CoordinateDistribution<scalar_t, lno_t,gno_t> **)realloc((void *)coordinateDistributions, (this->distributionCount + 1)* sizeof(CoordinateDistribution<scalar_t, lno_t,gno_t> *));
1018
1019 std::string distName = splittedStr[i++];
1020 gno_t np_ = 0;
1021 if(distName == "NORMAL"){
1022 int reqArg = 5;
1023 if (this->coordinate_dimension == 3){
1024 reqArg = 7;
1025 }
1026 if(i + reqArg > argCnt) {
1027 std::string tmp = Teuchos::toString<int>(reqArg);
1028 throw INVALID_SHAPE_ARG(distName, tmp);
1029 }
1030 np_ = fromString<gno_t>(splittedStr[i++]);
1031 CoordinatePoint<scalar_t> pp;
1032
1033 pp.x = fromString<scalar_t>(splittedStr[i++]);
1034 pp.y = fromString<scalar_t>(splittedStr[i++]);
1035 pp.z = 0;
1036 if(this->coordinate_dimension == 3){
1037 pp.z = fromString<scalar_t>(splittedStr[i++]);
1038 }
1039
1040 scalar_t sd_x = fromString<scalar_t>(splittedStr[i++]);
1041 scalar_t sd_y = fromString<scalar_t>(splittedStr[i++]);
1042 scalar_t sd_z = 0;
1043 if(this->coordinate_dimension == 3){
1044 sd_z = fromString<scalar_t>(splittedStr[i++]);
1045 }
1046 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 );
1047
1048 } else if(distName == "UNIFORM" ){
1049 int reqArg = 5;
1050 if (this->coordinate_dimension == 3){
1051 reqArg = 7;
1052 }
1053 if(i + reqArg > argCnt) {
1054 std::string tmp = Teuchos::toString<int>(reqArg);
1055 throw INVALID_SHAPE_ARG(distName, tmp);
1056 }
1057 np_ = fromString<gno_t>(splittedStr[i++]);
1058 scalar_t l_x = fromString<scalar_t>(splittedStr[i++]);
1059 scalar_t r_x = fromString<scalar_t>(splittedStr[i++]);
1060 scalar_t l_y = fromString<scalar_t>(splittedStr[i++]);
1061 scalar_t r_y = fromString<scalar_t>(splittedStr[i++]);
1062
1063 scalar_t l_z = 0, r_z = 0;
1064
1065 if(this->coordinate_dimension == 3){
1066 l_z = fromString<scalar_t>(splittedStr[i++]);
1067 r_z = fromString<scalar_t>(splittedStr[i++]);
1068 }
1069
1070 this->coordinateDistributions[this->distributionCount++] = new CoordinateUniformDistribution<scalar_t, lno_t,gno_t>( np_, this->coordinate_dimension, l_x, r_x, l_y, r_y, l_z, r_z, this->worldSize );
1071 } else if (distName == "GRID"){
1072 int reqArg = 6;
1073 if(this->coordinate_dimension == 3){
1074 reqArg = 9;
1075 }
1076 if(i + reqArg > argCnt) {
1077 std::string tmp = Teuchos::toString<int>(reqArg);
1078 throw INVALID_SHAPE_ARG(distName, tmp);
1079 }
1080
1081 gno_t np_x = fromString<gno_t>(splittedStr[i++]);
1082 gno_t np_y = fromString<gno_t>(splittedStr[i++]);
1083 gno_t np_z = 1;
1084
1085
1086 if(this->coordinate_dimension == 3){
1087 np_z = fromString<gno_t>(splittedStr[i++]);
1088 }
1089
1090 np_ = np_x * np_y * np_z;
1091 scalar_t l_x = fromString<scalar_t>(splittedStr[i++]);
1092 scalar_t r_x = fromString<scalar_t>(splittedStr[i++]);
1093 scalar_t l_y = fromString<scalar_t>(splittedStr[i++]);
1094 scalar_t r_y = fromString<scalar_t>(splittedStr[i++]);
1095
1096 scalar_t l_z = 0, r_z = 0;
1097
1098 if(this->coordinate_dimension == 3){
1099 l_z = fromString<scalar_t>(splittedStr[i++]);
1100 r_z = fromString<scalar_t>(splittedStr[i++]);
1101 }
1102
1103 if(np_x < 1 || np_z < 1 || np_y < 1 ){
1104 throw "Provide at least 1 point along each dimension for grid test.";
1105 }
1106 //std::cout << "ly:" << l_y << " ry:" << r_y << std::endl;
1107 this->coordinateDistributions[this->distributionCount++] = new CoordinateGridDistribution<scalar_t, lno_t,gno_t>
1108 (np_x, np_y,np_z, this->coordinate_dimension, l_x, r_x,l_y, r_y, l_z, r_z , this->myRank, this->worldSize);
1109
1110 }
1111 else {
1112 std::string tmp = Teuchos::toString<int>(this->coordinate_dimension);
1113 throw INVALIDSHAPE(distName, tmp);
1114 }
1115 this->numGlobalCoords += (gno_t) np_;
1116 }
1117 delete []splittedStr;
1118 }
1119
1120 void getProcLoadDistributions(std::string proc_load_distributions){
1121
1122 this->loadDistributions = new float[this->worldSize];
1123 if(proc_load_distributions == ""){
1124 float r = 1.0 / this->worldSize;
1125 for(int i = 0; i < this->worldSize; ++i){
1126 this->loadDistributions[i] = r;
1127 }
1128 } else{
1129
1130
1131 int argCnt = this->countChar(proc_load_distributions, ',') + 1;
1132 if(argCnt != this->worldSize) {
1133 throw "Invalid parameter count load distributions. Given " + Teuchos::toString<int>(argCnt) + " processor size is " + Teuchos::toString<int>(this->worldSize);
1134 }
1135 std::string *splittedStr = new std::string[argCnt];
1136 splitString(proc_load_distributions, ',', splittedStr);
1137 for(int i = 0; i < argCnt; ++i){
1138 this->loadDistributions[i] = fromString<float>(splittedStr[i]);
1139 }
1140 delete []splittedStr;
1141
1142
1143 float sum = 0;
1144 for(int i = 0; i < this->worldSize; ++i){
1145 sum += this->loadDistributions[i];
1146 }
1147 if (fabs(sum - 1.0) > 10*std::numeric_limits<float>::epsilon()){
1148 throw "Processor load ratios do not sum to 1.0.";
1149 }
1150 }
1151
1152 }
1153
1154 void getHoles(std::string holeDescription){
1155 if(holeDescription == ""){
1156 return;
1157 }
1158 this->holes = (Hole<scalar_t> **) malloc(sizeof (Hole <scalar_t>*));
1159 int argCnt = this->countChar(holeDescription, ',') + 1;
1160 std::string *splittedStr = new std::string[argCnt];
1161 splitString(holeDescription, ',', splittedStr);
1162
1163 for(int i = 0; i < argCnt; ){
1164 holes = (Hole<scalar_t> **)realloc((void *)holes, (this->holeCount + 1)* sizeof(Hole<scalar_t> *));
1165
1166 std::string shapeName = splittedStr[i++];
1167 if(shapeName == "SQUARE" && this->coordinate_dimension == 2){
1168 if(i + 3 > argCnt) {
1169 throw INVALID_SHAPE_ARG(shapeName, "3");
1170 }
1171 CoordinatePoint<scalar_t> pp;
1172 pp.x = fromString<scalar_t>(splittedStr[i++]);
1173 pp.y = fromString<scalar_t>(splittedStr[i++]);
1174 scalar_t edge = fromString<scalar_t>(splittedStr[i++]);
1175 this->holes[this->holeCount++] = new SquareHole<scalar_t>(pp, edge);
1176 } else if(shapeName == "RECTANGLE" && this->coordinate_dimension == 2){
1177 if(i + 4 > argCnt) {
1178 throw INVALID_SHAPE_ARG(shapeName, "4");
1179 }
1180 CoordinatePoint<scalar_t> pp;
1181 pp.x = fromString<scalar_t>(splittedStr[i++]);
1182 pp.y = fromString<scalar_t>(splittedStr[i++]);
1183 scalar_t edgex = fromString<scalar_t>(splittedStr[i++]);
1184 scalar_t edgey = fromString<scalar_t>(splittedStr[i++]);
1185
1186 this->holes[this->holeCount++] = new RectangleHole<scalar_t>(pp, edgex,edgey);
1187 } else if(shapeName == "CIRCLE" && this->coordinate_dimension == 2){
1188 if(i + 3 > argCnt) {
1189 throw INVALID_SHAPE_ARG(shapeName, "3");
1190 }
1191 CoordinatePoint<scalar_t> pp;
1192 pp.x = fromString<scalar_t>(splittedStr[i++]);
1193 pp.y = fromString<scalar_t>(splittedStr[i++]);
1194 scalar_t r = fromString<scalar_t>(splittedStr[i++]);
1195 this->holes[this->holeCount++] = new CircleHole<scalar_t>(pp, r);
1196 } else if(shapeName == "CUBE" && this->coordinate_dimension == 3){
1197 if(i + 4 > argCnt) {
1198 throw INVALID_SHAPE_ARG(shapeName, "4");
1199 }
1200 CoordinatePoint<scalar_t> pp;
1201 pp.x = fromString<scalar_t>(splittedStr[i++]);
1202 pp.y = fromString<scalar_t>(splittedStr[i++]);
1203 pp.z = fromString<scalar_t>(splittedStr[i++]);
1204 scalar_t edge = fromString<scalar_t>(splittedStr[i++]);
1205 this->holes[this->holeCount++] = new CubeHole<scalar_t>(pp, edge);
1206 } else if(shapeName == "RECTANGULAR_PRISM" && this->coordinate_dimension == 3){
1207 if(i + 6 > argCnt) {
1208 throw INVALID_SHAPE_ARG(shapeName, "6");
1209 }
1210 CoordinatePoint<scalar_t> pp;
1211 pp.x = fromString<scalar_t>(splittedStr[i++]);
1212 pp.y = fromString<scalar_t>(splittedStr[i++]);
1213 pp.z = fromString<scalar_t>(splittedStr[i++]);
1214 scalar_t edgex = fromString<scalar_t>(splittedStr[i++]);
1215 scalar_t edgey = fromString<scalar_t>(splittedStr[i++]);
1216 scalar_t edgez = fromString<scalar_t>(splittedStr[i++]);
1217 this->holes[this->holeCount++] = new RectangularPrismHole<scalar_t>(pp, edgex, edgey, edgez);
1218
1219 } else if(shapeName == "SPHERE" && this->coordinate_dimension == 3){
1220 if(i + 4 > argCnt) {
1221 throw INVALID_SHAPE_ARG(shapeName, "4");
1222 }
1223 CoordinatePoint<scalar_t> pp;
1224 pp.x = fromString<scalar_t>(splittedStr[i++]);
1225 pp.y = fromString<scalar_t>(splittedStr[i++]);
1226 pp.z = fromString<scalar_t>(splittedStr[i++]);
1227 scalar_t r = fromString<scalar_t>(splittedStr[i++]);
1228 this->holes[this->holeCount++] = new SphereHole<scalar_t>(pp, r);
1229 } else {
1230 std::string tmp = Teuchos::toString<int>(this->coordinate_dimension);
1231 throw INVALIDSHAPE(shapeName, tmp);
1232 }
1233 }
1234 delete [] splittedStr;
1235 }
1236
1237 void getWeightDistribution(std::string *weight_distribution_arr, int wdimension){
1238 int wcount = 0;
1239
1240 this->wd = new WeightDistribution<scalar_t,scalar_t> *[wdimension];
1241 for(int ii = 0; ii < MAX_WEIGHT_DIM; ++ii){
1242 std::string weight_distribution = weight_distribution_arr[ii];
1243 if(weight_distribution == "") continue;
1244 if(wcount == wdimension) {
1245 throw "Weight Dimension is provided as " + Teuchos::toString<int>(wdimension) + ". More weight distribution is provided.";
1246 }
1247
1248 int count = this->countChar(weight_distribution, ' ');
1249 std::string *splittedStr = new string[count + 1];
1250 this->splitString(weight_distribution, ' ', splittedStr);
1251 //std::cout << count << std::endl;
1252 scalar_t c=1;
1253 scalar_t a1=0,a2=0,a3=0;
1254 scalar_t x1=0,y1=0,z1=0;
1255 scalar_t b1=1,b2=1,b3=1;
1256 scalar_t *steps = NULL;
1257 scalar_t *values= NULL;
1258 int stepCount = 0;
1259 int valueCount = 1;
1260
1261 for (int i = 1; i < count + 1; ++i){
1262 int assignmentCount = this->countChar(splittedStr[i], '=');
1263 if (assignmentCount != 1){
1264 throw "Error at the argument " + splittedStr[i];
1265 }
1266 std::string parameter_vs_val[2];
1267 this->splitString(splittedStr[i], '=', parameter_vs_val);
1268 std::string parameter = parameter_vs_val[0];
1269 std::string value = parameter_vs_val[1];
1270 //std::cout << "parameter:" << parameter << " value:" << value << std::endl;
1271
1272 if (parameter == "a1"){
1273 a1 = this->fromString<scalar_t>(value);
1274 }
1275 else if (parameter == "a2"){
1276 if(this->coordinate_dimension > 1){
1277 a2 = this->fromString<scalar_t>(value);
1278 }
1279 else {
1280 throw parameter+ " argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1281 }
1282
1283 }
1284 else if (parameter == "a3"){
1285 if(this->coordinate_dimension > 2){
1286 a3 = this->fromString<scalar_t>(value);
1287 }
1288 else {
1289 throw parameter+ " argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1290 }
1291 }
1292 else if (parameter == "b1"){
1293 b1 = this->fromString<scalar_t>(value);
1294 }
1295 else if (parameter == "b2"){
1296 if(this->coordinate_dimension > 1){
1297 b2 = this->fromString<scalar_t>(value);
1298 }
1299 else {
1300 throw parameter+ " argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1301 }
1302 }
1303 else if (parameter == "b3"){
1304
1305 if(this->coordinate_dimension > 2){
1306 b3 = this->fromString<scalar_t>(value);
1307 }
1308 else {
1309 throw parameter+ " argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1310 }
1311 }
1312 else if (parameter == "c"){
1313 c = this->fromString<scalar_t>(value);
1314 }
1315 else if (parameter == "x1"){
1316 x1 = this->fromString<scalar_t>(value);
1317 }
1318 else if (parameter == "y1"){
1319 if(this->coordinate_dimension > 1){
1320 y1 = this->fromString<scalar_t>(value);
1321 }
1322 else {
1323 throw parameter+ " argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1324 }
1325 }
1326 else if (parameter == "z1"){
1327 if(this->coordinate_dimension > 2){
1328 z1 = this->fromString<scalar_t>(value);
1329 }
1330 else {
1331 throw parameter+ " argument is not valid when dimension is " + Teuchos::toString<int>(this->coordinate_dimension);
1332 }
1333 }
1334 else if (parameter == "steps"){
1335 stepCount = this->countChar(value, ',') + 1;
1336 std::string *stepstr = new std::string[stepCount];
1337 this->splitString(value, ',', stepstr);
1338 steps = new scalar_t[stepCount];
1339 for (int j = 0; j < stepCount; ++j){
1340 steps[j] = this->fromString<scalar_t>(stepstr[j]);
1341 }
1342 delete [] stepstr;
1343 }
1344 else if (parameter == "values"){
1345 valueCount = this->countChar(value, ',') + 1;
1346 std::string *stepstr = new std::string[valueCount];
1347 this->splitString(value, ',', stepstr);
1348 values = new scalar_t[valueCount];
1349 for (int j = 0; j < valueCount; ++j){
1350 values[j] = this->fromString<scalar_t>(stepstr[j]);
1351 }
1352 delete [] stepstr;
1353 }
1354 else {
1355 throw "Invalid parameter name at " + splittedStr[i];
1356 }
1357 }
1358
1359 delete []splittedStr;
1360 if(stepCount + 1!= valueCount){
1361 throw "Step count: " + Teuchos::toString<int>(stepCount) + " must be 1 less than value count: " + Teuchos::toString<int>(valueCount);
1362 }
1363
1364
1365 this->wd[wcount++] = new SteppedEquation<scalar_t,scalar_t>(a1, a2, a3, b1, b2, b3, c, x1, y1, z1, steps, values, stepCount);
1366
1367 if(stepCount > 0){
1368 delete [] steps;
1369 delete [] values;
1370
1371 }
1372 }
1373 if(wcount != this->numWeightsPerCoord){
1374 throw "Weight Dimension is provided as " + Teuchos::toString<int>(wdimension) + ". But " + Teuchos::toString<int>(wcount)+" weight distributions are provided.";
1375 }
1376 }
1377
1378 void parseParams(const Teuchos::ParameterList & params){
1379 try {
1380 std::string holeDescription = "";
1381 std::string proc_load_distributions = "";
1382 std::string distinctDescription = "";
1383 std::string coordinate_distributions = "";
1384 std::string numWeightsPerCoord_parameters[MAX_WEIGHT_DIM];
1385 for (int i = 0; i < MAX_WEIGHT_DIM; ++i){
1386 numWeightsPerCoord_parameters[i] = "";
1387 }
1388
1389
1390 for (Teuchos::ParameterList::ConstIterator pit = params.begin(); pit != params.end(); ++pit ){
1391 const std::string &paramName = params.name(pit);
1392
1393 const Teuchos::ParameterEntry &pe = params.entry(pit);
1394
1395 if(paramName.find("hole-") == 0){
1396 if(holeDescription == ""){
1397 holeDescription = getParamVal<std::string>(pe, paramName);
1398 }
1399 else {
1400 holeDescription +=","+getParamVal<std::string>(pe, paramName);
1401 }
1402 }
1403 else if(paramName.find("distribution-") == 0){
1404 if(coordinate_distributions == "")
1405 coordinate_distributions = getParamVal<std::string>(pe, paramName);
1406 else
1407 coordinate_distributions += ","+getParamVal<std::string>(pe, paramName);
1408 //std::cout << coordinate_distributions << std::endl;
1409 //TODO coordinate distribution description
1410 }
1411
1412 else if (paramName.find(weight_distribution_string) == 0){
1413 std::string weight_dist_param = paramName.substr(weight_distribution_string.size());
1414 int dash_pos = weight_dist_param.find("-");
1415 std::string distribution_index_string = weight_dist_param.substr(0, dash_pos);
1416 int distribution_index = fromString<int>(distribution_index_string);
1417
1418 if(distribution_index >= MAX_WEIGHT_DIM){
1419 throw "Given distribution index:" + distribution_index_string + " larger than maximum allowed number of weights:" + Teuchos::toString<int>(MAX_WEIGHT_DIM);
1420 }
1421 numWeightsPerCoord_parameters[distribution_index] += " " + weight_dist_param.substr(dash_pos + 1)+ "="+ getParamVal<std::string>(pe, paramName);
1422 }
1423 else if(paramName == "dim"){
1424 int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
1425 if(dim < 2 || dim > 3){
1426 throw INVALID(paramName);
1427 } else {
1428 this->coordinate_dimension = dim;
1429 }
1430 }
1431 else if(paramName == "wdim"){
1432 int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
1433 if(dim < 1 && dim > MAX_WEIGHT_DIM){
1434 throw INVALID(paramName);
1435 } else {
1436 this->numWeightsPerCoord = dim;
1437 }
1438 }
1439 else if(paramName == "predistribution"){
1440 int pre_distribution = fromString<int>(getParamVal<std::string>(pe, paramName));
1441 if(pre_distribution < 0 || pre_distribution > 3){
1442 throw INVALID(paramName);
1443 } else {
1444 this->predistribution = pre_distribution;
1445 }
1446 }
1447 else if(paramName == "perturbation_ratio"){
1448 float perturbation = fromString<float>(getParamVal<std::string>(pe, paramName));
1449 if(perturbation < 0 && perturbation > 1){
1450 throw INVALID(paramName);
1451 } else {
1452 this->perturbation_ratio = perturbation;
1453 }
1454 }
1455
1456
1457 else if(paramName == "proc_load_distributions"){
1458 proc_load_distributions = getParamVal<std::string>(pe, paramName);
1459 this->loadDistSet = true;
1460 }
1461
1462 else if(paramName == "distinct_coordinates"){
1463 distinctDescription = getParamVal<std::string>(pe, paramName);
1464 if(distinctDescription == "T"){
1465 this->distinctCoordSet = true;
1466 } else if(distinctDescription == "F"){
1467 this->distinctCoordSet = true;
1468 } else {
1469 throw "Invalid parameter for distinct_coordinates: " + distinctDescription + ". Candidates: T or F." ;
1470 }
1471 }
1472
1473 else if(paramName == "out_file"){
1474 this->outfile = getParamVal<std::string>(pe, paramName);
1475 }
1476
1477 else {
1478 throw INVALID(paramName);
1479 }
1480 }
1481
1482
1483 if(this->coordinate_dimension == 0){
1484 throw "Dimension must be provided to coordinate generator.";
1485 }
1486 /*
1487 if(this->globalNumCoords == 0){
1488 throw "Number of coordinates must be provided to coordinate generator.";
1489 }
1490 */
1491 /*
1492 if(maxx <= minx ){
1493 throw "Error: maxx= "+ Teuchos::toString<scalar_t>(maxx)+ " and minx=" + Teuchos::toString<scalar_t>(minx);
1494 }
1495 if(maxy <= miny ){
1496 throw "Error: maxy= "+ Teuchos::toString<scalar_t>(maxy)+ " and miny=" + Teuchos::toString<scalar_t>(miny);
1497
1498 }
1499 if(this->dimension == 3 && maxz <= minz ){
1500 throw "Error: maxz= "+ Teuchos::toString<scalar_t>(maxz)+ " and minz=" + Teuchos::toString<scalar_t>(minz);
1501 }
1502 */
1503 if (this->loadDistSet && this->distinctCoordSet){
1504 throw "distinct_coordinates and proc_load_distributions parameters cannot be satisfied together.";
1505 }
1506 this->getHoles(holeDescription);
1507 //this->getDistinctCoordinateDescription(distinctDescription);
1508 this->getProcLoadDistributions(proc_load_distributions);
1509 this->getCoordinateDistributions(coordinate_distributions);
1510 this->getWeightDistribution(numWeightsPerCoord_parameters, this->numWeightsPerCoord);
1511 /*
1512 if(this->numGlobalCoords <= 0){
1513 throw "Must have at least 1 point";
1514 }
1515 */
1516 }
1517 catch(std::string &s){
1518 throw std::runtime_error(s);
1519 }
1520
1521 catch(const char *s){
1522 throw std::runtime_error(s);
1523 }
1524 }
1525public:
1526
1528 if(this->holes){
1529 for (int i = 0; i < this->holeCount; ++i){
1530 delete this->holes[i];
1531 }
1532 free (this->holes);
1533 }
1534
1535
1536 delete []loadDistributions; //sized as the number of processors, the load of each processor.
1537 //delete []distinctCoordinates; // if processors have different or same range for coordinates to be created.
1538 if(coordinateDistributions){
1539
1540 for (int i = 0; i < this->distributionCount; ++i){
1541 delete this->coordinateDistributions[i];
1542 }
1543 free (this->coordinateDistributions);
1544 }
1545 if (this->wd){
1546 for (int i = 0; i < this->numWeightsPerCoord; ++i){
1547 delete this->wd[i];
1548 }
1549 delete []this->wd;
1550 }
1551
1552 if(this->numWeightsPerCoord){
1553 for(int i = 0; i < this->numWeightsPerCoord; ++i)
1554 delete [] this->wghts[i];
1555 delete []this->wghts;
1556 }
1557 if(this->coordinate_dimension){
1558 for(int i = 0; i < this->coordinate_dimension; ++i)
1559 delete [] this->coords[i];
1560 delete [] this->coords;
1561 }
1562 //delete []this->points;
1563 }
1564
1566 std::cout <<"\nGeometric Generator Parameter File Format:" << std::endl;
1567 std::cout <<"- dim=Coordinate Dimension: 2 or 3" << std::endl;
1568 std::cout <<"- Available distributions:" << std::endl;
1569 std::cout <<"\tUNIFORM: -> distribution-1=UNIFORM,NUMCOORDINATES,XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << std::endl;
1570 std::cout <<"\tGRID: -> distribution-2=GRID,XLENGTH,YLENGTH{,ZLENGTH},XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << std::endl;
1571 std::cout <<"\tNORMAL: -> distribution-3=NORMAL,XCENTER,YCENTER{,ZCENTER},XSD,YSD,{,ZSD}" << std::endl;
1572 std::cout <<"- wdim=numWeightsPerCoord: There should be as many weight function as number of weights per coord." << std::endl;
1573 std::cout <<"- Weight Equation: w = (a1 * (x - x1)^b1) + (a2 * (y - y1)^b2) + (a3 * (z - z1)^b3) + c" << std::endl;
1574 std::cout << "Parameter settings:" << std::endl;
1575 std::cout << "\tWeightDistribution-1-a1=a1 " << std::endl;
1576 std::cout << "\tWeightDistribution-1-a2=a2 " << std::endl;
1577 std::cout << "\tWeightDistribution-1-a3=a3 " << std::endl;
1578 std::cout << "\tWeightDistribution-1-b1=b1 " << std::endl;
1579 std::cout << "\tWeightDistribution-1-b2=b2 " << std::endl;
1580 std::cout << "\tWeightDistribution-1-b3=b3 " << std::endl;
1581 std::cout << "\tWeightDistribution-1-x1=x1 " << std::endl;
1582 std::cout << "\tWeightDistribution-1-y1=y1 " << std::endl;
1583 std::cout << "\tWeightDistribution-1-z1=z1 " << std::endl;
1584 std::cout << "\tWeightDistribution-1-c=c " << std::endl;
1585 std::cout << "\tIt is possible to set step function to the result of weight equation." << std::endl;
1586 std::cout << "\tWeightDistribution-1-steps=step1,step2,step3:increasing order" << std::endl;
1587 std::cout << "\tWeightDistribution-1-values=val1,val2,val3,val4." << std::endl;
1588 std::cout << "\t\tIf w < step1 -> w = val1" << std::endl;
1589 std::cout << "\t\tElse if w < step2 -> w = val2" << std::endl;
1590 std::cout << "\t\tElse if w < step3 -> w = val3" << std::endl;
1591 std::cout << "\t\tElse -> w = val4" << std::endl;
1592 std::cout <<"- Holes:" << std::endl;
1593 std::cout << "\thole-1:SPHERE,XCENTER,YCENTER,ZCENTER,RADIUS (only for dim=3)" << std::endl;
1594 std::cout << "\thole-2:CUBE,XCENTER,YCENTER,ZCENTER,EDGE (only for dim=3)" << std::endl;
1595 std::cout << "\thole-3:RECTANGULAR_PRISM,XCENTER,YCENTER,ZCENTER,XEDGE,YEDGE,ZEDGE (only for dim=3)" << std::endl;
1596 std::cout << "\thole-4:SQUARE,XCENTER,YCENTER,EDGE (only for dim=2)" << std::endl;
1597 std::cout << "\thole-5:RECTANGLE,XCENTER,YCENTER,XEDGE,YEDGE (only for dim=2)" << std::endl;
1598 std::cout << "\thole-6:CIRCLE,XCENTER,YCENTER,RADIUS (only for dim=2)" << std::endl;
1599 std::cout << "- out_file:out_file_path : if provided output will be written to files." << std::endl;
1600 std::cout << "- proc_load_distributions:ratio_0, ratio_1, ratio_2....ratio_n. Loads of each processor, should be as many as MPI ranks and should sum up to 1." << std::endl;
1601 std::cout << "- predistribution:distribution_option. Predistribution of the coordinates to the processors. 0 for NONE, 1 RCB, 2 MJ, 3 BLOCK." << std::endl;
1602 std::cout << "- perturbation_ratio:the percent of the local data which will be perturbed in order to simulate the changes in the dynamic partitioning. Float value between 0 and 1." << std::endl;
1603 }
1604
1605 GeometricGenerator(Teuchos::ParameterList &params, const RCP<const Teuchos::Comm<int> > & comm_){
1606 this->wd = NULL;
1607 this->comm = comm_;
1608 this->holes = NULL; //to represent if there is any hole in the input
1609 this->coordinate_dimension = 0; //dimension of the geometry
1610 this->numWeightsPerCoord = 0;
1611 this->worldSize = comm_->getSize(); //comminication world object.
1612 this->numGlobalCoords = 0; //global number of coordinates requested to be created.
1613 this->loadDistributions = NULL; //sized as the number of processors, the load of each processor.
1614 //this->distinctCoordinates = NULL; // if processors have different or same range for coordinates to be created.
1615 this->coordinateDistributions = NULL;
1616 this->holeCount = 0;
1617 this->distributionCount = 0;
1618 this->outfile = "";
1619 this->predistribution = 0;
1620 this->perturbation_ratio = 0;
1621 //this->points = NULL;
1622
1623 /*
1624 this->minx = 0;
1625 this->maxx = 0;
1626 this->miny = 0;
1627 this->maxy = 0;
1628 this->minz = 0;
1629 this->maxz = 0;
1630 */
1631 this->loadDistSet = false;
1632 this->distinctCoordSet = false;
1633 this->myRank = comm_->getRank();
1634
1635 try {
1636 this->parseParams(params);
1637 }
1638 catch(std::string &s){
1639 if(myRank == 0){
1641 }
1642 throw std::runtime_error(s);
1643 }
1644
1645 catch(const char *s){
1646 if(myRank == 0){
1648 }
1649 throw std::runtime_error(s);
1650 }
1651
1652
1653 lno_t myPointCount = 0;
1654 this->numGlobalCoords = 0;
1655
1656 gno_t prefixSum = 0;
1657 for(int i = 0; i < this->distributionCount; ++i){
1658 for(int ii = 0; ii < this->worldSize; ++ii){
1659 lno_t increment = lno_t (this->coordinateDistributions[i]->numPoints * this->loadDistributions[ii]);
1660 if (gno_t(ii) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1661 increment += 1;
1662 }
1663 this->numGlobalCoords += increment;
1664 if(ii < myRank){
1665 prefixSum += increment;
1666 }
1667 }
1668 myPointCount += lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
1669 if (gno_t(myRank) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1670 myPointCount += 1;
1671 }
1672 }
1673
1674 this->coords = new scalar_t *[this->coordinate_dimension];
1675 for(int i = 0; i < this->coordinate_dimension; ++i){
1676 this->coords[i] = new scalar_t[myPointCount];
1677 }
1678
1679 for (int ii = 0; ii < this->coordinate_dimension; ++ii){
1680#ifdef HAVE_ZOLTAN2_OMP
1681#pragma omp parallel for
1682#endif
1683 for(lno_t i = 0; i < myPointCount; ++i){
1684 this->coords[ii][i] = 0;
1685 }
1686 }
1687
1688 this->numLocalCoords = 0;
1689 srand ((myRank + 1) * this->numLocalCoords);
1690 for (int i = 0; i < distributionCount; ++i){
1691
1692 lno_t requestedPointCount = lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
1693 if (gno_t(myRank) < this->coordinateDistributions[i]->numPoints % this->worldSize){
1694 requestedPointCount += 1;
1695 }
1696 //std::cout << "req:" << requestedPointCount << std::endl;
1697 //this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->points + this->numLocalCoords, this->holes, this->holeCount, this->loadDistributions, myRank);
1698 this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->coords, this->numLocalCoords, this->holes, this->holeCount, this->loadDistributions, myRank);
1699 this->numLocalCoords += requestedPointCount;
1700 }
1701
1702 /*
1703
1704 if (this->myRank >= 0){
1705 for(lno_t i = 0; i < this->numLocalCoords; ++i){
1706
1707 std::cout <<"me:" << this->myRank << " "<< this->coords[0][i];
1708 if(this->coordinate_dimension > 1){
1709 std::cout << " " << this->coords[1][i];
1710 }
1711 if(this->coordinate_dimension > 2){
1712 std::cout << " " << this->coords[2][i];
1713 }
1714 std::cout << std::endl;
1715 }
1716 }
1717 */
1718
1719
1720
1721 if (this->predistribution){
1722 redistribute();
1723 }
1724
1725
1726
1727 int scale = 3;
1728 if (this->perturbation_ratio > 0){
1729 this->perturb_data(this->perturbation_ratio, scale);
1730 }
1731 /*
1732 if (this->myRank >= 0){
1733 std::cout << "after distribution" << std::endl;
1734 for(lno_t i = 0; i < this->numLocalCoords; ++i){
1735
1736 std::cout <<"me:" << this->myRank << " " << this->coords[0][i];
1737 if(this->coordinate_dimension > 1){
1738 std::cout << " " << this->coords[1][i];
1739 }
1740 if(this->coordinate_dimension > 2){
1741 std::cout << " " << this->coords[2][i];
1742 }
1743 std::cout << std::endl;
1744 }
1745 }
1746
1747 */
1748
1749
1750 if (this->distinctCoordSet){
1751 //TODO: Partition and migration.
1752 }
1753
1754 if(this->outfile != ""){
1755
1756 std::ofstream myfile;
1757 myfile.open ((this->outfile + Teuchos::toString<int>(myRank)).c_str());
1758 for(lno_t i = 0; i < this->numLocalCoords; ++i){
1759
1760 myfile << this->coords[0][i];
1761 if(this->coordinate_dimension > 1){
1762 myfile << " " << this->coords[1][i];
1763 }
1764 if(this->coordinate_dimension > 2){
1765 myfile << " " << this->coords[2][i];
1766 }
1767 myfile << std::endl;
1768 }
1769 myfile.close();
1770
1771 if (this->myRank == 0){
1772 std::ofstream gnuplotfile("gnu.gnuplot");
1773 for(int i = 0; i < this->worldSize; ++i){
1774 string s = "splot";
1775 if (this->coordinate_dimension == 2){
1776 s = "plot";
1777 }
1778 if (i > 0){
1779 s = "replot";
1780 }
1781 gnuplotfile << s << " \"" << (this->outfile + Teuchos::toString<int>(i)) << "\"" << std::endl;
1782 }
1783 gnuplotfile << "pause -1" << std::endl;
1784 gnuplotfile.close();
1785 }
1786 }
1787
1788
1789
1790 /*
1791 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));
1792
1793 RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmVector2;
1794 Zoltan2::PartitioningSolution< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> > solution;
1795 xmv.applyPartitioningSolution<Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >(this->tmVector, &tmVector2, solution);
1796 */
1797 if (this->numWeightsPerCoord > 0){
1798 this->wghts = new scalar_t *[this->numWeightsPerCoord];
1799 for(int i = 0; i < this->numWeightsPerCoord; ++i){
1800 this->wghts[i] = new scalar_t[this->numLocalCoords];
1801 }
1802 }
1803
1804 for(int ii = 0; ii < this->numWeightsPerCoord; ++ii){
1805 switch(this->coordinate_dimension){
1806 case 1:
1807 {
1808#ifdef HAVE_ZOLTAN2_OMP
1809#pragma omp parallel for
1810#endif
1811 for (lno_t i = 0; i < this->numLocalCoords; ++i){
1812 this->wghts[ii][i] = this->wd[ii]->get1DWeight(this->coords[0][i]);
1813 }
1814 }
1815 break;
1816 case 2:
1817 {
1818#ifdef HAVE_ZOLTAN2_OMP
1819#pragma omp parallel for
1820#endif
1821 for (lno_t i = 0; i < this->numLocalCoords; ++i){
1822 this->wghts[ii][i] = this->wd[ii]->get2DWeight(this->coords[0][i], this->coords[1][i]);
1823 }
1824 }
1825 break;
1826 case 3:
1827 {
1828#ifdef HAVE_ZOLTAN2_OMP
1829#pragma omp parallel for
1830#endif
1831 for (lno_t i = 0; i < this->numLocalCoords; ++i){
1832 this->wghts[ii][i] = this->wd[ii]->get3DWeight(this->coords[0][i], this->coords[1][i], this->coords[2][i]);
1833 }
1834 }
1835 break;
1836 }
1837 }
1838 }
1839
1840 //############################################################//
1842 //############################################################//
1843 void perturb_data(float used_perturbation_ratio, int scale){
1844 scalar_t *maxCoords= new scalar_t [this->coordinate_dimension];
1845 scalar_t *minCoords= new scalar_t [this->coordinate_dimension];
1846 for (int dim = 0; dim < this->coordinate_dimension; ++dim){
1847 minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
1848 for (lno_t i = 1; i < this->numLocalCoords; ++i){
1849 if (minCoords[dim] > this->coords[dim][i]){
1850 minCoords[dim] = this->coords[dim][i];
1851 }
1852
1853 if (maxCoords[dim] < this->coords[dim][i]){
1854 maxCoords[dim] = this->coords[dim][i];
1855 }
1856 }
1857
1858
1859
1860
1861 scalar_t center = (maxCoords[dim] + minCoords[dim]) / 2;
1862
1863 minCoords[dim] = center - (center - minCoords[dim]) * scale;
1864 maxCoords[dim] = (maxCoords[dim] - center) * scale + center;
1865
1866 }
1867
1868 gno_t numLocalToPerturb = gno_t(this->numLocalCoords*used_perturbation_ratio);
1869 //std::cout << "numLocalToPerturb :" << numLocalToPerturb << std::endl;
1870 for (int dim = 0; dim < this->coordinate_dimension; ++dim){
1871 scalar_t range = maxCoords[dim] - minCoords[dim];
1872 for (gno_t i = 0; i < numLocalToPerturb; ++i){
1873 this->coords[dim][i] = (rand() / double (RAND_MAX)) * (range) + minCoords[dim];
1874
1875 }
1876 }
1877 delete []maxCoords;
1878 delete []minCoords;
1879 }
1880
1881 //############################################################//
1883 //############################################################//
1884
1885 //Returns the partitioning dimension as even as possible
1886 void getBestSurface (int remaining, int *dimProcs, int dim, int currentDim, int &bestSurface, int *bestDimProcs){
1887
1888 if (currentDim < dim - 1){
1889 int ipx = 1;
1890 while(ipx <= remaining) {
1891 if(remaining % ipx == 0) {
1892 int nremain = remaining / ipx;
1893 dimProcs[currentDim] = ipx;
1894 getBestSurface (nremain, dimProcs, dim, currentDim + 1, bestSurface, bestDimProcs);
1895 }
1896 ipx++;
1897 }
1898 }
1899 else {
1900 dimProcs[currentDim] = remaining;
1901 int surface = 0;
1902 for (int i = 0; i < dim; ++i) surface += dimProcs[i];
1903 if (surface < bestSurface){
1904 bestSurface = surface;
1905 for (int i = 0; i < dim; ++i) bestDimProcs[i] = dimProcs[i];
1906 }
1907 }
1908
1909 }
1910
1911 //returns min and max coordinates along each dimension
1912 void getMinMaxCoords(scalar_t *globalMaxCoords, scalar_t *globalMinCoords){
1913 scalar_t *maxCoords= new scalar_t [this->coordinate_dimension];
1914 scalar_t *minCoords= new scalar_t [this->coordinate_dimension];
1915 for (int dim = 0; dim < this->coordinate_dimension; ++dim){
1916 minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
1917 for (lno_t i = 1; i < this->numLocalCoords; ++i){
1918 if (minCoords[dim] > this->coords[dim][i]){
1919 minCoords[dim] = this->coords[dim][i];
1920 }
1921
1922 if (maxCoords[dim] < this->coords[dim][i]){
1923 maxCoords[dim] = this->coords[dim][i];
1924 }
1925 }
1926 }
1927
1928 reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MAX,
1929 this->coordinate_dimension,
1930 maxCoords,
1931 globalMaxCoords);
1932
1933
1934 reduceAll<int, scalar_t>( *(this->comm), Teuchos::REDUCE_MIN,
1935 this->coordinate_dimension,
1936 minCoords,
1937 globalMinCoords);
1938
1939 delete []maxCoords;
1940 delete []minCoords;
1941 }
1942
1943
1944 //performs a block partitioning.
1945 //then distributes the points of the overloaded parts to underloaded parts.
1946 void blockPartition(int *coordinate_grid_parts){
1947
1948#ifdef _MSC_VER
1949 typedef SSIZE_T ssize_t;
1950#endif
1951
1952 //############################################################//
1954 //############################################################//
1955
1956 scalar_t *maxCoords= new scalar_t [this->coordinate_dimension];
1957 scalar_t *minCoords= new scalar_t [this->coordinate_dimension];
1958 //global min and max coordinates in each dimension.
1959 this->getMinMaxCoords(maxCoords, minCoords);
1960
1961
1962 //############################################################//
1964 //############################################################//
1965 int remaining = this->worldSize;
1966 int coord_dim = this->coordinate_dimension;
1967 int *dimProcs = new int[coord_dim];
1968 int *bestDimProcs = new int[coord_dim];
1969 int currentDim = 0;
1970
1971 int bestSurface = 0;
1972 dimProcs[0] = remaining;
1973 for (int i = 1; i < coord_dim; ++i) dimProcs[i] = 1;
1974 for (int i = 0; i < coord_dim; ++i) bestSurface += dimProcs[i];
1975 for (int i = 0; i < coord_dim; ++i) bestDimProcs[i] = dimProcs[i];
1976 //divides the parts into dimensions as even as possible.
1977 getBestSurface ( remaining, dimProcs, coord_dim, currentDim, bestSurface, bestDimProcs);
1978
1979
1980 delete []dimProcs;
1981
1982 //############################################################//
1984 //############################################################//
1985 int *shiftProcCount = new int[coord_dim];
1986 //how the consecutive parts along a dimension
1987 //differs in the index.
1988
1989 int remainingProc = this->worldSize;
1990 for (int dim = 0; dim < coord_dim; ++dim){
1991 remainingProc = remainingProc / bestDimProcs[dim];
1992 shiftProcCount[dim] = remainingProc;
1993 }
1994
1995 scalar_t *dim_slices = new scalar_t[coord_dim];
1996 for (int dim = 0; dim < coord_dim; ++dim){
1997 scalar_t dim_range = maxCoords[dim] - minCoords[dim];
1998 scalar_t slice = dim_range / bestDimProcs[dim];
1999 dim_slices[dim] = slice;
2000 }
2001
2002 //############################################################//
2004 //############################################################//
2005
2006 gno_t *numPointsInParts = new gno_t[this->worldSize];
2007 gno_t *numGlobalPointsInParts = new gno_t[this->worldSize];
2008 gno_t *numPointsInPartsInclusiveUptoMyIndex = new gno_t[this->worldSize];
2009
2010 gno_t *partBegins = new gno_t [this->worldSize];
2011 gno_t *partNext = new gno_t[this->numLocalCoords];
2012
2013
2014 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2015 partNext[i] = -1;
2016 }
2017 for (int i = 0; i < this->worldSize; ++i) {
2018 partBegins[i] = 1;
2019 }
2020
2021 for (int i = 0; i < this->worldSize; ++i)
2022 numPointsInParts[i] = 0;
2023
2024 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2025 int partIndex = 0;
2026 for (int dim = 0; dim < coord_dim; ++dim){
2027 int shift = int ((this->coords[dim][i] - minCoords[dim]) / dim_slices[dim]);
2028 if (shift >= bestDimProcs[dim]){
2029 shift = bestDimProcs[dim] - 1;
2030 }
2031 shift = shift * shiftProcCount[dim];
2032 partIndex += shift;
2033 }
2034 numPointsInParts[partIndex] += 1;
2035 coordinate_grid_parts[i] = partIndex;
2036
2037 partNext[i] = partBegins[partIndex];
2038 partBegins[partIndex] = i;
2039
2040 }
2041
2042 //############################################################//
2044 //############################################################//
2045 reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2046 this->worldSize,
2047 numPointsInParts,
2048 numGlobalPointsInParts);
2049
2050
2051 Teuchos::scan<int,gno_t>(
2052 *(this->comm), Teuchos::REDUCE_SUM,
2053 this->worldSize,
2054 numPointsInParts,
2055 numPointsInPartsInclusiveUptoMyIndex
2056 );
2057
2058
2059
2060
2061 /*
2062 gno_t totalSize = 0;
2063 for (int i = 0; i < this->worldSize; ++i){
2064 totalSize += numPointsInParts[i];
2065 }
2066 if (totalSize != this->numLocalCoords){
2067 std::cout << "me:" << this->myRank << " ts:" << totalSize << " nl:" << this->numLocalCoords << std::endl;
2068 }
2069 */
2070
2071
2072 //std::cout << "me:" << this->myRank << " ilk print" << std::endl;
2073
2074 gno_t optimal_num = gno_t(this->numGlobalCoords/double(this->worldSize)+0.5);
2075#ifdef printparts
2076 if (this->myRank == 0){
2077 gno_t totalSize = 0;
2078 for (int i = 0; i < this->worldSize; ++i){
2079 std::cout << "me:" << this->myRank << " NumPoints in part:" << i << " is: " << numGlobalPointsInParts[i] << std::endl;
2080 totalSize += numGlobalPointsInParts[i];
2081 }
2082 std::cout << "Total:" << totalSize << " ng:" << this->numGlobalCoords << std::endl;
2083 std::cout << "optimal_num:" << optimal_num << std::endl;
2084 }
2085#endif
2086 ssize_t *extraInPart = new ssize_t [this->worldSize];
2087
2088 ssize_t extraExchange = 0;
2089 for (int i = 0; i < this->worldSize; ++i){
2090 extraInPart[i] = numGlobalPointsInParts[i] - optimal_num;
2091 extraExchange += extraInPart[i];
2092 }
2093 if (extraExchange != 0){
2094 int addition = -1;
2095 if (extraExchange < 0) addition = 1;
2096 for (ssize_t i = 0; i < extraExchange; ++i){
2097 extraInPart[i % this->worldSize] += addition;
2098 }
2099 }
2100
2101 //############################################################//
2103 //############################################################//
2104
2105 int overloadedPartCount = 0;
2106 int *overloadedPartIndices = new int [this->worldSize];
2107
2108
2109 int underloadedPartCount = 0;
2110 int *underloadedPartIndices = new int [this->worldSize];
2111
2112 for (int i = 0; i < this->worldSize; ++i){
2113 if(extraInPart[i] > 0){
2114 overloadedPartIndices[overloadedPartCount++] = i;
2115 }
2116 else if(extraInPart[i] < 0){
2117 underloadedPartIndices[underloadedPartCount++] = i;
2118 }
2119 }
2120
2121 int underloadpartindex = underloadedPartCount - 1;
2122
2123
2124 int numPartsISendFrom = 0;
2125 int *mySendFromParts = new int[this->worldSize * 2];
2126 gno_t *mySendFromPartsCounts = new gno_t[this->worldSize * 2];
2127
2128 int numPartsISendTo = 0;
2129 int *mySendParts = new int[this->worldSize * 2];
2130 gno_t *mySendCountsToParts = new gno_t[this->worldSize * 2];
2131
2132
2133 //############################################################//
2138 //############################################################//
2139 for (int i = overloadedPartCount - 1; i >= 0; --i){
2140 //get the overloaded part
2141 //the overload
2142 int overloadPart = overloadedPartIndices[i];
2143 gno_t overload = extraInPart[overloadPart];
2144 gno_t myload = numPointsInParts[overloadPart];
2145
2146
2147 //the inclusive load of the processors up to me
2148 gno_t inclusiveLoadUpToMe = numPointsInPartsInclusiveUptoMyIndex[overloadPart];
2149
2150 //the exclusive load of the processors up to me
2151 gno_t exclusiveLoadUptoMe = inclusiveLoadUpToMe - myload;
2152
2153
2154 if (exclusiveLoadUptoMe >= overload){
2155 //this processor does not have to convert anything.
2156 //for this overloaded part.
2157 //set the extra for this processor to zero.
2158 overloadedPartIndices[i] = -1;
2159 extraInPart[overloadPart] = 0;
2160 //then consume underloaded parts.
2161 while (overload > 0){
2162 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2163 ssize_t underload = extraInPart[nextUnderLoadedPart];
2164 ssize_t left = overload + underload;
2165
2166 if(left >= 0){
2167 extraInPart[nextUnderLoadedPart] = 0;
2168 underloadedPartIndices[underloadpartindex--] = -1;
2169 }
2170 else {
2171 extraInPart[nextUnderLoadedPart] = left;
2172 }
2173 overload = left;
2174 }
2175 }
2176 else if (exclusiveLoadUptoMe < overload){
2177 //if the previous processors load is not enough
2178 //then this processor should convert some of its elements.
2179 gno_t mySendCount = overload - exclusiveLoadUptoMe;
2180 //how much more needed.
2181 gno_t sendAfterMe = 0;
2182 //if my load is not enough
2183 //then the next processor will continue to convert.
2184 if (mySendCount > myload){
2185 sendAfterMe = mySendCount - myload;
2186 mySendCount = myload;
2187 }
2188
2189
2190 //this processors will convert from overloaded part
2191 //as many as mySendCount items.
2192 mySendFromParts[numPartsISendFrom] = overloadPart;
2193 mySendFromPartsCounts[numPartsISendFrom++] = mySendCount;
2194
2195 //first consume underloaded parts for the previous processors.
2196 while (exclusiveLoadUptoMe > 0){
2197 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2198 ssize_t underload = extraInPart[nextUnderLoadedPart];
2199 ssize_t left = exclusiveLoadUptoMe + underload;
2200
2201 if(left >= 0){
2202 extraInPart[nextUnderLoadedPart] = 0;
2203 underloadedPartIndices[underloadpartindex--] = -1;
2204 }
2205 else {
2206 extraInPart[nextUnderLoadedPart] = left;
2207 }
2208 exclusiveLoadUptoMe = left;
2209 }
2210
2211 //consume underloaded parts for my load.
2212 while (mySendCount > 0){
2213 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2214 ssize_t underload = extraInPart[nextUnderLoadedPart];
2215 ssize_t left = mySendCount + underload;
2216
2217 if(left >= 0){
2218 mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2219 mySendCountsToParts[numPartsISendTo++] = -underload;
2220
2221 extraInPart[nextUnderLoadedPart] = 0;
2222 underloadedPartIndices[underloadpartindex--] = -1;
2223 }
2224 else {
2225 extraInPart[nextUnderLoadedPart] = left;
2226
2227 mySendParts[numPartsISendTo] = nextUnderLoadedPart;
2228 mySendCountsToParts[numPartsISendTo++] = mySendCount;
2229
2230 }
2231 mySendCount = left;
2232 }
2233 //consume underloaded parts for the load of the processors after my index.
2234 while (sendAfterMe > 0){
2235 int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
2236 ssize_t underload = extraInPart[nextUnderLoadedPart];
2237 ssize_t left = sendAfterMe + underload;
2238
2239 if(left >= 0){
2240 extraInPart[nextUnderLoadedPart] = 0;
2241 underloadedPartIndices[underloadpartindex--] = -1;
2242 }
2243 else {
2244 extraInPart[nextUnderLoadedPart] = left;
2245 }
2246 sendAfterMe = left;
2247 }
2248 }
2249 }
2250
2251
2252 //############################################################//
2254 //############################################################//
2255 for (int i = 0 ; i < numPartsISendFrom; ++i){
2256
2257 //get the part from which the elements will be converted.
2258 int sendFromPart = mySendFromParts[i];
2259 ssize_t sendCount = mySendFromPartsCounts[i];
2260 while(sendCount > 0){
2261 int partToSendIndex = numPartsISendTo - 1;
2262 int partToSend = mySendParts[partToSendIndex];
2263
2264 int sendCountToThisPart = mySendCountsToParts[partToSendIndex];
2265
2266 //determine which part i should convert to
2267 //and how many to this part.
2268 if (sendCountToThisPart <= sendCount){
2269 mySendParts[partToSendIndex] = 0;
2270 mySendCountsToParts[partToSendIndex] = 0;
2271 --numPartsISendTo;
2272 sendCount -= sendCountToThisPart;
2273 }
2274 else {
2275 mySendCountsToParts[partToSendIndex] = sendCountToThisPart - sendCount;
2276 sendCountToThisPart = sendCount;
2277 sendCount = 0;
2278 }
2279
2280
2281 gno_t toChange = partBegins[sendFromPart];
2282 gno_t previous_begin = partBegins[partToSend];
2283
2284 //do the conversion.
2285 for (int k = 0; k < sendCountToThisPart - 1; ++k){
2286 coordinate_grid_parts[toChange] = partToSend;
2287 toChange = partNext[toChange];
2288 }
2289 coordinate_grid_parts[toChange] = partToSend;
2290
2291 gno_t newBegin = partNext[toChange];
2292 partNext[toChange] = previous_begin;
2293 partBegins[partToSend] = partBegins[sendFromPart];
2294 partBegins[sendFromPart] = newBegin;
2295 }
2296 }
2297
2298 //if (this->myRank == 0) std::cout << "4" << std::endl;
2299
2300#ifdef printparts
2301
2302
2303 for (int i = 0; i < this->worldSize; ++i) numPointsInParts[i] = 0;
2304
2305 for (int i = 0; i < this->numLocalCoords; ++i){
2306 numPointsInParts[coordinate_grid_parts[i]] += 1;
2307 }
2308
2309 reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
2310 this->worldSize,
2311 numPointsInParts,
2312 numGlobalPointsInParts);
2313 if (this->myRank == 0){
2314 std::cout << "reassigning" << std::endl;
2315 gno_t totalSize = 0;
2316 for (int i = 0; i < this->worldSize; ++i){
2317 std::cout << "me:" << this->myRank << " NumPoints in part:" << i << " is: " << numGlobalPointsInParts[i] << std::endl;
2318 totalSize += numGlobalPointsInParts[i];
2319 }
2320 std::cout << "Total:" << totalSize << " ng:" << this->numGlobalCoords << std::endl;
2321 }
2322#endif
2323 delete []mySendCountsToParts;
2324 delete []mySendParts;
2325 delete []mySendFromPartsCounts;
2326 delete []mySendFromParts;
2327 delete []underloadedPartIndices;
2328 delete []overloadedPartIndices;
2329 delete []extraInPart;
2330 delete []partNext;
2331 delete []partBegins;
2332 delete []numPointsInPartsInclusiveUptoMyIndex;
2333 delete []numPointsInParts;
2334 delete []numGlobalPointsInParts;
2335
2336 delete []shiftProcCount;
2337 delete []bestDimProcs;
2338 delete []dim_slices;
2339 delete []minCoords;
2340 delete []maxCoords;
2341 }
2342
2343 //given the part numbers for each local coordinate,
2344 //distributes the coordinates to the corresponding processors.
2345 void distribute_points(int *coordinate_grid_parts){
2346
2347 Tpetra::Distributor distributor(comm);
2348 ArrayView<const int> pIds( coordinate_grid_parts, this->numLocalCoords);
2349 gno_t numMyNewGnos = distributor.createFromSends(pIds);
2350
2351
2352 Kokkos::View<scalar_t*, Kokkos::HostSpace> recvBuf2(
2353 Kokkos::ViewAllocateWithoutInitializing("recvBuf2"),
2354 numMyNewGnos);
2355
2356 for (int i = 0; i < this->coordinate_dimension; ++i){
2357 Kokkos::View<scalar_t*, Kokkos::HostSpace> s;
2358 if (this->numLocalCoords > 0)
2359 s = Kokkos::View<scalar_t *, Kokkos::HostSpace>(
2360 this->coords[i], this->numLocalCoords); //unmanaged
2361
2362 distributor.doPostsAndWaits(s, 1, recvBuf2);
2363
2364 delete [] this->coords[i];
2365 this->coords[i] = new scalar_t[numMyNewGnos];
2366 for (lno_t j = 0; j < numMyNewGnos; ++j){
2367 this->coords[i][j] = recvBuf2[j];
2368 }
2369
2370 }
2371 this->numLocalCoords = numMyNewGnos;
2372 }
2373
2374 //calls MJ for p = numProcs
2375 int predistributeMJ(int *coordinate_grid_parts){
2376 int coord_dim = this->coordinate_dimension;
2377
2378 lno_t numLocalPoints = this->numLocalCoords;
2379 gno_t numGlobalPoints = this->numGlobalCoords;
2380
2381
2382 //T **weight = NULL;
2383 //typedef Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> tMVector_t;
2384 RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
2385 new Tpetra::Map<lno_t, gno_t, node_t> (numGlobalPoints, numLocalPoints, 0, comm));
2386
2387 Teuchos::Array<Teuchos::ArrayView<const scalar_t> > coordView(coord_dim);
2388
2389
2390
2391 for (int i=0; i < coord_dim; i++){
2392 if(numLocalPoints > 0){
2393 Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalPoints);
2394 coordView[i] = a;
2395 } else{
2396 Teuchos::ArrayView<const scalar_t> a;
2397 coordView[i] = a;
2398 }
2399 }
2400
2401 RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >tmVector = RCP< Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> >(
2402 new Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t>( mp, coordView.view(0, coord_dim), coord_dim));
2403
2404
2405 RCP<const tMVector_t> coordsConst = Teuchos::rcp_const_cast<const tMVector_t>(tmVector);
2406 std::vector<const scalar_t *> weights;
2407 std::vector <int> stride;
2408
2409 typedef Zoltan2::XpetraMultiVectorAdapter<tMVector_t> inputAdapter_t;
2410 //inputAdapter_t ia(coordsConst);
2411 inputAdapter_t ia(coordsConst,weights, stride);
2412
2413 Teuchos::RCP <Teuchos::ParameterList> params ;
2414 params =RCP <Teuchos::ParameterList> (new Teuchos::ParameterList, true);
2415
2416
2417 params->set("algorithm", "multijagged");
2418 params->set("num_global_parts", this->worldSize);
2419
2420 //TODO we need to fix the setting parts.
2421 //Although MJ sets the parts with
2422 //currently the part setting is not correct when migration is done.
2423 //params->set("migration_check_option", 2);
2424
2425
2427
2428
2429 try {
2430#ifdef HAVE_ZOLTAN2_MPI
2431 problem = new Zoltan2::PartitioningProblem<inputAdapter_t>(&ia, params.getRawPtr(),
2432 MPI_COMM_WORLD);
2433#else
2434 problem = new Zoltan2::PartitioningProblem<inputAdapter_t>(&ia, params.getRawPtr());
2435#endif
2436 }
2437 CATCH_EXCEPTIONS("PartitioningProblem()")
2438
2439 try {
2440 problem->solve();
2441 }
2442 CATCH_EXCEPTIONS("solve()")
2443
2444 const typename inputAdapter_t::part_t *partIds = problem->getSolution().getPartListView();
2445
2446 for (lno_t i = 0; i < this->numLocalCoords;++i){
2447 coordinate_grid_parts[i] = partIds[i];
2448 //std::cout << "me:" << this->myRank << " i:" << i << " goes to:" << partIds[i] << std::endl;
2449 }
2450 delete problem;
2451 return 0;
2452 }
2453
2454 //calls RCP for p = numProcs
2455 int predistributeRCB(int *coordinate_grid_parts){
2456 int rank = this->myRank;
2457 int nprocs = this->worldSize;
2458 DOTS<tMVector_t> dots_;
2459
2460 MEMORY_CHECK(rank==0 || rank==nprocs-1, "After initializing MPI");
2461
2462
2463 int nWeights = 0;
2464 int debugLevel=0;
2465 string memoryOn("memoryOn");
2466 string memoryOff("memoryOff");
2467 bool doMemory=false;
2468 int numGlobalParts = nprocs;
2469 int dummyTimer=0;
2470 bool remap=0;
2471
2472 string balanceCount("balance_object_count");
2473 string balanceWeight("balance_object_weight");
2474 string mcnorm1("multicriteria_minimize_total_weight");
2475 string mcnorm2("multicriteria_balance_total_maximum");
2476 string mcnorm3("multicriteria_minimize_maximum_weight");
2477 string objective(balanceWeight); // default
2478
2479 // Process command line input
2480 CommandLineProcessor commandLine(false, true);
2481 //commandLine.setOption("size", &numGlobalCoords,
2482 // "Approximate number of global coordinates.");
2483 int input_option = 0;
2484 commandLine.setOption("input_option", &input_option,
2485 "whether to use mesh creation, geometric generator, or file input");
2486 string inputFile = "";
2487
2488 commandLine.setOption("input_file", &inputFile,
2489 "the input file for geometric generator or file input");
2490
2491
2492 commandLine.setOption("size", &numGlobalCoords,
2493 "Approximate number of global coordinates.");
2494 commandLine.setOption("numParts", &numGlobalParts,
2495 "Number of parts (default is one per proc).");
2496 commandLine.setOption("nWeights", &nWeights,
2497 "Number of weights per coordinate, zero implies uniform weights.");
2498 commandLine.setOption("debug", &debugLevel, "Zoltan1 debug level");
2499 commandLine.setOption("remap", "no-remap", &remap,
2500 "Zoltan1 REMAP parameter; disabled by default for scalability testing");
2501 commandLine.setOption("timers", &dummyTimer, "ignored");
2502 commandLine.setOption(memoryOn.c_str(), memoryOff.c_str(), &doMemory,
2503 "do memory profiling");
2504
2505 string doc(balanceCount);
2506 doc.append(": ignore weights\n");
2507 doc.append(balanceWeight);
2508 doc.append(": balance on first weight\n");
2509 doc.append(mcnorm1);
2510 doc.append(": given multiple weights, balance their total.\n");
2511 doc.append(mcnorm3);
2512 doc.append(": given multiple weights, "
2513 "balance the maximum for each coordinate.\n");
2514 doc.append(mcnorm2);
2515 doc.append(": given multiple weights, balance the L2 norm of the weights.\n");
2516 commandLine.setOption("objective", &objective, doc.c_str());
2517
2518 CommandLineProcessor::EParseCommandLineReturn rc =
2519 commandLine.parse(0, NULL);
2520
2521
2522
2523 if (rc != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
2524 if (rc == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) {
2525 if (rank==0) std::cout << "PASS" << std::endl;
2526 return 1;
2527 }
2528 else {
2529 if (rank==0) std::cout << "FAIL" << std::endl;
2530 return 0;
2531 }
2532 }
2533
2534 //MEMORY_CHECK(doMemory && rank==0, "After processing parameters");
2535
2536 // Create the data structure
2537
2538 int coord_dim = this->coordinate_dimension;
2539
2540
2541 RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
2542 new Tpetra::Map<lno_t, gno_t, node_t> (this->numGlobalCoords, this->numLocalCoords, 0, this->comm));
2543
2544 Teuchos::Array<Teuchos::ArrayView<const scalar_t> > coordView(coord_dim);
2545 for (int i=0; i < coord_dim; i++){
2546 if(numLocalCoords > 0){
2547 Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalCoords);
2548 coordView[i] = a;
2549 } else{
2550 Teuchos::ArrayView<const scalar_t> a;
2551 coordView[i] = a;
2552 }
2553 }
2554
2555 tMVector_t *tmVector = new tMVector_t( mp, coordView.view(0, coord_dim), coord_dim);
2556
2557 dots_.coordinates = tmVector;
2558 dots_.weights.resize(nWeights);
2559
2560
2561 MEMORY_CHECK(doMemory && rank==0, "After creating input");
2562
2563 // Now call Zoltan to partition the problem.
2564
2565 float ver;
2566 int aok = Zoltan_Initialize(0,NULL, &ver);
2567
2568 if (aok != 0){
2569 printf("Zoltan_Initialize failed\n");
2570 exit(0);
2571 }
2572
2573 struct Zoltan_Struct *zz;
2574 zz = Zoltan_Create(MPI_COMM_WORLD);
2575
2576 Zoltan_Set_Param(zz, "LB_METHOD", "RCB");
2577 Zoltan_Set_Param(zz, "LB_APPROACH", "PARTITION");
2578 Zoltan_Set_Param(zz, "CHECK_GEOM", "0");
2579 Zoltan_Set_Param(zz, "NUM_GID_ENTRIES", "1");
2580 Zoltan_Set_Param(zz, "NUM_LID_ENTRIES", "0");
2581 Zoltan_Set_Param(zz, "RETURN_LISTS", "PART");
2582 std::ostringstream oss;
2583 oss << numGlobalParts;
2584 Zoltan_Set_Param(zz, "NUM_GLOBAL_PARTS", oss.str().c_str());
2585 oss.str("");
2586 oss << debugLevel;
2587 Zoltan_Set_Param(zz, "DEBUG_LEVEL", oss.str().c_str());
2588
2589 if (remap)
2590 Zoltan_Set_Param(zz, "REMAP", "1");
2591 else
2592 Zoltan_Set_Param(zz, "REMAP", "0");
2593
2594 if (objective != balanceCount){
2595 oss.str("");
2596 oss << nWeights;
2597 Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", oss.str().c_str());
2598
2599 if (objective == mcnorm1)
2600 Zoltan_Set_Param(zz, "RCB_MULTICRITERIA_NORM", "1");
2601 else if (objective == mcnorm2)
2602 Zoltan_Set_Param(zz, "RCB_MULTICRITERIA_NORM", "2");
2603 else if (objective == mcnorm3)
2604 Zoltan_Set_Param(zz, "RCB_MULTICRITERIA_NORM", "3");
2605 }
2606 else{
2607 Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", "0");
2608 }
2609
2610 Zoltan_Set_Num_Obj_Fn(zz, getNumObj<tMVector_t>, &dots_);
2611 Zoltan_Set_Obj_List_Fn(zz, getObjList<tMVector_t>, &dots_);
2612 Zoltan_Set_Num_Geom_Fn(zz, getDim<tMVector_t>, &dots_);
2613 Zoltan_Set_Geom_Multi_Fn(zz, getCoords<tMVector_t>, &dots_);
2614
2615 int changes, numGidEntries, numLidEntries, numImport, numExport;
2616 ZOLTAN_ID_PTR importGlobalGids, importLocalGids;
2617 ZOLTAN_ID_PTR exportGlobalGids, exportLocalGids;
2618 int *importProcs, *importToPart, *exportProcs, *exportToPart;
2619
2620 MEMORY_CHECK(doMemory && rank==0, "Before Zoltan_LB_Partition");
2621
2622
2623 aok = Zoltan_LB_Partition(zz, &changes, &numGidEntries, &numLidEntries,
2624 &numImport, &importGlobalGids, &importLocalGids,
2625 &importProcs, &importToPart,
2626 &numExport, &exportGlobalGids, &exportLocalGids,
2627 &exportProcs, &exportToPart);
2628
2629
2630 MEMORY_CHECK(doMemory && rank==0, "After Zoltan_LB_Partition");
2631
2632 for (lno_t i = 0; i < numLocalCoords; i++)
2633 coordinate_grid_parts[i] = exportToPart[i];
2634 Zoltan_Destroy(&zz);
2635 MEMORY_CHECK(doMemory && rank==0, "After Zoltan_Destroy");
2636
2637 delete dots_.coordinates;
2638 return 0;
2639}
2641 int *coordinate_grid_parts = new int[this->numLocalCoords];
2642 switch (this->predistribution){
2643 case 1:
2644 this->predistributeRCB(coordinate_grid_parts);
2645 break;
2646 case 2:
2647
2648 this->predistributeMJ(coordinate_grid_parts);
2649 break;
2650 case 3:
2651 //block
2652 blockPartition(coordinate_grid_parts);
2653 break;
2654 }
2655 this->distribute_points(coordinate_grid_parts);
2656
2657 delete []coordinate_grid_parts;
2658
2659
2660 }
2661
2662 //############################################################//
2664 //############################################################//
2665
2666
2668 return this->numWeightsPerCoord;
2669 }
2671 return this->coordinate_dimension;
2672 }
2674 return this->numLocalCoords;
2675 }
2677 return this->numGlobalCoords;
2678 }
2679
2681 return this->coords;
2682 }
2683
2685 return this->wghts;
2686 }
2687
2688 void getLocalCoordinatesCopy( scalar_t ** c){
2689 for(int ii = 0; ii < this->coordinate_dimension; ++ii){
2690#ifdef HAVE_ZOLTAN2_OMP
2691#pragma omp parallel for
2692#endif
2693 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2694 c[ii][i] = this->coords[ii][i];
2695 }
2696 }
2697 }
2698
2699 void getLocalWeightsCopy(scalar_t **w){
2700 for(int ii = 0; ii < this->numWeightsPerCoord; ++ii){
2701#ifdef HAVE_ZOLTAN2_OMP
2702#pragma omp parallel for
2703#endif
2704 for (lno_t i = 0; i < this->numLocalCoords; ++i){
2705 w[ii][i] = this->wghts[ii][i];
2706 }
2707 }
2708 }
2709};
2710}
2711
2712#endif /* GEOMETRICGENERATOR */
#define INVALID_SHAPE_ARG(SHAPE, REQUIRED)
#define MAX_ITER_ALLOWED
#define INVALIDSHAPE(STR, DIM)
#define INVALID(STR)
#define CATCH_EXCEPTIONS(pp)
#define MAX_WEIGHT_DIM
Defines the PartitioningProblem class.
Defines the PartitioningSolution class.
#define MEMORY_CHECK(iPrint, msg)
Defines the XpetraMultiVectorAdapter.
virtual bool isInArea(CoordinatePoint< T > dot)
CircleHole(CoordinatePoint< T > center_, T edge_)
CoordinateDistribution(gno_t np_, int dim, int wSize)
void GetPoints(lno_t requestedPointcount, T **coords, lno_t tindex, Hole< T > **holes, lno_t holeCount, float *sharedRatios_, int myRank)
void GetPoints(lno_t requestedPointcount, CoordinatePoint< T > *points, Hole< T > **holes, lno_t holeCount, float *sharedRatios_, int myRank)
virtual CoordinatePoint< T > getPoint(gno_t point_index, unsigned int &state)=0
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 CoordinatePoint< T > getPoint(gno_t pindex, unsigned int &state)
CoordinateNormalDistribution(gno_t np_, int dim, CoordinatePoint< T > center_, T sd_x, T sd_y, T sd_z, int wSize)
virtual CoordinatePoint< T > getPoint(gno_t pindex, unsigned int &state)
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)
virtual CoordinatePoint< T > getPoint(gno_t pindex, unsigned int &state)
CubeHole(CoordinatePoint< T > center_, T edge_)
virtual bool isInArea(CoordinatePoint< T > dot)
std::vector< std::vector< float > > weights
int predistributeMJ(int *coordinate_grid_parts)
void distribute_points(int *coordinate_grid_parts)
int predistributeRCB(int *coordinate_grid_parts)
void perturb_data(float used_perturbation_ratio, int scale)
void getBestSurface(int remaining, int *dimProcs, int dim, int currentDim, int &bestSurface, int *bestDimProcs)
GeometricGenerator(Teuchos::ParameterList &params, const RCP< const Teuchos::Comm< int > > &comm_)
void getMinMaxCoords(scalar_t *globalMaxCoords, scalar_t *globalMinCoords)
void blockPartition(int *coordinate_grid_parts)
virtual bool isInArea(CoordinatePoint< T > dot)=0
CoordinatePoint< T > center
Hole(CoordinatePoint< T > center_, T edge1_, T edge2_, T edge3_)
virtual bool isInArea(CoordinatePoint< T > dot)
RectangleHole(CoordinatePoint< T > center_, T edge_x_, T edge_y_)
RectangularPrismHole(CoordinatePoint< T > center_, T edge_x_, T edge_y_, T edge_z_)
virtual bool isInArea(CoordinatePoint< T > dot)
virtual bool isInArea(CoordinatePoint< T > dot)
SphereHole(CoordinatePoint< T > center_, T edge_)
SquareHole(CoordinatePoint< T > center_, T edge_)
virtual bool isInArea(CoordinatePoint< T > dot)
Expression is a generic following method.
virtual weighttype get1DWeight(T x)
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_)
virtual weighttype get2DWeight(T x, T y)
virtual weighttype getWeight(CoordinatePoint< T > p)
virtual weighttype get3DWeight(T x, T y, T z)
virtual weighttype getWeight(CoordinatePoint< T > P)=0
virtual weighttype get1DWeight(T x)=0
virtual weighttype get2DWeight(T x, T y)=0
virtual weighttype get3DWeight(T x, T y, T z)=0
PartitioningProblem sets up partitioning problems for the user.
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
void solve(bool updateInputData=true)
Direct the problem to create a solution.
map_t::local_ordinal_type lno_t
map_t::global_ordinal_type gno_t
const std::string shapes[]
void getCoords(void *data, int numGid, int numLid, int numObj, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids, int dim, double *coords_, 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)
int getDim(void *data, int *ierr)
int getNumObj(void *data, int *ierr)
const std::string weight_distribution_string
static ArrayRCP< ArrayRCP< zscalar_t > > weights
Tpetra::MultiVector< zscalar_t, zlno_t, zgno_t, znode_t > tMVector_t