Zoltan 2 Version 0.5
GeometricGenerator.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //   Zoltan2: A package of combinatorial algorithms for scientific computing
00006 //                  Copyright 2012 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Karen Devine      (kddevin@sandia.gov)
00039 //                    Erik Boman        (egboman@sandia.gov)
00040 //                    Siva Rajamanickam (srajama@sandia.gov)
00041 //
00042 // ***********************************************************************
00043 //
00044 // @HEADER
00045 #include <Teuchos_Comm.hpp>
00046 #include <Teuchos_ParameterList.hpp>
00047 #include <Teuchos_FilteredIterator.hpp>
00048 #include <Teuchos_ParameterEntry.hpp>
00049 #include <iostream>
00050 #include <ctime>
00051 #include <limits>
00052 #include <climits>
00053 #include <string>
00054 #include <cstdlib>
00055 #include <sstream>
00056 #include <fstream>
00057 #include <Tpetra_MultiVector_decl.hpp>
00058 #include <Zoltan2_XpetraMultiVectorInput.hpp>
00059 #include <Zoltan2_PartitioningSolution.hpp>
00060 #include <Teuchos_ArrayViewDecl.hpp>
00061 
00062 enum shape {SQUARE, RECTANGLE, CIRCLE, CUBE, RECTANGULAR_PRISM, SPHERE};
00063 const std::string shapes[] = {"SQUARE", "RECTANGLE", "CIRCLE", "CUBE", "RECTANGULAR_PRISM", "SPHERE"};
00064 #define SHAPE_COUNT 6
00065 
00066 enum distribution {normal, uniform};
00067 const std::string distribution[] = {"distribution", "uniform"};
00068 #define DISTRIBUTION_COUNT 2
00069 
00070 #define HOLE_ALLOC_STEP  10
00071 #define MAX_WEIGHT_DIM  10
00072 #define INVALID(STR) "Invalid argument at " + STR
00073 #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"
00074 
00075 #define INVALID_SHAPE_ARG(SHAPE, REQUIRED) "Invalid argument count for shape " + SHAPE + ". Requires " + REQUIRED + " argument(s)."
00076 #define MAX_ITER_ALLOWED 500
00077 
00078 const std::string weight_distribution = "WeightDistribution";
00079 const std::string weight_distribution_string = weight_distribution + "-";
00080 template <typename T>
00081 struct CoordinatePoint {
00082   T x;
00083   T y;
00084   T z;
00085   /*
00086   bool isInArea(int dim, T *minCoords, T *maxCoords){
00087     dim = min(3, dim);
00088     for (int i = 0; i < dim; ++i){
00089       bool maybe = true;
00090       switch(i){
00091       case 0:
00092         maybe = x < maxCoords[i] && x > maxCoords[i];
00093         break;
00094       case 1:
00095         maybe = y < maxCoords[i] && y > maxCoords[i];
00096         break;
00097       case 2:
00098         maybe = z < maxCoords[i] && z > maxCoords[i];
00099         break;
00100       }
00101       if (!maybe){
00102         return false;
00103       }
00104     }
00105     return true;
00106   }
00107    */
00108   CoordinatePoint(){
00109     x = 0;y=0;z=0;
00110   }
00111 };
00112 
00113 template <typename T>
00114 class Hole{
00115 public:
00116   CoordinatePoint<T> center;
00117   T edge1, edge2, edge3;
00118   Hole(CoordinatePoint<T> center_, T edge1_, T edge2_, T edge3_){
00119     this->center.x = center_.x;
00120     this->center.y = center_.y;
00121     this->center.z = center_.z;
00122     this->edge1 = edge1_;
00123     this->edge2 = edge2_;
00124     this->edge3 = edge3_;
00125   }
00126   virtual bool isInArea(CoordinatePoint<T> dot) = 0;
00127 };
00128 
00129 template <typename T>
00130 class SquareHole:public Hole<T>{
00131 public:
00132   SquareHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00133 
00134   virtual bool isInArea(CoordinatePoint<T> dot){
00135     return fabs(dot.x - this->center.x) < this->edge1 / 2 && fabs(dot.y - this->center.y) < this->edge1 / 2;
00136   }
00137 };
00138 
00139 template <typename T>
00140 class RectangleHole:public Hole<T>{
00141 public:
00142   RectangleHole(CoordinatePoint<T> center_  , T edge_x_, T edge_y_): Hole<T>(center_, edge_x_,  edge_y_, 0){}
00143   virtual bool isInArea(CoordinatePoint<T> dot){
00144     return fabs(dot.x - this->center.x) < this->edge1 / 2 && fabs(dot.y - this->center.y) < this->edge2 / 2;
00145   }
00146 };
00147 
00148 template <typename T>
00149 class CircleHole:public Hole<T>{
00150 public:
00151   CircleHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00152   virtual bool isInArea(CoordinatePoint<T> dot){
00153     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;
00154   }
00155 };
00156 
00157 template <typename T>
00158 class CubeHole:public Hole<T>{
00159 public:
00160   CubeHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00161   virtual bool isInArea(CoordinatePoint<T> dot){
00162     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;
00163   }
00164 };
00165 
00166 template <typename T>
00167 class RectangularPrismHole:public Hole<T>{
00168 public:
00169   RectangularPrismHole(CoordinatePoint<T> center_  , T edge_x_, T edge_y_, T edge_z_): Hole<T>(center_, edge_x_,  edge_y_, edge_z_){}
00170   virtual bool isInArea(CoordinatePoint<T> dot){
00171     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;
00172   }
00173 };
00174 
00175 template <typename T>
00176 class SphereHole:public Hole<T>{
00177 public:
00178   SphereHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00179   virtual bool isInArea(CoordinatePoint<T> dot){
00180     return  fabs((dot.x - this->center.x) * (dot.x - this->center.x) * (dot.x - this->center.x)) +
00181         fabs((dot.y - this->center.y) * (dot.y - this->center.y) * (dot.y - this->center.y)) +
00182         fabs((dot.z - this->center.z) * (dot.z - this->center.z) * (dot.z - this->center.z))
00183         <
00184         this->edge1 * this->edge1 * this->edge1;
00185   }
00186 };
00187 
00188 template <typename T, typename weighttype>
00189 class WeightDistribution{
00190 public:
00191   virtual weighttype getWeight(CoordinatePoint<T> P)=0;
00192   virtual weighttype get1DWeight(T x)=0;
00193   virtual weighttype get2DWeight(T x, T y)=0;
00194   virtual weighttype get3DWeight(T x, T y, T z)=0;
00195   WeightDistribution(){};
00196   virtual ~WeightDistribution(){};
00197 };
00198 
00217 template <typename T, typename weighttype>
00218 class SteppedEquation:public WeightDistribution<T, weighttype>{
00219   T a1,a2,a3;
00220   T b1,b2,b3;
00221   T c;
00222   T x1,y1,z1;
00223 
00224   T *steps;
00225   weighttype *values;
00226   int stepCount;
00227 public:
00228   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>(){
00229     this->a1 = a1_;
00230     this->a2 = a2_;
00231     this->a3 = a3_;
00232     this->b1 = b1_;
00233     this->b2 = b2_;
00234     this->b3 = b3_;
00235     this->c = c_;
00236     this->x1 = x1_;
00237     this->y1 = y1_;
00238     this->z1 = z1_;
00239 
00240 
00241     this->stepCount = stepCount_;
00242     if(this->stepCount > 0){
00243       this->steps = new T[this->stepCount];
00244       this->values = new T[this->stepCount + 1];
00245 
00246       for (int i = 0; i < this->stepCount; ++i){
00247         this->steps[i] = steps_[i];
00248         this->values[i] = values_[i];
00249       }
00250       this->values[this->stepCount] = values_[this->stepCount];
00251     }
00252 
00253   }
00254 
00255   virtual ~SteppedEquation(){
00256     if(this->stepCount > 0){
00257       delete [] this->steps;
00258       delete [] this->values;
00259     }
00260   }
00261 
00262 
00263   virtual weighttype get1DWeight(T x){
00264     T expressionRes = this->a1 * pow( (x - this->x1), b1) + c;
00265     if(this->stepCount > 0){
00266       for (int i = 0; i < this->stepCount; ++i){
00267         if (expressionRes < this->steps[i]) return this->values[i];
00268       }
00269       return this->values[this->stepCount];
00270     }
00271     else {
00272       return weighttype(expressionRes);
00273     }
00274   }
00275 
00276   virtual weighttype get2DWeight(T x, T y){
00277     T expressionRes = this->a1 * pow( (x - this->x1), b1) + this->a2 * pow( (y - this->y1), b2) + c;
00278     if(this->stepCount > 0){
00279       for (int i = 0; i < this->stepCount; ++i){
00280         if (expressionRes < this->steps[i]) return this->values[i];
00281       }
00282       return this->values[this->stepCount];
00283     }
00284     else {
00285       return weighttype(expressionRes);
00286     }
00287   }
00288 
00289   void print (T x, T y, T z){
00290     cout << this->a1 << "*" <<  "math.pow( (" << x  << "- " <<  this->x1 << " ), " << b1 <<")" << "+" <<  this->a2<< "*"<<  "math.pow( (" << y << "-" <<  this->y1 << "), " <<
00291         b2 << " ) + " << this->a3 << " * math.pow( (" << z << "-" <<  this->z1 << "), " << b3 << ")+ " << c << " == " <<
00292         this->a1 * pow( (x - this->x1), b1) + this->a2 * pow( (y - this->y1), b2) + this->a3 * pow( (z - this->z1), b3) + c << endl;
00293 
00294   }
00295 
00296   virtual weighttype get3DWeight(T x, T y, T z){
00297     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;
00298 
00299     //this->print(x,y,z);
00300     if(this->stepCount > 0){
00301       for (int i = 0; i < this->stepCount; ++i){
00302         if (expressionRes < this->steps[i]) {
00303           //cout << "0exp:" << expressionRes << " step:" << steps[i] << " value:" << values[i] << endl;
00304           return this->values[i];
00305         }
00306       }
00307 
00308       //cout << "1exp:" << expressionRes << " step:" << steps[stepCount] << " value:" << values[stepCount] << endl;
00309       return this->values[this->stepCount];
00310     }
00311     else {
00312       return weighttype(expressionRes);
00313     }
00314   }
00315 
00316   virtual weighttype getWeight(CoordinatePoint<T> p){
00317     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);
00318     if(this->stepCount > 0){
00319       for (int i = 0; i < this->stepCount; ++i){
00320         if (expressionRes < this->steps[i]) return this->values[i];
00321       }
00322       return this->values[this->stepCount];
00323     }
00324     else {
00325       return weighttype(expressionRes);
00326     }
00327   }
00328 };
00329 
00330 template <typename T, typename lno_t, typename gno_t>
00331 class CoordinateDistribution{
00332 public:
00333   gno_t numPoints;
00334   int dimension;
00335   lno_t requested;
00336   gno_t assignedPrevious;
00337   int worldSize;
00338   virtual ~CoordinateDistribution(){}
00339 
00340   CoordinateDistribution(gno_t np_, int dim, int worldSize):numPoints(np_), dimension(dim), requested(0), assignedPrevious(0), worldSize(worldSize){}
00341   virtual CoordinatePoint<T> getPoint(gno_t point_index, unsigned int &state) = 0;
00342   virtual T getXCenter() = 0;
00343   virtual T getXRadius() =0;
00344 
00345   void GetPoints(lno_t requestedPointcount, CoordinatePoint<T> *points /*preallocated sized numPoints*/,
00346       Hole<T> **holes, lno_t holeCount,
00347       float *sharedRatios_, int myRank){
00348 
00349     for (int i = 0; i < myRank; ++i){
00350       //cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< endl;
00351       this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
00352       if (i < this->numPoints % this->worldSize){
00353         this->assignedPrevious += 1;
00354       }
00355     }
00356 
00357     this->requested = requestedPointcount;
00358 
00359     unsigned int slice =  UINT_MAX/(this->worldSize);
00360     unsigned int stateBegin = myRank * slice;
00361 
00362 //#ifdef HAVE_ZOLTAN2_OMP
00363 //#pragma omp parallel for
00364 //#endif
00365 #ifdef HAVE_ZOLTAN2_OMP
00366 #pragma omp parallel
00367 #endif
00368     {
00369       int me = 0;
00370       int tsize = 1;
00371 #ifdef HAVE_ZOLTAN2_OMP
00372       me = omp_get_thread_num();
00373       tsize = omp_get_num_threads();
00374 #endif
00375       unsigned int state = stateBegin + me * slice/(tsize);
00376 
00377 #ifdef HAVE_ZOLTAN2_OMP
00378 #pragma omp for
00379 #endif
00380       for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
00381         lno_t iteration = 0;
00382         while(1){
00383           if(++iteration > MAX_ITER_ALLOWED) {
00384             throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00385           }
00386           CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, &state);
00387 
00388           bool isInHole = false;
00389           for(lno_t i = 0; i < holeCount; ++i){
00390             if(holes[i][0].isInArea(p)){
00391               isInHole = true;
00392               break;
00393             }
00394           }
00395           if(isInHole) continue;
00396           points[cnt].x = p.x;
00397 
00398           points[cnt].y = p.y;
00399           points[cnt].z = p.z;
00400           break;
00401         }
00402       }
00403     }
00404 //#pragma omp parallel
00405       /*
00406     {
00407 
00408       lno_t cnt = 0;
00409       lno_t iteration = 0;
00410       while (cnt < requestedPointcount){
00411         if(++iteration > MAX_ITER_ALLOWED) {
00412           throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00413         }
00414         CoordinatePoint <T> p = this->getPoint();
00415 
00416         bool isInHole = false;
00417         for(lno_t i = 0; i < holeCount; ++i){
00418           if(holes[i][0].isInArea(p)){
00419             isInHole = true;
00420             break;
00421           }
00422         }
00423         if(isInHole) continue;
00424         iteration = 0;
00425         points[cnt].x = p.x;
00426         points[cnt].y = p.y;
00427         points[cnt].z = p.z;
00428         ++cnt;
00429       }
00430     }
00431     */
00432   }
00433 
00434   void GetPoints(lno_t requestedPointcount, T **coords/*preallocated sized numPoints*/, lno_t tindex,
00435       Hole<T> **holes, lno_t holeCount,
00436       float *sharedRatios_, int myRank){
00437 
00438     for (int i = 0; i < myRank; ++i){
00439       //cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< endl;
00440       this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
00441       if (i < this->numPoints % this->worldSize){
00442         this->assignedPrevious += 1;
00443       }
00444     }
00445 
00446     this->requested = requestedPointcount;
00447 
00448     unsigned int slice =  UINT_MAX/(this->worldSize);
00449     unsigned int stateBegin = myRank * slice;
00450 #ifdef HAVE_ZOLTAN2_OMP
00451 #pragma omp parallel
00452 #endif
00453     {
00454       int me = 0;
00455       int tsize = 1;
00456 #ifdef HAVE_ZOLTAN2_OMP
00457       me = omp_get_thread_num();
00458       tsize = omp_get_num_threads();
00459 #endif
00460       unsigned int state = stateBegin + me * (slice/(tsize));
00461       /*
00462 #pragma omp critical
00463       {
00464 
00465         cout << "myRank:" << me << " stateBeg:" << stateBegin << " tsize:" << tsize << " state:" << state <<  " slice: " << slice / tsize <<  endl;
00466       }
00467       */
00468 #ifdef HAVE_ZOLTAN2_OMP
00469 #pragma omp for
00470 #endif
00471       for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
00472         lno_t iteration = 0;
00473         while(1){
00474           if(++iteration > MAX_ITER_ALLOWED) {
00475             throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00476           }
00477           CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, state);
00478 
00479           bool isInHole = false;
00480           for(lno_t i = 0; i < holeCount; ++i){
00481             if(holes[i][0].isInArea(p)){
00482               isInHole = true;
00483               break;
00484             }
00485           }
00486           if(isInHole) continue;
00487           coords[0][cnt + tindex] = p.x;
00488           if(this->dimension > 1){
00489             coords[1][cnt + tindex] = p.y;
00490             if(this->dimension > 2){
00491               coords[2][cnt + tindex] = p.z;
00492             }
00493           }
00494           break;
00495         }
00496       }
00497     }
00498   }
00499 };
00500 
00501 template <typename T, typename lno_t, typename gno_t>
00502 class CoordinateNormalDistribution:public CoordinateDistribution<T,lno_t,gno_t>{
00503 public:
00504   CoordinatePoint<T> center;
00505   T standartDevx;
00506   T standartDevy;
00507   T standartDevz;
00508 
00509 
00510   virtual T getXCenter(){
00511     return center.x;
00512   }
00513   virtual T getXRadius(){
00514     return standartDevx;
00515   }
00516 
00517   CoordinateNormalDistribution(gno_t np_, int dim, CoordinatePoint<T> center_ , T sd_x, T sd_y, T sd_z, int worldSize): CoordinateDistribution<T,lno_t,gno_t>(np_,dim,worldSize),
00518       standartDevx(sd_x), standartDevy(sd_y), standartDevz(sd_z){
00519     this->center.x = center_.x;
00520     this->center.y = center_.y;
00521     this->center.z = center_.z;
00522   }
00523 
00524   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00525 
00526     pindex = 0; // not used in normal distribution.
00527     CoordinatePoint <T> p;
00528 
00529     for(int i = 0; i < this->dimension; ++i){
00530       switch(i){
00531       case 0:
00532         p.x = normalDist(this->center.x, this->standartDevx, state);
00533         break;
00534       case 1:
00535         p.y = normalDist(this->center.y, this->standartDevy, state);
00536         break;
00537       case 2:
00538         p.z = normalDist(this->center.z, this->standartDevz, state);
00539         break;
00540       default:
00541         throw "unsupported dimension";
00542       }
00543     }
00544     return p;
00545   }
00546 
00547   virtual ~CoordinateNormalDistribution(){};
00548 private:
00549   T normalDist(T center, T sd, unsigned int &state) {
00550     static bool derived=false;
00551     static T storedDerivation;
00552     T polarsqrt, normalsquared, normal1, normal2;
00553     if (!derived) {
00554       do {
00555         normal1=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
00556         normal2=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
00557         normalsquared=normal1*normal1+normal2*normal2;
00558       } while ( normalsquared>=1.0 || normalsquared == 0.0);
00559 
00560       polarsqrt=sqrt(-2.0*log(normalsquared)/normalsquared);
00561       storedDerivation=normal1*polarsqrt;
00562       derived=true;
00563       return normal2*polarsqrt*sd + center;
00564     }
00565     else {
00566       derived=false;
00567       return storedDerivation*sd + center;
00568     }
00569   }
00570 };
00571 
00572 template <typename T, typename lno_t, typename gno_t>
00573 class CoordinateUniformDistribution:public CoordinateDistribution<T,lno_t, gno_t>{
00574 public:
00575   T leftMostx;
00576   T rightMostx;
00577   T leftMosty;
00578   T rightMosty;
00579   T leftMostz;
00580   T rightMostz;
00581 
00582 
00583   virtual T getXCenter(){
00584     return (rightMostx - leftMostx)/2  + leftMostx;
00585   }
00586   virtual T getXRadius(){
00587     return (rightMostx - leftMostx)/2;
00588   }
00589 
00590 
00591   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 worldSize ): CoordinateDistribution<T,lno_t,gno_t>(np_,dim,worldSize),
00592       leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y), leftMostz(l_z), rightMostz(r_z){}
00593 
00594   virtual ~CoordinateUniformDistribution(){};
00595   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00596 
00597 
00598     pindex = 0; //not used in uniform dist.
00599     CoordinatePoint <T> p;
00600     for(int i = 0; i < this->dimension; ++i){
00601       switch(i){
00602       case 0:
00603         p.x = uniformDist(this->leftMostx, this->rightMostx, state);
00604         break;
00605       case 1:
00606         p.y = uniformDist(this->leftMosty, this->rightMosty, state);
00607         break;
00608       case 2:
00609         p.z = uniformDist(this->leftMostz, this->rightMostz, state);
00610         break;
00611       default:
00612         throw "unsupported dimension";
00613       }
00614     }
00615     return p;
00616   }
00617 
00618 private:
00619 
00620   T uniformDist(T a, T b, unsigned int &state)
00621   {
00622     return (b-a)*(rand_r(&state) / double(RAND_MAX)) + a;
00623   }
00624 };
00625 
00626 template <typename T, typename lno_t, typename gno_t>
00627 class CoordinateGridDistribution:public CoordinateDistribution<T,lno_t,gno_t>{
00628 public:
00629   T leftMostx;
00630   T rightMostx;
00631   T leftMosty;
00632   T rightMosty;
00633   T leftMostz;
00634   T rightMostz;
00635   gno_t along_X, along_Y, along_Z;
00636   //T currentX, currentY, currentZ;
00637   T processCnt;
00638   int myRank;
00639   T xstep, ystep, zstep;
00640   gno_t xshift, yshift, zshift;
00641 
00642   virtual T getXCenter(){
00643     return (rightMostx - leftMostx)/2  + leftMostx;
00644   }
00645   virtual T getXRadius(){
00646     return (rightMostx - leftMostx)/2;
00647   }
00648 
00649 
00650   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 worldSize): CoordinateDistribution<T,lno_t,gno_t>(alongX * alongY * alongZ,dim,worldSize),
00651       leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y), leftMostz(l_z), rightMostz(r_z), myRank(myRank_){
00652     //currentX = leftMostx, currentY = leftMosty, currentZ = leftMostz;
00653     this->processCnt = 0;
00654     this->along_X = alongX; this->along_Y = alongY; this->along_Z = alongZ;
00655 
00656     if(this->along_X > 1)
00657       this->xstep = (rightMostx - leftMostx) / (alongX - 1);
00658     else
00659       this->xstep = 1;
00660     if(this->along_Y > 1)
00661       this->ystep = (rightMosty - leftMosty) / (alongY - 1);
00662     else
00663       this->ystep = 1;
00664     if(this->along_Z > 1)
00665       this->zstep = (rightMostz - leftMostz) / (alongZ - 1);
00666     else
00667       this->zstep = 1;
00668     xshift = 0; yshift = 0; zshift = 0;
00669   }
00670 
00671   virtual ~CoordinateGridDistribution(){};
00672   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00673     //lno_t before = processCnt + this->assignedPrevious;
00674     //cout << "before:" << processCnt << " " << this->assignedPrevious << endl;
00675     //lno_t xshift = 0, yshift = 0, zshift = 0;
00676 
00677     //lno_t tmp = before % (this->along_X * this->along_Y);
00678     //xshift  = tmp % this->along_X;
00679     //yshift = tmp / this->along_X;
00680     //zshift = before / (this->along_X * this->along_Y);
00681 
00682     state = 0; //not used here
00683     this->zshift = pindex / (along_X * along_Y);
00684     this->yshift = (pindex % (along_X * along_Y)) / along_X;
00685     this->xshift = (pindex % (along_X * along_Y)) % along_X;
00686 
00687     //this->xshift = pindex / (along_Z * along_Y);
00688    // this->zshift = (pindex % (along_Z * along_Y)) / along_Y;
00689     //this->yshift = (pindex % (along_Z * along_Y)) % along_Y;
00690 
00691     CoordinatePoint <T> p;
00692     p.x = xshift * this->xstep + leftMostx;
00693     p.y = yshift * this->ystep + leftMosty;
00694     p.z = zshift * this->zstep + leftMostz;
00695 /*
00696     ++xshift;
00697     if(xshift == this->along_X){
00698       ++yshift;
00699       xshift = 0;
00700       if(yshift == this->along_Y){
00701         ++zshift;
00702         yshift = 0;
00703       }
00704     }
00705 */
00706     /*
00707     if(this->processCnt == 0){
00708       this->xshift = this->assignedPrevious / (along_Z * along_Y);
00709       //this->yshift = (this->assignedPrevious % (along_X * along_Y)) / along_X;
00710       this->zshift = (this->assignedPrevious % (along_Z * along_Y)) / along_Y;
00711       //this->xshift = (this->assignedPrevious % (along_X * along_Y)) % along_X;
00712       this->yshift = (this->assignedPrevious % (along_Z * along_Y)) % along_Y;
00713       ++this->processCnt;
00714     }
00715 
00716     CoordinatePoint <T> p;
00717     p.x = xshift * this->xstep + leftMostx;
00718     p.y = yshift * this->ystep + leftMosty;
00719     p.z = zshift * this->zstep + leftMostz;
00720 
00721     ++yshift;
00722     if(yshift == this->along_Y){
00723       ++zshift;
00724       yshift = 0;
00725       if(zshift == this->along_Z){
00726         ++xshift;
00727         zshift = 0;
00728       }
00729 
00730     }
00731     */
00732     /*
00733     if(this->requested - 1 > this->processCnt){
00734       this->processCnt++;
00735     } else {
00736       cout << "req:" << this->requested << " pr:" <<this->processCnt << endl;
00737       cout << "p:" << p.x << " " << p.y << " " << p.z <<  endl ;
00738       cout << "s:" << xshift << " " << yshift << " " << zshift <<  endl ;
00739       cout << "st:" << this->xstep << " " << this->ystep << " " << this->zstep <<  endl ;
00740     }
00741     */
00742     return p;
00743   }
00744 
00745 private:
00746 
00747 };
00748 
00749 template <typename T, typename lno_t, typename gno_t, typename node_t>
00750 class GeometricGenerator {
00751 private:
00752   Hole<T> **holes; //to represent if there is any hole in the input
00753   int holeCount;
00754   int coordinate_dimension;  //dimension of the geometry
00755   gno_t numGlobalCoords;  //global number of coordinates requested to be created.
00756   lno_t numLocalCoords;
00757   float *loadDistributions; //sized as the number of processors, the load of each processor.
00758   bool loadDistSet;
00759   bool distinctCoordSet;
00760   CoordinateDistribution<T, lno_t,gno_t> **coordinateDistributions;
00761   int distributionCount;
00762   //CoordinatePoint<T> *points;
00763   T **coords;
00764   T **wghts;
00765   WeightDistribution<T,T> **wd;
00766   int weight_dimension;  //dimension of the geometry
00767 
00768   //RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmVector;
00769   //RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmwVector;
00770   int worldSize;
00771   int myRank;
00772   T minx;
00773   T maxx;
00774   T miny;
00775   T maxy;
00776   T minz;
00777   T maxz;
00778   std::string outfile;
00779 
00780 
00781 
00782   template <typename tt>
00783   tt getParamVal( const Teuchos::ParameterEntry& pe, const std::string &paramname){
00784     tt returnVal;
00785     try {
00786       returnVal = Teuchos::getValue<tt>(pe);
00787     }
00788     catch (...){
00789       throw INVALID(paramname);
00790     }
00791     return returnVal;
00792   }
00793 
00794   int countChar (std::string inStr, char countChar){
00795     int cnt = 0;
00796     for (unsigned int i = 0; i < inStr.size(); ++i){
00797       if (inStr[i] == countChar) {
00798         cnt++;
00799       }
00800     }
00801     return cnt;
00802   }
00803 
00804   template <typename tt>
00805   std::string toString(tt obj){
00806     std::stringstream ss (std::stringstream::in |std::stringstream::out);
00807     ss << obj;
00808     std::string tmp = "";
00809     ss >> tmp;
00810     return tmp;
00811   }
00812 
00813   template <typename tt>
00814   tt fromString(std::string obj){
00815     std::stringstream ss (std::stringstream::in | std::stringstream::out);
00816     ss << obj;
00817     tt tmp;
00818     ss >> tmp;
00819     if (ss.fail()){
00820       throw "Cannot convert string " + obj;
00821     }
00822     return tmp;
00823   }
00824 
00825 
00826   void splitString(std::string inStr, char splitChar, std::string *outSplittedStr){
00827     std::stringstream ss (std::stringstream::in | std::stringstream::out);
00828     ss << inStr;
00829     int cnt = 0;
00830     while (!ss.eof()){
00831       std::string tmp = "";
00832       std::getline(ss, tmp ,splitChar);
00833       outSplittedStr[cnt++] = tmp;
00834     }
00835   }
00836 
00837   /*
00838   void getDistinctCoordinateDescription(std::string distinctDescription){
00839 
00840     this->distinctCoordinates = new bool[this->dimension];
00841     if(distinctDescription != ""){
00842       int argCnt = this->countChar(distinctDescription, ',') + 1;
00843       if(argCnt != this->dimension) {
00844         throw "Invalid parameter count for distinct_coordinates. Size should be equal to dimension.";
00845       }
00846 
00847       std::string *splittedStr = new std::string[argCnt];
00848       splitString(holeDescription, ',', splittedStr);
00849 
00850       for(int i = 0; i < argCnt; ++i){
00851         if(splittedStr[i] == "T"){
00852           distinctCoordinates[i] = true;
00853         } else if(splittedStr[i] == "F"){
00854           distinctCoordinates[i] = false;
00855         } else {
00856           throw "Invalid parameter " + splittedStr[i] + " for distinct_coordinates.";
00857         }
00858       }
00859       delete []splittedStr;
00860     }
00861     else {
00862       for(int i = 0; i < this->dimension; ++i){
00863         distinctCoordinates[i] = false;
00864       }
00865     }
00866   }
00867    */
00868 
00869 
00870   void getCoordinateDistributions(std::string coordinate_distributions){
00871     if(coordinate_distributions == ""){
00872       throw "At least one distribution function must be provided to coordinate generator.";
00873     }
00874 
00875     int argCnt = this->countChar(coordinate_distributions, ',') + 1;
00876     std::string *splittedStr = new std::string[argCnt];
00877     splitString(coordinate_distributions, ',', splittedStr);
00878     coordinateDistributions = (CoordinateDistribution<T, lno_t,gno_t> **) malloc(sizeof (CoordinateDistribution<T, lno_t,gno_t> *) * 1);
00879     for(int i = 0; i < argCnt; ){
00880       coordinateDistributions = (CoordinateDistribution<T, lno_t,gno_t> **)realloc((void *)coordinateDistributions, (this->distributionCount + 1)* sizeof(CoordinateDistribution<T, lno_t,gno_t> *));
00881 
00882       std::string distName = splittedStr[i++];
00883       gno_t np_ = 0;
00884       if(distName == "NORMAL"){
00885         int reqArg = 5;
00886         if (this->coordinate_dimension == 3){
00887           reqArg = 7;
00888         }
00889         if(i + reqArg > argCnt) {
00890           std::string tmp = toString<int>(reqArg);
00891           throw INVALID_SHAPE_ARG(distName, tmp);
00892         }
00893         np_ = fromString<gno_t>(splittedStr[i++]);
00894         CoordinatePoint<T> pp;
00895 
00896         pp.x = fromString<T>(splittedStr[i++]);
00897         pp.y = fromString<T>(splittedStr[i++]);
00898         pp.z = 0;
00899         if(this->coordinate_dimension == 3){
00900           pp.z = fromString<T>(splittedStr[i++]);
00901         }
00902 
00903         T sd_x = fromString<T>(splittedStr[i++]);
00904         T sd_y = fromString<T>(splittedStr[i++]);
00905         T sd_z = 0;
00906         if(this->coordinate_dimension == 3){
00907           sd_z = fromString<T>(splittedStr[i++]);
00908         }
00909         this->coordinateDistributions[this->distributionCount++] = new CoordinateNormalDistribution<T, lno_t,gno_t>(np_, this->coordinate_dimension, pp , sd_x, sd_y, sd_z, worldSize );
00910 
00911       } else if(distName == "UNIFORM" ){
00912         int reqArg = 5;
00913         if (this->coordinate_dimension == 3){
00914           reqArg = 7;
00915         }
00916         if(i + reqArg > argCnt) {
00917           std::string tmp = toString<int>(reqArg);
00918           throw INVALID_SHAPE_ARG(distName, tmp);
00919         }
00920         np_ = fromString<gno_t>(splittedStr[i++]);
00921         T l_x = fromString<T>(splittedStr[i++]);
00922         T r_x = fromString<T>(splittedStr[i++]);
00923         T l_y = fromString<T>(splittedStr[i++]);
00924         T r_y = fromString<T>(splittedStr[i++]);
00925 
00926         T l_z = 0, r_z = 0;
00927 
00928         if(this->coordinate_dimension == 3){
00929           l_z = fromString<T>(splittedStr[i++]);
00930           r_z = fromString<T>(splittedStr[i++]);
00931         }
00932 
00933         this->coordinateDistributions[this->distributionCount++] = new CoordinateUniformDistribution<T, lno_t,gno_t>( np_,  this->coordinate_dimension, l_x, r_x, l_y, r_y, l_z, r_z, worldSize );
00934       } else if (distName == "GRID"){
00935         int reqArg = 6;
00936         if(this->coordinate_dimension == 3){
00937           reqArg = 9;
00938         }
00939         if(i + reqArg > argCnt) {
00940           std::string tmp = toString<int>(reqArg);
00941           throw INVALID_SHAPE_ARG(distName, tmp);
00942         }
00943 
00944         gno_t np_x = fromString<gno_t>(splittedStr[i++]);
00945         gno_t np_y = fromString<gno_t>(splittedStr[i++]);
00946         gno_t np_z = 1;
00947 
00948 
00949         if(this->coordinate_dimension == 3){
00950           np_z = fromString<gno_t>(splittedStr[i++]);
00951         }
00952 
00953         np_ = np_x * np_y * np_z;
00954         T l_x = fromString<T>(splittedStr[i++]);
00955         T r_x = fromString<T>(splittedStr[i++]);
00956         T l_y = fromString<T>(splittedStr[i++]);
00957         T r_y = fromString<T>(splittedStr[i++]);
00958 
00959         T l_z = 0, r_z = 0;
00960 
00961         if(this->coordinate_dimension == 3){
00962           l_z = fromString<T>(splittedStr[i++]);
00963           r_z = fromString<T>(splittedStr[i++]);
00964         }
00965 
00966         if(np_x < 1 || np_z < 1 || np_y < 1 ){
00967           throw "Provide at least 1 point along each dimension for grid test.";
00968         }
00969         //cout << "ly:" << l_y << " ry:" << r_y << endl;
00970         this->coordinateDistributions[this->distributionCount++] = new CoordinateGridDistribution<T, lno_t,gno_t>
00971         (np_x, np_y,np_z, this->coordinate_dimension, l_x, r_x,l_y, r_y, l_z, r_z , this->myRank, worldSize);
00972 
00973       }
00974       else {
00975         std::string tmp = toString<int>(this->coordinate_dimension);
00976         throw INVALIDSHAPE(distName, tmp);
00977       }
00978       this->numGlobalCoords += (gno_t) np_;
00979     }
00980     delete []splittedStr;
00981   }
00982 
00983   void getProcLoadDistributions(std::string proc_load_distributions){
00984 
00985     this->loadDistributions = new float[this->worldSize];
00986     if(proc_load_distributions == ""){
00987       float r = 1.0 / this->worldSize;
00988       for(int i = 0; i < this->worldSize; ++i){
00989         this->loadDistributions[i] = r;
00990       }
00991     } else{
00992 
00993 
00994       int argCnt = this->countChar(proc_load_distributions, ',') + 1;
00995       if(argCnt != this->worldSize) {
00996         throw "Invalid parameter count load distributions. Given " + toString<int>(argCnt) + " processor size is " + toString<int>(worldSize);
00997       }
00998       std::string *splittedStr = new std::string[argCnt];
00999       splitString(proc_load_distributions, ',', splittedStr);
01000       for(int i = 0; i < argCnt; ++i){
01001         this->loadDistributions[i] = fromString<float>(splittedStr[i]);
01002       }
01003       delete []splittedStr;
01004 
01005 
01006       float sum = 0;
01007       for(int i = 0; i < this->worldSize; ++i){
01008         sum += this->loadDistributions[i];
01009       }
01010       if (fabs(sum - 1.0) > 10*std::numeric_limits<float>::epsilon()){
01011         throw "Processor load ratios do not sum to 1.0.";
01012       }
01013     }
01014 
01015   }
01016 
01017   void getHoles(std::string holeDescription){
01018     if(holeDescription == ""){
01019       return;
01020     }
01021     this->holes =  (Hole<T> **) malloc(sizeof (Hole <T>*));
01022     int argCnt = this->countChar(holeDescription, ',') + 1;
01023     std::string *splittedStr = new std::string[argCnt];
01024     splitString(holeDescription, ',', splittedStr);
01025 
01026     for(int i = 0; i < argCnt; ){
01027       holes = (Hole<T> **)realloc((void *)holes, (this->holeCount + 1)* sizeof(Hole<T> *));
01028 
01029       std::string shapeName = splittedStr[i++];
01030       if(shapeName == "SQUARE" && this->coordinate_dimension == 2){
01031         if(i + 3 > argCnt) {
01032           throw INVALID_SHAPE_ARG(shapeName, "3");
01033         }
01034         CoordinatePoint<T> pp;
01035         pp.x = fromString<T>(splittedStr[i++]);
01036         pp.y = fromString<T>(splittedStr[i++]);
01037         T edge = fromString<T>(splittedStr[i++]);
01038         this->holes[this->holeCount++] = new SquareHole<T>(pp, edge);
01039       } else if(shapeName == "RECTANGLE" && this->coordinate_dimension == 2){
01040         if(i + 4 > argCnt) {
01041           throw INVALID_SHAPE_ARG(shapeName, "4");
01042         }
01043         CoordinatePoint<T> pp;
01044         pp.x = fromString<T>(splittedStr[i++]);
01045         pp.y = fromString<T>(splittedStr[i++]);
01046         T edgex = fromString<T>(splittedStr[i++]);
01047         T edgey = fromString<T>(splittedStr[i++]);
01048 
01049         this->holes[this->holeCount++] = new RectangleHole<T>(pp, edgex,edgey);
01050       } else if(shapeName == "CIRCLE" && this->coordinate_dimension == 2){
01051         if(i + 3 > argCnt) {
01052           throw INVALID_SHAPE_ARG(shapeName, "3");
01053         }
01054         CoordinatePoint<T> pp;
01055         pp.x = fromString<T>(splittedStr[i++]);
01056         pp.y = fromString<T>(splittedStr[i++]);
01057         T r = fromString<T>(splittedStr[i++]);
01058         this->holes[this->holeCount++] = new CircleHole<T>(pp, r);
01059       }  else if(shapeName == "CUBE" && this->coordinate_dimension == 3){
01060         if(i + 4 > argCnt) {
01061           throw INVALID_SHAPE_ARG(shapeName, "4");
01062         }
01063         CoordinatePoint<T> pp;
01064         pp.x = fromString<T>(splittedStr[i++]);
01065         pp.y = fromString<T>(splittedStr[i++]);
01066         pp.z = fromString<T>(splittedStr[i++]);
01067         T edge = fromString<T>(splittedStr[i++]);
01068         this->holes[this->holeCount++] = new CubeHole<T>(pp, edge);
01069       }  else if(shapeName == "RECTANGULAR_PRISM" && this->coordinate_dimension == 3){
01070         if(i + 6 > argCnt) {
01071           throw INVALID_SHAPE_ARG(shapeName, "6");
01072         }
01073         CoordinatePoint<T> pp;
01074         pp.x = fromString<T>(splittedStr[i++]);
01075         pp.y = fromString<T>(splittedStr[i++]);
01076         pp.z = fromString<T>(splittedStr[i++]);
01077         T edgex = fromString<T>(splittedStr[i++]);
01078         T edgey = fromString<T>(splittedStr[i++]);
01079         T edgez = fromString<T>(splittedStr[i++]);
01080         this->holes[this->holeCount++] = new RectangularPrismHole<T>(pp, edgex, edgey, edgez);
01081 
01082       }  else if(shapeName == "SPHERE" && this->coordinate_dimension == 3){
01083         if(i + 4 > argCnt) {
01084           throw INVALID_SHAPE_ARG(shapeName, "4");
01085         }
01086         CoordinatePoint<T> pp;
01087         pp.x = fromString<T>(splittedStr[i++]);
01088         pp.y = fromString<T>(splittedStr[i++]);
01089         pp.z = fromString<T>(splittedStr[i++]);
01090         T r = fromString<T>(splittedStr[i++]);
01091         this->holes[this->holeCount++] = new SphereHole<T>(pp, r);
01092       }  else {
01093         std::string tmp = toString<int>(this->coordinate_dimension);
01094         throw INVALIDSHAPE(shapeName, tmp);
01095       }
01096     }
01097     delete [] splittedStr;
01098   }
01099 
01100   void getWeightDistribution(std::string *weight_distribution_arr, int wdimension){
01101     int wcount = 0;
01102 
01103     this->wd = new WeightDistribution<T,T> *[wdimension];
01104     for(int ii = 0; ii < MAX_WEIGHT_DIM; ++ii){
01105       std::string weight_distribution = weight_distribution_arr[ii];
01106       if(weight_distribution == "") continue;
01107       if(wcount == wdimension) {
01108         throw "Weight Dimension is provided as " + toString<int>(wdimension) + ". More weight distribution is provided.";
01109       }
01110 
01111       int count = this->countChar(weight_distribution, ' ');
01112       std::string *splittedStr = new string[count + 1];
01113       this->splitString(weight_distribution, ' ', splittedStr);
01114       //cout << count << endl;
01115       T c=1;
01116       T a1=0,a2=0,a3=0;
01117       T x1=0,y1=0,z1=0;
01118       T b1=1,b2=1,b3=1;
01119       T *steps = NULL;
01120       T *values= NULL;
01121       int stepCount = 0;
01122       int valueCount = 1;
01123 
01124       for (int i = 1; i < count + 1; ++i){
01125         int assignmentCount = this->countChar(splittedStr[i], '=');
01126         if (assignmentCount != 1){
01127           throw "Error at the argument " + splittedStr[i];
01128         }
01129         std::string parameter_vs_val[2];
01130         this->splitString(splittedStr[i], '=', parameter_vs_val);
01131         std::string parameter = parameter_vs_val[0];
01132         std::string value = parameter_vs_val[1];
01133         //cout << "parameter:" << parameter << " value:" << value << endl;
01134 
01135         if (parameter == "a1"){
01136           a1 = this->fromString<T>(value);
01137         }
01138         else if (parameter == "a2"){
01139           if(this->coordinate_dimension > 1){
01140             a2 = this->fromString<T>(value);
01141           }
01142           else {
01143             throw  parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01144           }
01145 
01146         }
01147         else if (parameter == "a3"){
01148           if(this->coordinate_dimension > 2){
01149             a3 = this->fromString<T>(value);
01150           }
01151           else {
01152             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01153           }
01154         }
01155         else if (parameter == "b1"){
01156           b1 = this->fromString<T>(value);
01157         }
01158         else if (parameter == "b2"){
01159           if(this->coordinate_dimension > 1){
01160             b2 = this->fromString<T>(value);
01161           }
01162           else {
01163             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01164           }
01165         }
01166         else if (parameter == "b3"){
01167 
01168           if(this->coordinate_dimension > 2){
01169             b3 = this->fromString<T>(value);
01170           }
01171           else {
01172             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01173           }
01174         }
01175         else if (parameter == "c"){
01176           c = this->fromString<T>(value);
01177         }
01178         else if (parameter == "x1"){
01179           x1 = this->fromString<T>(value);
01180         }
01181         else if (parameter == "y1"){
01182           if(this->coordinate_dimension > 1){
01183             y1 = this->fromString<T>(value);
01184           }
01185           else {
01186             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01187           }
01188         }
01189         else if (parameter == "z1"){
01190           if(this->coordinate_dimension > 2){
01191             z1 = this->fromString<T>(value);
01192           }
01193           else {
01194             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01195           }
01196         }
01197         else if (parameter == "steps"){
01198           stepCount = this->countChar(value, ',') + 1;
01199           std::string *stepstr = new std::string[stepCount];
01200           this->splitString(value, ',', stepstr);
01201           steps = new T[stepCount];
01202           for (int i = 0; i < stepCount; ++i){
01203             steps[i] = this->fromString<T>(stepstr[i]);
01204           }
01205           delete [] stepstr;
01206         }
01207         else if (parameter == "values"){
01208           valueCount = this->countChar(value, ',') + 1;
01209           std::string *stepstr = new std::string[valueCount];
01210           this->splitString(value, ',', stepstr);
01211           values = new T[valueCount];
01212           for (int i = 0; i < valueCount; ++i){
01213             values[i] = this->fromString<T>(stepstr[i]);
01214           }
01215           delete [] stepstr;
01216         }
01217         else {
01218           throw "Invalid parameter name at " + splittedStr[i];
01219         }
01220       }
01221 
01222       delete []splittedStr;
01223       if(stepCount + 1!= valueCount){
01224         throw "Step count: " + this->toString<int>(stepCount) + " must be 1 less than value count: " + this->toString<int>(valueCount);
01225       }
01226 
01227 
01228       this->wd[wcount++] =  new SteppedEquation<T,T>(a1, a2,  a3,  b1,  b2,  b3,  c,  x1,  y1,  z1, steps, values, stepCount);
01229 
01230       if(stepCount > 0){
01231         delete [] steps;
01232         delete [] values;
01233 
01234       }
01235     }
01236     if(wcount != this->weight_dimension){
01237       throw "Weight Dimension is provided as " + toString<int>(wdimension) + ". But " + toString<int>(wcount)+" weight distributions are provided.";
01238     }
01239   }
01240 
01241   void parseParams(Teuchos::ParameterList params){
01242     try {
01243       std::string holeDescription  = "";
01244       std::string proc_load_distributions = "";
01245       std::string distinctDescription = "";
01246       std::string coordinate_distributions = "";
01247       std::string outfile = "";
01248       std::string weight_dimension_parameters[MAX_WEIGHT_DIM];
01249       for (int i = 0; i < MAX_WEIGHT_DIM; ++i){
01250         weight_dimension_parameters[i] = "";
01251       }
01252 
01253 
01254       for (Teuchos::ParameterList::ConstIterator pit = params.begin(); pit != params.end(); ++pit ){
01255         const std::string &paramName = params.name(pit);
01256 
01257         const Teuchos::ParameterEntry &pe = params.entry(pit);
01258 
01259         if(paramName.find("hole-") == 0){
01260           if(holeDescription == ""){
01261             holeDescription = getParamVal<std::string>(pe, paramName);
01262           }
01263           else {
01264             holeDescription +=","+getParamVal<std::string>(pe, paramName);
01265           }
01266         }
01267         else if(paramName.find("distribution-") == 0){
01268           if(coordinate_distributions == "")
01269             coordinate_distributions = getParamVal<std::string>(pe, paramName);
01270           else
01271             coordinate_distributions +=  ","+getParamVal<std::string>(pe, paramName);
01272           //cout << coordinate_distributions << endl;
01273           //TODO coordinate distribution description
01274         }
01275 
01276         else if (paramName.find(weight_distribution_string) == 0){
01277           std::string weight_dist_param = paramName.substr(weight_distribution_string.size());
01278           int dash_pos = weight_dist_param.find("-");
01279           std::string distribution_index_string = weight_dist_param.substr(0, dash_pos);
01280           int distribution_index = fromString<int>(distribution_index_string);
01281 
01282           if(distribution_index >= MAX_WEIGHT_DIM){
01283             throw "Given distribution index:" + distribution_index_string + " larger than maximum allowed weight dimension:" + toString<int>(MAX_WEIGHT_DIM);
01284           }
01285           weight_dimension_parameters[distribution_index] +=  " " + weight_dist_param.substr(dash_pos + 1)+ "="+ getParamVal<std::string>(pe, paramName);
01286         }
01287         else if(paramName == "dim"){
01288           int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
01289           if(dim < 2 && dim > 3){
01290             throw INVALID(paramName);
01291           } else {
01292             this->coordinate_dimension = dim;
01293           }
01294         }
01295         else if(paramName == "wdim"){
01296           int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
01297           if(dim < 1 && dim > MAX_WEIGHT_DIM){
01298             throw INVALID(paramName);
01299           } else {
01300             this->weight_dimension = dim;
01301           }
01302         }
01303 
01304         else if(paramName == "proc_load_distributions"){
01305           proc_load_distributions = getParamVal<std::string>(pe, paramName);
01306           this->loadDistSet = true;
01307         }
01308 
01309         else if(paramName == "distinct_coordinates"){
01310           distinctDescription = getParamVal<std::string>(pe, paramName);
01311           if(distinctDescription == "T"){
01312             this->distinctCoordSet = true;
01313           } else if(distinctDescription == "F"){
01314             this->distinctCoordSet = true;
01315           } else {
01316             throw "Invalid parameter for distinct_coordinates: " + distinctDescription + ". Candiates: T or F." ;
01317           }
01318         }
01319 
01320         else if(paramName == "out_file"){
01321           this->outfile = getParamVal<std::string>(pe, paramName);
01322         }
01323 
01324         else {
01325           throw INVALID(paramName);
01326         }
01327       }
01328 
01329 
01330       if(this->coordinate_dimension == 0){
01331         throw "Dimension must be provided to coordinate generator.";
01332       }
01333       /*
01334       if(this->globalNumCoords == 0){
01335         throw "Number of coordinates must be provided to coordinate generator.";
01336       }
01337        */
01338       /*
01339       if(maxx <= minx ){
01340         throw "Error: maxx= "+ toString<T>(maxx)+ " and minx=" + toString<T>(minx);
01341       }
01342       if(maxy <= miny ){
01343         throw "Error: maxy= "+ toString<T>(maxy)+ " and miny=" + toString<T>(miny);
01344 
01345       }
01346       if(this->dimension == 3 && maxz <= minz ){
01347         throw "Error: maxz= "+ toString<T>(maxz)+ " and minz=" + toString<T>(minz);
01348       }
01349        */
01350       if (this->loadDistSet && this->distinctCoordSet){
01351         throw "distinct_coordinates and proc_load_distributions parameters cannot be satisfied together.";
01352       }
01353       this->getHoles(holeDescription);
01354       //this->getDistinctCoordinateDescription(distinctDescription);
01355       this->getProcLoadDistributions(proc_load_distributions);
01356       this->getCoordinateDistributions(coordinate_distributions);
01357       this->getWeightDistribution(weight_dimension_parameters, this->weight_dimension);
01358       /*
01359       if(this->numGlobalCoords <= 0){
01360         throw "Must have at least 1 point";
01361       }
01362        */
01363     }
01364     catch(std::string s){
01365       throw s;
01366     }
01367 
01368     catch(char * s){
01369       throw s;
01370     }
01371 
01372     catch(char const* s){
01373       throw s;
01374     }
01375 
01376   }
01377 public:
01378 
01379   ~GeometricGenerator(){
01380     if(this->holes){
01381       for (int i = 0; i < this->holeCount; ++i){
01382         delete this->holes[i];
01383       }
01384       free (this->holes);
01385     }
01386 
01387 
01388     delete []loadDistributions; //sized as the number of processors, the load of each processor.
01389     //delete []distinctCoordinates; // if processors have different or same range for coordinates to be created.
01390     if(coordinateDistributions){
01391 
01392       for (int i = 0; i < this->distributionCount; ++i){
01393         delete this->coordinateDistributions[i];
01394       }
01395       free (this->coordinateDistributions);
01396     }
01397     if (this->wd){
01398       for (int i = 0; i < this->weight_dimension; ++i){
01399         delete this->wd[i];
01400       }
01401       delete []this->wd;
01402     }
01403 
01404     if(this->weight_dimension){
01405       for(int i = 0; i < this->weight_dimension; ++i)
01406       delete [] this->wghts[i];
01407       delete []this->wghts;
01408     }
01409     if(this->coordinate_dimension){
01410       for(int i = 0; i < this->coordinate_dimension; ++i)
01411       delete [] this->coords[i];
01412       delete [] this->coords;
01413     }
01414     //delete []this->points;
01415   }
01416 
01417   void print_description(){
01418     cout <<"\nGeometric Generator Parameter File Format:" << endl;
01419     cout <<"- dim=Coordinate Dimension: 2 or 3" << endl;
01420     cout <<"- Available distributions:" << endl;
01421     cout <<"\tUNIFORM: -> distribution-1=UNIFORM,NUMCOORDINATES,XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << endl;
01422     cout <<"\tGRID: -> distribution-2=GRID,XLENGTH,YLENGTH{,ZLENGTH},XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << endl;
01423     cout <<"\tNORMAL: -> distribution-3=NORMAL,XCENTER,YCENTER{,ZCENTER},XSD,YSD,{,ZSD}" << endl;
01424     cout <<"- wdim=weight_dimension: weight dimension >= 0. There should be as many weight function as weight dimension." << endl;
01425     cout <<"- Weight Equation: w = (a1 * (x - x1)^b1) + (a2 * (y - y1)^b2) + (a3 * (z - z1)^b3) + c" << endl;
01426     cout << "Parameter settings:" << endl;
01427     cout << "\tWeightDistribution-1-a1=a1 " << endl;
01428     cout << "\tWeightDistribution-1-a2=a2 " << endl;
01429     cout << "\tWeightDistribution-1-a3=a3 " << endl;
01430     cout << "\tWeightDistribution-1-b1=b1 " << endl;
01431     cout << "\tWeightDistribution-1-b2=b2 " << endl;
01432     cout << "\tWeightDistribution-1-b3=b3 " << endl;
01433     cout << "\tWeightDistribution-1-x1=x1 " << endl;
01434     cout << "\tWeightDistribution-1-y1=y1 " << endl;
01435     cout << "\tWeightDistribution-1-z1=z1 " << endl;
01436     cout << "\tWeightDistribution-1-c=c " << endl;
01437     cout << "\tIt is possible to set step function to the result of weight equation." << endl;
01438     cout << "\tWeightDistribution-1-steps=step1,step2,step3:increasing order" << endl;
01439     cout << "\tWeightDistribution-1-values=val1,val2,val3,val4." << endl;
01440     cout << "\t\tIf w < step1 -> w = val1" << endl;
01441     cout << "\t\tElse if w < step2 -> w = val2" << endl;
01442     cout << "\t\tElse if w < step3 -> w = val3" << endl;
01443     cout << "\t\tElse  -> w = val4" << endl;
01444     cout <<"- Holes:" << endl;
01445     cout << "\thole-1:SPHERE,XCENTER,YCENTER,ZCENTER,RADIUS (only for dim=3)" << endl;
01446     cout << "\thole-2:CUBE,XCENTER,YCENTER,ZCENTER,EDGE (only for dim=3)" << endl;
01447     cout << "\thole-3:RECTANGULAR_PRISM,XCENTER,YCENTER,ZCENTER,XEDGE,YEDGE,ZEDGE (only for dim=3)" << endl;
01448     cout << "\thole-4:SQUARE,XCENTER,YCENTER,EDGE (only for dim=2)" << endl;
01449     cout << "\thole-5:RECTANGLE,XCENTER,YCENTER,XEDGE,YEDGE (only for dim=2)" << endl;
01450     cout << "\thole-6:CIRCLE,XCENTER,YCENTER,RADIUS (only for dim=2)" << endl;
01451     cout << "- out_file:out_file_path : if provided output will be written to files." << endl;
01452     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." << endl;
01453 
01454   }
01455 
01456   GeometricGenerator(Teuchos::ParameterList &params, const RCP<const Teuchos::Comm<int> > & comm_){
01457     this->wd = NULL;
01458     this->holes = NULL; //to represent if there is any hole in the input
01459     this->coordinate_dimension = 0;  //dimension of the geometry
01460     this->weight_dimension = 0;  //dimension of the geometry
01461     this->worldSize = comm_->getSize(); //comminication world object.
01462     this->numGlobalCoords = 0;  //global number of coordinates requested to be created.
01463     this->loadDistributions = NULL; //sized as the number of processors, the load of each processor.
01464     //this->distinctCoordinates = NULL; // if processors have different or same range for coordinates to be created.
01465     this->coordinateDistributions = NULL;
01466     this->holeCount = 0;
01467     this->distributionCount = 0;
01468     this->outfile = "";
01469     //this->points =  NULL;
01470 
01471     /*
01472     this->minx = 0;
01473     this->maxx = 0;
01474     this->miny = 0;
01475     this->maxy = 0;
01476     this->minz = 0;
01477     this->maxz = 0;
01478      */
01479     this->loadDistSet = false;
01480     this->distinctCoordSet = false;
01481     this->myRank = comm_->getRank();
01482 
01483     try {
01484       this->parseParams(params);
01485     }
01486     catch(std::string s){
01487       if(myRank == 0){
01488         print_description();
01489       }
01490       throw s;
01491     }
01492 
01493 
01494 
01495 
01496     lno_t myPointCount = 0;
01497     this->numGlobalCoords = 0;
01498 
01499     gno_t prefixSum = 0;
01500     for(int i = 0; i < this->distributionCount; ++i){
01501       for(int ii = 0; ii < worldSize; ++ii){
01502         lno_t increment  = lno_t (this->coordinateDistributions[i]->numPoints * this->loadDistributions[ii]);
01503         if (ii < this->coordinateDistributions[i]->numPoints % worldSize){
01504           increment += 1;
01505         }
01506         this->numGlobalCoords += increment;
01507         if(ii < myRank){
01508           prefixSum += increment;
01509         }
01510       }
01511       myPointCount += lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
01512       if (myRank < this->coordinateDistributions[i]->numPoints % worldSize){
01513         myPointCount += 1;
01514       }
01515     }
01516 
01517     this->coords = new T *[this->coordinate_dimension];
01518     for(int i = 0; i < this->coordinate_dimension; ++i){
01519       this->coords[i] = new T[myPointCount];
01520     }
01521 
01522     for (int ii = 0; ii < this->coordinate_dimension; ++ii){
01523 #ifdef HAVE_ZOLTAN2_OMP
01524 #pragma omp parallel for
01525 #endif
01526       for(int i = 0; i < myPointCount; ++i){
01527         this->coords[ii][i] = 0;
01528       }
01529     }
01530 
01531     this->numLocalCoords = 0;
01532     srand (myRank + 1);
01533     for (int i = 0; i < distributionCount; ++i){
01534 
01535       lno_t requestedPointCount = lno_t(this->coordinateDistributions[i]->numPoints *  this->loadDistributions[myRank]);
01536       if (myRank < this->coordinateDistributions[i]->numPoints % worldSize){
01537         requestedPointCount += 1;
01538       }
01539       //cout << "req:" << requestedPointCount << endl;
01540       //this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->points + this->numLocalCoords, this->holes, this->holeCount,  this->loadDistributions, myRank);
01541       this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->coords, this->numLocalCoords, this->holes, this->holeCount,  this->loadDistributions, myRank);
01542       this->numLocalCoords += requestedPointCount;
01543     }
01544 
01545 
01546 
01547     if (this->distinctCoordSet){
01548       //TODO: Partition and migration.
01549     }
01550 
01551     if(this->outfile != ""){
01552 
01553       std::ofstream myfile;
01554       myfile.open ((this->outfile + toString<int>(myRank)).c_str());
01555       for(lno_t i = 0; i < this->numLocalCoords; ++i){
01556 
01557         myfile << this->coords[0][i];
01558         if(this->coordinate_dimension > 1){
01559           myfile << " " << this->coords[1][i];
01560         }
01561         if(this->coordinate_dimension > 2){
01562           myfile << " " << this->coords[2][i];
01563         }
01564         myfile << std::endl;
01565       }
01566       myfile.close();
01567     }
01568 
01569 
01570 
01571     /*
01572     Zoltan2::XpetraMultiVectorInput < Tpetra::MultiVector<T, lno_t, gno_t, node_t> > xmv (RCP < Tpetra::MultiVector<T, lno_t, gno_t, node_t> > (tmVector));
01573 
01574     RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmVector2;
01575     Zoltan2::PartitioningSolution< Tpetra::MultiVector<T, lno_t, gno_t, node_t> > solution;
01576     xmv.applyPartitioningSolution<Tpetra::MultiVector<T, lno_t, gno_t, node_t> >(this->tmVector, &tmVector2, solution);
01577      */
01578 
01579     this->wghts = new T *[this->weight_dimension];
01580     for(int i = 0; i < this->weight_dimension; ++i){
01581       this->wghts[i] = new T[this->numLocalCoords];
01582     }
01583 
01584     for(int ii = 0; ii < this->weight_dimension; ++ii){
01585       switch(this->coordinate_dimension){
01586       case 1:
01587 #ifdef HAVE_ZOLTAN2_OMP
01588 #pragma omp parallel for
01589 #endif
01590         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01591           this->wghts[ii][i] = this->wd[ii]->get1DWeight(this->coords[0][i]);
01592         }
01593         break;
01594       case 2:
01595 #ifdef HAVE_ZOLTAN2_OMP
01596 #pragma omp parallel for
01597 #endif
01598         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01599           this->wghts[ii][i] = this->wd[ii]->get2DWeight(this->coords[0][i], this->coords[1][i]);
01600         }
01601         break;
01602       case 3:
01603 #ifdef HAVE_ZOLTAN2_OMP
01604 #pragma omp parallel for
01605 #endif
01606         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01607           this->wghts[ii][i] = this->wd[ii]->get3DWeight(this->coords[0][i], this->coords[1][i], this->coords[2][i]);
01608         }
01609         break;
01610       }
01611     }
01612   }
01613 
01614   int getWeightDimension(){
01615     return this->weight_dimension;
01616   }
01617   int getCoordinateDimension(){
01618     return this->coordinate_dimension;
01619   }
01620   lno_t getNumLocalCoords(){
01621     return this->numLocalCoords;
01622   }
01623   gno_t getNumGlobalCoords(){
01624     return this->numGlobalCoords;
01625   }
01626 
01627   T **getLocalCoordinatesView(){
01628     return this->coords;
01629   }
01630 
01631   T **getLocalWeightsView(){
01632     return this->wghts;
01633   }
01634 
01635   void getLocalCoordinatesCopy( T ** c){
01636     for(int ii = 0; ii < this->coordinate_dimension; ++ii){
01637 #ifdef HAVE_ZOLTAN2_OMP
01638 #pragma omp parallel for
01639 #endif
01640       for (lno_t i = 0; i < this->numLocalCoords; ++i){
01641         c[ii][i] = this->coords[ii][i];
01642       }
01643     }
01644   }
01645 
01646   void getLocalWeightsCopy(T **w){
01647     for(int ii = 0; ii < this->weight_dimension; ++ii){
01648 #ifdef HAVE_ZOLTAN2_OMP
01649 #pragma omp parallel for
01650 #endif
01651       for (lno_t i = 0; i < this->numLocalCoords; ++i){
01652         w[ii][i] = this->wghts[ii][i];
01653       }
01654     }
01655   }
01656 };