Zoltan2 Version of the Day
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 // KDD const std::string weight_distribution = "WeightDistribution";
00079 const std::string weight_distribution_string = "WeightDistribution-";
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 wSize) :
00341     numPoints(np_), dimension(dim), requested(0), assignedPrevious(0),
00342     worldSize(wSize){}
00343 
00344   virtual CoordinatePoint<T> getPoint(gno_t point_index, unsigned int &state) = 0;
00345   virtual T getXCenter() = 0;
00346   virtual T getXRadius() =0;
00347 
00348   void GetPoints(lno_t requestedPointcount, CoordinatePoint<T> *points /*preallocated sized numPoints*/,
00349       Hole<T> **holes, lno_t holeCount,
00350       float *sharedRatios_, int myRank){
00351 
00352     for (int i = 0; i < myRank; ++i){
00353       //cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< endl;
00354       this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
00355       if (i < this->numPoints % this->worldSize){
00356         this->assignedPrevious += 1;
00357       }
00358     }
00359 
00360     this->requested = requestedPointcount;
00361 
00362     unsigned int slice =  UINT_MAX/(this->worldSize);
00363     unsigned int stateBegin = myRank * slice;
00364 
00365 //#ifdef HAVE_ZOLTAN2_OMP
00366 //#pragma omp parallel for
00367 //#endif
00368 #ifdef HAVE_ZOLTAN2_OMP
00369 #pragma omp parallel
00370 #endif
00371     {
00372       int me = 0;
00373       int tsize = 1;
00374 #ifdef HAVE_ZOLTAN2_OMP
00375       me = omp_get_thread_num();
00376       tsize = omp_get_num_threads();
00377 #endif
00378       unsigned int state = stateBegin + me * slice/(tsize);
00379 
00380 #ifdef HAVE_ZOLTAN2_OMP
00381 #pragma omp for
00382 #endif
00383       for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
00384         lno_t iteration = 0;
00385         while(1){
00386           if(++iteration > MAX_ITER_ALLOWED) {
00387             throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00388           }
00389           CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, &state);
00390 
00391           bool isInHole = false;
00392           for(lno_t i = 0; i < holeCount; ++i){
00393             if(holes[i][0].isInArea(p)){
00394               isInHole = true;
00395               break;
00396             }
00397           }
00398           if(isInHole) continue;
00399           points[cnt].x = p.x;
00400 
00401           points[cnt].y = p.y;
00402           points[cnt].z = p.z;
00403           break;
00404         }
00405       }
00406     }
00407 //#pragma omp parallel
00408       /*
00409     {
00410 
00411       lno_t cnt = 0;
00412       lno_t iteration = 0;
00413       while (cnt < requestedPointcount){
00414         if(++iteration > MAX_ITER_ALLOWED) {
00415           throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00416         }
00417         CoordinatePoint <T> p = this->getPoint();
00418 
00419         bool isInHole = false;
00420         for(lno_t i = 0; i < holeCount; ++i){
00421           if(holes[i][0].isInArea(p)){
00422             isInHole = true;
00423             break;
00424           }
00425         }
00426         if(isInHole) continue;
00427         iteration = 0;
00428         points[cnt].x = p.x;
00429         points[cnt].y = p.y;
00430         points[cnt].z = p.z;
00431         ++cnt;
00432       }
00433     }
00434     */
00435   }
00436 
00437   void GetPoints(lno_t requestedPointcount, T **coords/*preallocated sized numPoints*/, lno_t tindex,
00438       Hole<T> **holes, lno_t holeCount,
00439       float *sharedRatios_, int myRank){
00440 
00441     for (int i = 0; i < myRank; ++i){
00442       //cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< endl;
00443       this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
00444       if (i < this->numPoints % this->worldSize){
00445         this->assignedPrevious += 1;
00446       }
00447     }
00448 
00449     this->requested = requestedPointcount;
00450 
00451     unsigned int slice =  UINT_MAX/(this->worldSize);
00452     unsigned int stateBegin = myRank * slice;
00453 #ifdef HAVE_ZOLTAN2_OMP
00454 #pragma omp parallel
00455 #endif
00456     {
00457       int me = 0;
00458       int tsize = 1;
00459 #ifdef HAVE_ZOLTAN2_OMP
00460       me = omp_get_thread_num();
00461       tsize = omp_get_num_threads();
00462 #endif
00463       unsigned int state = stateBegin + me * (slice/(tsize));
00464       /*
00465 #pragma omp critical
00466       {
00467 
00468         cout << "myRank:" << me << " stateBeg:" << stateBegin << " tsize:" << tsize << " state:" << state <<  " slice: " << slice / tsize <<  endl;
00469       }
00470       */
00471 #ifdef HAVE_ZOLTAN2_OMP
00472 #pragma omp for
00473 #endif
00474       for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
00475         lno_t iteration = 0;
00476         while(1){
00477           if(++iteration > MAX_ITER_ALLOWED) {
00478             throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00479           }
00480           CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, state);
00481 
00482           bool isInHole = false;
00483           for(lno_t i = 0; i < holeCount; ++i){
00484             if(holes[i][0].isInArea(p)){
00485               isInHole = true;
00486               break;
00487             }
00488           }
00489           if(isInHole) continue;
00490           coords[0][cnt + tindex] = p.x;
00491           if(this->dimension > 1){
00492             coords[1][cnt + tindex] = p.y;
00493             if(this->dimension > 2){
00494               coords[2][cnt + tindex] = p.z;
00495             }
00496           }
00497           break;
00498         }
00499       }
00500     }
00501   }
00502 };
00503 
00504 template <typename T, typename lno_t, typename gno_t>
00505 class CoordinateNormalDistribution:public CoordinateDistribution<T,lno_t,gno_t>{
00506 public:
00507   CoordinatePoint<T> center;
00508   T standartDevx;
00509   T standartDevy;
00510   T standartDevz;
00511 
00512 
00513   virtual T getXCenter(){
00514     return center.x;
00515   }
00516   virtual T getXRadius(){
00517     return standartDevx;
00518   }
00519 
00520   CoordinateNormalDistribution(gno_t np_, int dim, CoordinatePoint<T> center_ ,
00521                                T sd_x, T sd_y, T sd_z, int wSize) : 
00522     CoordinateDistribution<T,lno_t,gno_t>(np_,dim,wSize),
00523     standartDevx(sd_x), standartDevy(sd_y), standartDevz(sd_z)
00524   {
00525     this->center.x = center_.x;
00526     this->center.y = center_.y;
00527     this->center.z = center_.z;
00528   }
00529 
00530   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00531 
00532     pindex = 0; // not used in normal distribution.
00533     CoordinatePoint <T> p;
00534 
00535     for(int i = 0; i < this->dimension; ++i){
00536       switch(i){
00537       case 0:
00538         p.x = normalDist(this->center.x, this->standartDevx, state);
00539         break;
00540       case 1:
00541         p.y = normalDist(this->center.y, this->standartDevy, state);
00542         break;
00543       case 2:
00544         p.z = normalDist(this->center.z, this->standartDevz, state);
00545         break;
00546       default:
00547         throw "unsupported dimension";
00548       }
00549     }
00550     return p;
00551   }
00552 
00553   virtual ~CoordinateNormalDistribution(){};
00554 private:
00555   T normalDist(T center_, T sd, unsigned int &state) {
00556     static bool derived=false;
00557     static T storedDerivation;
00558     T polarsqrt, normalsquared, normal1, normal2;
00559     if (!derived) {
00560       do {
00561         normal1=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
00562         normal2=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
00563         normalsquared=normal1*normal1+normal2*normal2;
00564       } while ( normalsquared>=1.0 || normalsquared == 0.0);
00565 
00566       polarsqrt=sqrt(-2.0*log(normalsquared)/normalsquared);
00567       storedDerivation=normal1*polarsqrt;
00568       derived=true;
00569       return normal2*polarsqrt*sd + center_;
00570     }
00571     else {
00572       derived=false;
00573       return storedDerivation*sd + center_;
00574     }
00575   }
00576 };
00577 
00578 template <typename T, typename lno_t, typename gno_t>
00579 class CoordinateUniformDistribution:public CoordinateDistribution<T,lno_t, gno_t>{
00580 public:
00581   T leftMostx;
00582   T rightMostx;
00583   T leftMosty;
00584   T rightMosty;
00585   T leftMostz;
00586   T rightMostz;
00587 
00588 
00589   virtual T getXCenter(){
00590     return (rightMostx - leftMostx)/2  + leftMostx;
00591   }
00592   virtual T getXRadius(){
00593     return (rightMostx - leftMostx)/2;
00594   }
00595 
00596 
00597   CoordinateUniformDistribution(gno_t np_, int dim, T l_x, T r_x, T l_y, T r_y,
00598                                 T l_z, T r_z, int wSize ) :
00599       CoordinateDistribution<T,lno_t,gno_t>(np_,dim,wSize),
00600       leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y),
00601       leftMostz(l_z), rightMostz(r_z){}
00602 
00603   virtual ~CoordinateUniformDistribution(){};
00604   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00605 
00606 
00607     pindex = 0; //not used in uniform dist.
00608     CoordinatePoint <T> p;
00609     for(int i = 0; i < this->dimension; ++i){
00610       switch(i){
00611       case 0:
00612         p.x = uniformDist(this->leftMostx, this->rightMostx, state);
00613         break;
00614       case 1:
00615         p.y = uniformDist(this->leftMosty, this->rightMosty, state);
00616         break;
00617       case 2:
00618         p.z = uniformDist(this->leftMostz, this->rightMostz, state);
00619         break;
00620       default:
00621         throw "unsupported dimension";
00622       }
00623     }
00624     return p;
00625   }
00626 
00627 private:
00628 
00629   T uniformDist(T a, T b, unsigned int &state)
00630   {
00631     return (b-a)*(rand_r(&state) / double(RAND_MAX)) + a;
00632   }
00633 };
00634 
00635 template <typename T, typename lno_t, typename gno_t>
00636 class CoordinateGridDistribution:public CoordinateDistribution<T,lno_t,gno_t>{
00637 public:
00638   T leftMostx;
00639   T rightMostx;
00640   T leftMosty;
00641   T rightMosty;
00642   T leftMostz;
00643   T rightMostz;
00644   gno_t along_X, along_Y, along_Z;
00645   //T currentX, currentY, currentZ;
00646   T processCnt;
00647   int myRank;
00648   T xstep, ystep, zstep;
00649   gno_t xshift, yshift, zshift;
00650 
00651   virtual T getXCenter(){
00652     return (rightMostx - leftMostx)/2  + leftMostx;
00653   }
00654   virtual T getXRadius(){
00655     return (rightMostx - leftMostx)/2;
00656   }
00657 
00658 
00659   CoordinateGridDistribution(gno_t alongX, gno_t alongY, gno_t alongZ, int dim, 
00660                              T l_x, T r_x, T l_y, T r_y, T l_z, T r_z , 
00661                              int myRank_, int wSize) :
00662       CoordinateDistribution<T,lno_t,gno_t>(alongX * alongY * alongZ,dim,wSize),
00663       leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y), leftMostz(l_z), rightMostz(r_z), myRank(myRank_){
00664     //currentX = leftMostx, currentY = leftMosty, currentZ = leftMostz;
00665     this->processCnt = 0;
00666     this->along_X = alongX; this->along_Y = alongY; this->along_Z = alongZ;
00667 
00668     if(this->along_X > 1)
00669       this->xstep = (rightMostx - leftMostx) / (alongX - 1);
00670     else
00671       this->xstep = 1;
00672     if(this->along_Y > 1)
00673       this->ystep = (rightMosty - leftMosty) / (alongY - 1);
00674     else
00675       this->ystep = 1;
00676     if(this->along_Z > 1)
00677       this->zstep = (rightMostz - leftMostz) / (alongZ - 1);
00678     else
00679       this->zstep = 1;
00680     xshift = 0; yshift = 0; zshift = 0;
00681   }
00682 
00683   virtual ~CoordinateGridDistribution(){};
00684   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00685     //lno_t before = processCnt + this->assignedPrevious;
00686     //cout << "before:" << processCnt << " " << this->assignedPrevious << endl;
00687     //lno_t xshift = 0, yshift = 0, zshift = 0;
00688 
00689     //lno_t tmp = before % (this->along_X * this->along_Y);
00690     //xshift  = tmp % this->along_X;
00691     //yshift = tmp / this->along_X;
00692     //zshift = before / (this->along_X * this->along_Y);
00693 
00694     state = 0; //not used here
00695     this->zshift = pindex / (along_X * along_Y);
00696     this->yshift = (pindex % (along_X * along_Y)) / along_X;
00697     this->xshift = (pindex % (along_X * along_Y)) % along_X;
00698 
00699     //this->xshift = pindex / (along_Z * along_Y);
00700    // this->zshift = (pindex % (along_Z * along_Y)) / along_Y;
00701     //this->yshift = (pindex % (along_Z * along_Y)) % along_Y;
00702 
00703     CoordinatePoint <T> p;
00704     p.x = xshift * this->xstep + leftMostx;
00705     p.y = yshift * this->ystep + leftMosty;
00706     p.z = zshift * this->zstep + leftMostz;
00707 /*
00708     ++xshift;
00709     if(xshift == this->along_X){
00710       ++yshift;
00711       xshift = 0;
00712       if(yshift == this->along_Y){
00713         ++zshift;
00714         yshift = 0;
00715       }
00716     }
00717 */
00718     /*
00719     if(this->processCnt == 0){
00720       this->xshift = this->assignedPrevious / (along_Z * along_Y);
00721       //this->yshift = (this->assignedPrevious % (along_X * along_Y)) / along_X;
00722       this->zshift = (this->assignedPrevious % (along_Z * along_Y)) / along_Y;
00723       //this->xshift = (this->assignedPrevious % (along_X * along_Y)) % along_X;
00724       this->yshift = (this->assignedPrevious % (along_Z * along_Y)) % along_Y;
00725       ++this->processCnt;
00726     }
00727 
00728     CoordinatePoint <T> p;
00729     p.x = xshift * this->xstep + leftMostx;
00730     p.y = yshift * this->ystep + leftMosty;
00731     p.z = zshift * this->zstep + leftMostz;
00732 
00733     ++yshift;
00734     if(yshift == this->along_Y){
00735       ++zshift;
00736       yshift = 0;
00737       if(zshift == this->along_Z){
00738         ++xshift;
00739         zshift = 0;
00740       }
00741 
00742     }
00743     */
00744     /*
00745     if(this->requested - 1 > this->processCnt){
00746       this->processCnt++;
00747     } else {
00748       cout << "req:" << this->requested << " pr:" <<this->processCnt << endl;
00749       cout << "p:" << p.x << " " << p.y << " " << p.z <<  endl ;
00750       cout << "s:" << xshift << " " << yshift << " " << zshift <<  endl ;
00751       cout << "st:" << this->xstep << " " << this->ystep << " " << this->zstep <<  endl ;
00752     }
00753     */
00754     return p;
00755   }
00756 
00757 private:
00758 
00759 };
00760 
00761 template <typename T, typename lno_t, typename gno_t, typename node_t>
00762 class GeometricGenerator {
00763 private:
00764   Hole<T> **holes; //to represent if there is any hole in the input
00765   int holeCount;
00766   int coordinate_dimension;  //dimension of the geometry
00767   gno_t numGlobalCoords;  //global number of coordinates requested to be created.
00768   lno_t numLocalCoords;
00769   float *loadDistributions; //sized as the number of processors, the load of each processor.
00770   bool loadDistSet;
00771   bool distinctCoordSet;
00772   CoordinateDistribution<T, lno_t,gno_t> **coordinateDistributions;
00773   int distributionCount;
00774   //CoordinatePoint<T> *points;
00775   T **coords;
00776   T **wghts;
00777   WeightDistribution<T,T> **wd;
00778   int weight_dimension;  //dimension of the geometry
00779 
00780   //RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmVector;
00781   //RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmwVector;
00782   int worldSize;
00783   int myRank;
00784   T minx;
00785   T maxx;
00786   T miny;
00787   T maxy;
00788   T minz;
00789   T maxz;
00790   std::string outfile;
00791 
00792 
00793 
00794   template <typename tt>
00795   tt getParamVal( const Teuchos::ParameterEntry& pe, const std::string &paramname){
00796     tt returnVal;
00797     try {
00798       returnVal = Teuchos::getValue<tt>(pe);
00799     }
00800     catch (...){
00801       throw INVALID(paramname);
00802     }
00803     return returnVal;
00804   }
00805 
00806   int countChar (std::string inStr, char cntChar){
00807     int cnt = 0;
00808     for (unsigned int i = 0; i < inStr.size(); ++i){
00809       if (inStr[i] == cntChar) {
00810         cnt++;
00811       }
00812     }
00813     return cnt;
00814   }
00815 
00816   template <typename tt>
00817   std::string toString(tt obj){
00818     std::stringstream ss (std::stringstream::in |std::stringstream::out);
00819     ss << obj;
00820     std::string tmp = "";
00821     ss >> tmp;
00822     return tmp;
00823   }
00824 
00825   template <typename tt>
00826   tt fromString(std::string obj){
00827     std::stringstream ss (std::stringstream::in | std::stringstream::out);
00828     ss << obj;
00829     tt tmp;
00830     ss >> tmp;
00831     if (ss.fail()){
00832       throw "Cannot convert string " + obj;
00833     }
00834     return tmp;
00835   }
00836 
00837 
00838   void splitString(std::string inStr, char splitChar, std::string *outSplittedStr){
00839     std::stringstream ss (std::stringstream::in | std::stringstream::out);
00840     ss << inStr;
00841     int cnt = 0;
00842     while (!ss.eof()){
00843       std::string tmp = "";
00844       std::getline(ss, tmp ,splitChar);
00845       outSplittedStr[cnt++] = tmp;
00846     }
00847   }
00848 
00849   /*
00850   void getDistinctCoordinateDescription(std::string distinctDescription){
00851 
00852     this->distinctCoordinates = new bool[this->dimension];
00853     if(distinctDescription != ""){
00854       int argCnt = this->countChar(distinctDescription, ',') + 1;
00855       if(argCnt != this->dimension) {
00856         throw "Invalid parameter count for distinct_coordinates. Size should be equal to dimension.";
00857       }
00858 
00859       std::string *splittedStr = new std::string[argCnt];
00860       splitString(holeDescription, ',', splittedStr);
00861 
00862       for(int i = 0; i < argCnt; ++i){
00863         if(splittedStr[i] == "T"){
00864           distinctCoordinates[i] = true;
00865         } else if(splittedStr[i] == "F"){
00866           distinctCoordinates[i] = false;
00867         } else {
00868           throw "Invalid parameter " + splittedStr[i] + " for distinct_coordinates.";
00869         }
00870       }
00871       delete []splittedStr;
00872     }
00873     else {
00874       for(int i = 0; i < this->dimension; ++i){
00875         distinctCoordinates[i] = false;
00876       }
00877     }
00878   }
00879    */
00880 
00881 
00882   void getCoordinateDistributions(std::string coordinate_distributions){
00883     if(coordinate_distributions == ""){
00884       throw "At least one distribution function must be provided to coordinate generator.";
00885     }
00886 
00887     int argCnt = this->countChar(coordinate_distributions, ',') + 1;
00888     std::string *splittedStr = new std::string[argCnt];
00889     splitString(coordinate_distributions, ',', splittedStr);
00890     coordinateDistributions = (CoordinateDistribution<T, lno_t,gno_t> **) malloc(sizeof (CoordinateDistribution<T, lno_t,gno_t> *) * 1);
00891     for(int i = 0; i < argCnt; ){
00892       coordinateDistributions = (CoordinateDistribution<T, lno_t,gno_t> **)realloc((void *)coordinateDistributions, (this->distributionCount + 1)* sizeof(CoordinateDistribution<T, lno_t,gno_t> *));
00893 
00894       std::string distName = splittedStr[i++];
00895       gno_t np_ = 0;
00896       if(distName == "NORMAL"){
00897         int reqArg = 5;
00898         if (this->coordinate_dimension == 3){
00899           reqArg = 7;
00900         }
00901         if(i + reqArg > argCnt) {
00902           std::string tmp = toString<int>(reqArg);
00903           throw INVALID_SHAPE_ARG(distName, tmp);
00904         }
00905         np_ = fromString<gno_t>(splittedStr[i++]);
00906         CoordinatePoint<T> pp;
00907 
00908         pp.x = fromString<T>(splittedStr[i++]);
00909         pp.y = fromString<T>(splittedStr[i++]);
00910         pp.z = 0;
00911         if(this->coordinate_dimension == 3){
00912           pp.z = fromString<T>(splittedStr[i++]);
00913         }
00914 
00915         T sd_x = fromString<T>(splittedStr[i++]);
00916         T sd_y = fromString<T>(splittedStr[i++]);
00917         T sd_z = 0;
00918         if(this->coordinate_dimension == 3){
00919           sd_z = fromString<T>(splittedStr[i++]);
00920         }
00921         this->coordinateDistributions[this->distributionCount++] = new CoordinateNormalDistribution<T, lno_t,gno_t>(np_, this->coordinate_dimension, pp , sd_x, sd_y, sd_z, this->worldSize );
00922 
00923       } else if(distName == "UNIFORM" ){
00924         int reqArg = 5;
00925         if (this->coordinate_dimension == 3){
00926           reqArg = 7;
00927         }
00928         if(i + reqArg > argCnt) {
00929           std::string tmp = toString<int>(reqArg);
00930           throw INVALID_SHAPE_ARG(distName, tmp);
00931         }
00932         np_ = fromString<gno_t>(splittedStr[i++]);
00933         T l_x = fromString<T>(splittedStr[i++]);
00934         T r_x = fromString<T>(splittedStr[i++]);
00935         T l_y = fromString<T>(splittedStr[i++]);
00936         T r_y = fromString<T>(splittedStr[i++]);
00937 
00938         T l_z = 0, r_z = 0;
00939 
00940         if(this->coordinate_dimension == 3){
00941           l_z = fromString<T>(splittedStr[i++]);
00942           r_z = fromString<T>(splittedStr[i++]);
00943         }
00944 
00945         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, this->worldSize );
00946       } else if (distName == "GRID"){
00947         int reqArg = 6;
00948         if(this->coordinate_dimension == 3){
00949           reqArg = 9;
00950         }
00951         if(i + reqArg > argCnt) {
00952           std::string tmp = toString<int>(reqArg);
00953           throw INVALID_SHAPE_ARG(distName, tmp);
00954         }
00955 
00956         gno_t np_x = fromString<gno_t>(splittedStr[i++]);
00957         gno_t np_y = fromString<gno_t>(splittedStr[i++]);
00958         gno_t np_z = 1;
00959 
00960 
00961         if(this->coordinate_dimension == 3){
00962           np_z = fromString<gno_t>(splittedStr[i++]);
00963         }
00964 
00965         np_ = np_x * np_y * np_z;
00966         T l_x = fromString<T>(splittedStr[i++]);
00967         T r_x = fromString<T>(splittedStr[i++]);
00968         T l_y = fromString<T>(splittedStr[i++]);
00969         T r_y = fromString<T>(splittedStr[i++]);
00970 
00971         T l_z = 0, r_z = 0;
00972 
00973         if(this->coordinate_dimension == 3){
00974           l_z = fromString<T>(splittedStr[i++]);
00975           r_z = fromString<T>(splittedStr[i++]);
00976         }
00977 
00978         if(np_x < 1 || np_z < 1 || np_y < 1 ){
00979           throw "Provide at least 1 point along each dimension for grid test.";
00980         }
00981         //cout << "ly:" << l_y << " ry:" << r_y << endl;
00982         this->coordinateDistributions[this->distributionCount++] = new CoordinateGridDistribution<T, lno_t,gno_t>
00983         (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);
00984 
00985       }
00986       else {
00987         std::string tmp = toString<int>(this->coordinate_dimension);
00988         throw INVALIDSHAPE(distName, tmp);
00989       }
00990       this->numGlobalCoords += (gno_t) np_;
00991     }
00992     delete []splittedStr;
00993   }
00994 
00995   void getProcLoadDistributions(std::string proc_load_distributions){
00996 
00997     this->loadDistributions = new float[this->worldSize];
00998     if(proc_load_distributions == ""){
00999       float r = 1.0 / this->worldSize;
01000       for(int i = 0; i < this->worldSize; ++i){
01001         this->loadDistributions[i] = r;
01002       }
01003     } else{
01004 
01005 
01006       int argCnt = this->countChar(proc_load_distributions, ',') + 1;
01007       if(argCnt != this->worldSize) {
01008         throw "Invalid parameter count load distributions. Given " + toString<int>(argCnt) + " processor size is " + toString<int>(this->worldSize);
01009       }
01010       std::string *splittedStr = new std::string[argCnt];
01011       splitString(proc_load_distributions, ',', splittedStr);
01012       for(int i = 0; i < argCnt; ++i){
01013         this->loadDistributions[i] = fromString<float>(splittedStr[i]);
01014       }
01015       delete []splittedStr;
01016 
01017 
01018       float sum = 0;
01019       for(int i = 0; i < this->worldSize; ++i){
01020         sum += this->loadDistributions[i];
01021       }
01022       if (fabs(sum - 1.0) > 10*std::numeric_limits<float>::epsilon()){
01023         throw "Processor load ratios do not sum to 1.0.";
01024       }
01025     }
01026 
01027   }
01028 
01029   void getHoles(std::string holeDescription){
01030     if(holeDescription == ""){
01031       return;
01032     }
01033     this->holes =  (Hole<T> **) malloc(sizeof (Hole <T>*));
01034     int argCnt = this->countChar(holeDescription, ',') + 1;
01035     std::string *splittedStr = new std::string[argCnt];
01036     splitString(holeDescription, ',', splittedStr);
01037 
01038     for(int i = 0; i < argCnt; ){
01039       holes = (Hole<T> **)realloc((void *)holes, (this->holeCount + 1)* sizeof(Hole<T> *));
01040 
01041       std::string shapeName = splittedStr[i++];
01042       if(shapeName == "SQUARE" && this->coordinate_dimension == 2){
01043         if(i + 3 > argCnt) {
01044           throw INVALID_SHAPE_ARG(shapeName, "3");
01045         }
01046         CoordinatePoint<T> pp;
01047         pp.x = fromString<T>(splittedStr[i++]);
01048         pp.y = fromString<T>(splittedStr[i++]);
01049         T edge = fromString<T>(splittedStr[i++]);
01050         this->holes[this->holeCount++] = new SquareHole<T>(pp, edge);
01051       } else if(shapeName == "RECTANGLE" && this->coordinate_dimension == 2){
01052         if(i + 4 > argCnt) {
01053           throw INVALID_SHAPE_ARG(shapeName, "4");
01054         }
01055         CoordinatePoint<T> pp;
01056         pp.x = fromString<T>(splittedStr[i++]);
01057         pp.y = fromString<T>(splittedStr[i++]);
01058         T edgex = fromString<T>(splittedStr[i++]);
01059         T edgey = fromString<T>(splittedStr[i++]);
01060 
01061         this->holes[this->holeCount++] = new RectangleHole<T>(pp, edgex,edgey);
01062       } else if(shapeName == "CIRCLE" && this->coordinate_dimension == 2){
01063         if(i + 3 > argCnt) {
01064           throw INVALID_SHAPE_ARG(shapeName, "3");
01065         }
01066         CoordinatePoint<T> pp;
01067         pp.x = fromString<T>(splittedStr[i++]);
01068         pp.y = fromString<T>(splittedStr[i++]);
01069         T r = fromString<T>(splittedStr[i++]);
01070         this->holes[this->holeCount++] = new CircleHole<T>(pp, r);
01071       }  else if(shapeName == "CUBE" && this->coordinate_dimension == 3){
01072         if(i + 4 > argCnt) {
01073           throw INVALID_SHAPE_ARG(shapeName, "4");
01074         }
01075         CoordinatePoint<T> pp;
01076         pp.x = fromString<T>(splittedStr[i++]);
01077         pp.y = fromString<T>(splittedStr[i++]);
01078         pp.z = fromString<T>(splittedStr[i++]);
01079         T edge = fromString<T>(splittedStr[i++]);
01080         this->holes[this->holeCount++] = new CubeHole<T>(pp, edge);
01081       }  else if(shapeName == "RECTANGULAR_PRISM" && this->coordinate_dimension == 3){
01082         if(i + 6 > argCnt) {
01083           throw INVALID_SHAPE_ARG(shapeName, "6");
01084         }
01085         CoordinatePoint<T> pp;
01086         pp.x = fromString<T>(splittedStr[i++]);
01087         pp.y = fromString<T>(splittedStr[i++]);
01088         pp.z = fromString<T>(splittedStr[i++]);
01089         T edgex = fromString<T>(splittedStr[i++]);
01090         T edgey = fromString<T>(splittedStr[i++]);
01091         T edgez = fromString<T>(splittedStr[i++]);
01092         this->holes[this->holeCount++] = new RectangularPrismHole<T>(pp, edgex, edgey, edgez);
01093 
01094       }  else if(shapeName == "SPHERE" && this->coordinate_dimension == 3){
01095         if(i + 4 > argCnt) {
01096           throw INVALID_SHAPE_ARG(shapeName, "4");
01097         }
01098         CoordinatePoint<T> pp;
01099         pp.x = fromString<T>(splittedStr[i++]);
01100         pp.y = fromString<T>(splittedStr[i++]);
01101         pp.z = fromString<T>(splittedStr[i++]);
01102         T r = fromString<T>(splittedStr[i++]);
01103         this->holes[this->holeCount++] = new SphereHole<T>(pp, r);
01104       }  else {
01105         std::string tmp = toString<int>(this->coordinate_dimension);
01106         throw INVALIDSHAPE(shapeName, tmp);
01107       }
01108     }
01109     delete [] splittedStr;
01110   }
01111 
01112   void getWeightDistribution(std::string *weight_distribution_arr, int wdimension){
01113     int wcount = 0;
01114 
01115     this->wd = new WeightDistribution<T,T> *[wdimension];
01116     for(int ii = 0; ii < MAX_WEIGHT_DIM; ++ii){
01117       std::string weight_distribution = weight_distribution_arr[ii];
01118       if(weight_distribution == "") continue;
01119       if(wcount == wdimension) {
01120         throw "Weight Dimension is provided as " + toString<int>(wdimension) + ". More weight distribution is provided.";
01121       }
01122 
01123       int count = this->countChar(weight_distribution, ' ');
01124       std::string *splittedStr = new string[count + 1];
01125       this->splitString(weight_distribution, ' ', splittedStr);
01126       //cout << count << endl;
01127       T c=1;
01128       T a1=0,a2=0,a3=0;
01129       T x1=0,y1=0,z1=0;
01130       T b1=1,b2=1,b3=1;
01131       T *steps = NULL;
01132       T *values= NULL;
01133       int stepCount = 0;
01134       int valueCount = 1;
01135 
01136       for (int i = 1; i < count + 1; ++i){
01137         int assignmentCount = this->countChar(splittedStr[i], '=');
01138         if (assignmentCount != 1){
01139           throw "Error at the argument " + splittedStr[i];
01140         }
01141         std::string parameter_vs_val[2];
01142         this->splitString(splittedStr[i], '=', parameter_vs_val);
01143         std::string parameter = parameter_vs_val[0];
01144         std::string value = parameter_vs_val[1];
01145         //cout << "parameter:" << parameter << " value:" << value << endl;
01146 
01147         if (parameter == "a1"){
01148           a1 = this->fromString<T>(value);
01149         }
01150         else if (parameter == "a2"){
01151           if(this->coordinate_dimension > 1){
01152             a2 = this->fromString<T>(value);
01153           }
01154           else {
01155             throw  parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01156           }
01157 
01158         }
01159         else if (parameter == "a3"){
01160           if(this->coordinate_dimension > 2){
01161             a3 = this->fromString<T>(value);
01162           }
01163           else {
01164             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01165           }
01166         }
01167         else if (parameter == "b1"){
01168           b1 = this->fromString<T>(value);
01169         }
01170         else if (parameter == "b2"){
01171           if(this->coordinate_dimension > 1){
01172             b2 = this->fromString<T>(value);
01173           }
01174           else {
01175             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01176           }
01177         }
01178         else if (parameter == "b3"){
01179 
01180           if(this->coordinate_dimension > 2){
01181             b3 = this->fromString<T>(value);
01182           }
01183           else {
01184             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01185           }
01186         }
01187         else if (parameter == "c"){
01188           c = this->fromString<T>(value);
01189         }
01190         else if (parameter == "x1"){
01191           x1 = this->fromString<T>(value);
01192         }
01193         else if (parameter == "y1"){
01194           if(this->coordinate_dimension > 1){
01195             y1 = this->fromString<T>(value);
01196           }
01197           else {
01198             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01199           }
01200         }
01201         else if (parameter == "z1"){
01202           if(this->coordinate_dimension > 2){
01203             z1 = this->fromString<T>(value);
01204           }
01205           else {
01206             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01207           }
01208         }
01209         else if (parameter == "steps"){
01210           stepCount = this->countChar(value, ',') + 1;
01211           std::string *stepstr = new std::string[stepCount];
01212           this->splitString(value, ',', stepstr);
01213           steps = new T[stepCount];
01214           for (int j = 0; j < stepCount; ++j){
01215             steps[j] = this->fromString<T>(stepstr[j]);
01216           }
01217           delete [] stepstr;
01218         }
01219         else if (parameter == "values"){
01220           valueCount = this->countChar(value, ',') + 1;
01221           std::string *stepstr = new std::string[valueCount];
01222           this->splitString(value, ',', stepstr);
01223           values = new T[valueCount];
01224           for (int j = 0; j < valueCount; ++j){
01225             values[j] = this->fromString<T>(stepstr[j]);
01226           }
01227           delete [] stepstr;
01228         }
01229         else {
01230           throw "Invalid parameter name at " + splittedStr[i];
01231         }
01232       }
01233 
01234       delete []splittedStr;
01235       if(stepCount + 1!= valueCount){
01236         throw "Step count: " + this->toString<int>(stepCount) + " must be 1 less than value count: " + this->toString<int>(valueCount);
01237       }
01238 
01239 
01240       this->wd[wcount++] =  new SteppedEquation<T,T>(a1, a2,  a3,  b1,  b2,  b3,  c,  x1,  y1,  z1, steps, values, stepCount);
01241 
01242       if(stepCount > 0){
01243         delete [] steps;
01244         delete [] values;
01245 
01246       }
01247     }
01248     if(wcount != this->weight_dimension){
01249       throw "Weight Dimension is provided as " + toString<int>(wdimension) + ". But " + toString<int>(wcount)+" weight distributions are provided.";
01250     }
01251   }
01252 
01253   void parseParams(Teuchos::ParameterList params){
01254     try {
01255       std::string holeDescription  = "";
01256       std::string proc_load_distributions = "";
01257       std::string distinctDescription = "";
01258       std::string coordinate_distributions = "";
01259       std::string weight_dimension_parameters[MAX_WEIGHT_DIM];
01260       for (int i = 0; i < MAX_WEIGHT_DIM; ++i){
01261         weight_dimension_parameters[i] = "";
01262       }
01263 
01264 
01265       for (Teuchos::ParameterList::ConstIterator pit = params.begin(); pit != params.end(); ++pit ){
01266         const std::string &paramName = params.name(pit);
01267 
01268         const Teuchos::ParameterEntry &pe = params.entry(pit);
01269 
01270         if(paramName.find("hole-") == 0){
01271           if(holeDescription == ""){
01272             holeDescription = getParamVal<std::string>(pe, paramName);
01273           }
01274           else {
01275             holeDescription +=","+getParamVal<std::string>(pe, paramName);
01276           }
01277         }
01278         else if(paramName.find("distribution-") == 0){
01279           if(coordinate_distributions == "")
01280             coordinate_distributions = getParamVal<std::string>(pe, paramName);
01281           else
01282             coordinate_distributions +=  ","+getParamVal<std::string>(pe, paramName);
01283           //cout << coordinate_distributions << endl;
01284           //TODO coordinate distribution description
01285         }
01286 
01287         else if (paramName.find(weight_distribution_string) == 0){
01288           std::string weight_dist_param = paramName.substr(weight_distribution_string.size());
01289           int dash_pos = weight_dist_param.find("-");
01290           std::string distribution_index_string = weight_dist_param.substr(0, dash_pos);
01291           int distribution_index = fromString<int>(distribution_index_string);
01292 
01293           if(distribution_index >= MAX_WEIGHT_DIM){
01294             throw "Given distribution index:" + distribution_index_string + " larger than maximum allowed weight dimension:" + toString<int>(MAX_WEIGHT_DIM);
01295           }
01296           weight_dimension_parameters[distribution_index] +=  " " + weight_dist_param.substr(dash_pos + 1)+ "="+ getParamVal<std::string>(pe, paramName);
01297         }
01298         else if(paramName == "dim"){
01299           int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
01300           if(dim < 2 && dim > 3){
01301             throw INVALID(paramName);
01302           } else {
01303             this->coordinate_dimension = dim;
01304           }
01305         }
01306         else if(paramName == "wdim"){
01307           int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
01308           if(dim < 1 && dim > MAX_WEIGHT_DIM){
01309             throw INVALID(paramName);
01310           } else {
01311             this->weight_dimension = dim;
01312           }
01313         }
01314 
01315         else if(paramName == "proc_load_distributions"){
01316           proc_load_distributions = getParamVal<std::string>(pe, paramName);
01317           this->loadDistSet = true;
01318         }
01319 
01320         else if(paramName == "distinct_coordinates"){
01321           distinctDescription = getParamVal<std::string>(pe, paramName);
01322           if(distinctDescription == "T"){
01323             this->distinctCoordSet = true;
01324           } else if(distinctDescription == "F"){
01325             this->distinctCoordSet = true;
01326           } else {
01327             throw "Invalid parameter for distinct_coordinates: " + distinctDescription + ". Candiates: T or F." ;
01328           }
01329         }
01330 
01331         else if(paramName == "out_file"){
01332           this->outfile = getParamVal<std::string>(pe, paramName);
01333         }
01334 
01335         else {
01336           throw INVALID(paramName);
01337         }
01338       }
01339 
01340 
01341       if(this->coordinate_dimension == 0){
01342         throw "Dimension must be provided to coordinate generator.";
01343       }
01344       /*
01345       if(this->globalNumCoords == 0){
01346         throw "Number of coordinates must be provided to coordinate generator.";
01347       }
01348        */
01349       /*
01350       if(maxx <= minx ){
01351         throw "Error: maxx= "+ toString<T>(maxx)+ " and minx=" + toString<T>(minx);
01352       }
01353       if(maxy <= miny ){
01354         throw "Error: maxy= "+ toString<T>(maxy)+ " and miny=" + toString<T>(miny);
01355 
01356       }
01357       if(this->dimension == 3 && maxz <= minz ){
01358         throw "Error: maxz= "+ toString<T>(maxz)+ " and minz=" + toString<T>(minz);
01359       }
01360        */
01361       if (this->loadDistSet && this->distinctCoordSet){
01362         throw "distinct_coordinates and proc_load_distributions parameters cannot be satisfied together.";
01363       }
01364       this->getHoles(holeDescription);
01365       //this->getDistinctCoordinateDescription(distinctDescription);
01366       this->getProcLoadDistributions(proc_load_distributions);
01367       this->getCoordinateDistributions(coordinate_distributions);
01368       this->getWeightDistribution(weight_dimension_parameters, this->weight_dimension);
01369       /*
01370       if(this->numGlobalCoords <= 0){
01371         throw "Must have at least 1 point";
01372       }
01373        */
01374     }
01375     catch(std::string s){
01376       throw s;
01377     }
01378 
01379     catch(char * s){
01380       throw s;
01381     }
01382 
01383     catch(char const* s){
01384       throw s;
01385     }
01386 
01387   }
01388 public:
01389 
01390   ~GeometricGenerator(){
01391     if(this->holes){
01392       for (int i = 0; i < this->holeCount; ++i){
01393         delete this->holes[i];
01394       }
01395       free (this->holes);
01396     }
01397 
01398 
01399     delete []loadDistributions; //sized as the number of processors, the load of each processor.
01400     //delete []distinctCoordinates; // if processors have different or same range for coordinates to be created.
01401     if(coordinateDistributions){
01402 
01403       for (int i = 0; i < this->distributionCount; ++i){
01404         delete this->coordinateDistributions[i];
01405       }
01406       free (this->coordinateDistributions);
01407     }
01408     if (this->wd){
01409       for (int i = 0; i < this->weight_dimension; ++i){
01410         delete this->wd[i];
01411       }
01412       delete []this->wd;
01413     }
01414 
01415     if(this->weight_dimension){
01416       for(int i = 0; i < this->weight_dimension; ++i)
01417       delete [] this->wghts[i];
01418       delete []this->wghts;
01419     }
01420     if(this->coordinate_dimension){
01421       for(int i = 0; i < this->coordinate_dimension; ++i)
01422       delete [] this->coords[i];
01423       delete [] this->coords;
01424     }
01425     //delete []this->points;
01426   }
01427 
01428   void print_description(){
01429     cout <<"\nGeometric Generator Parameter File Format:" << endl;
01430     cout <<"- dim=Coordinate Dimension: 2 or 3" << endl;
01431     cout <<"- Available distributions:" << endl;
01432     cout <<"\tUNIFORM: -> distribution-1=UNIFORM,NUMCOORDINATES,XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << endl;
01433     cout <<"\tGRID: -> distribution-2=GRID,XLENGTH,YLENGTH{,ZLENGTH},XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << endl;
01434     cout <<"\tNORMAL: -> distribution-3=NORMAL,XCENTER,YCENTER{,ZCENTER},XSD,YSD,{,ZSD}" << endl;
01435     cout <<"- wdim=weight_dimension: weight dimension >= 0. There should be as many weight function as weight dimension." << endl;
01436     cout <<"- Weight Equation: w = (a1 * (x - x1)^b1) + (a2 * (y - y1)^b2) + (a3 * (z - z1)^b3) + c" << endl;
01437     cout << "Parameter settings:" << endl;
01438     cout << "\tWeightDistribution-1-a1=a1 " << endl;
01439     cout << "\tWeightDistribution-1-a2=a2 " << endl;
01440     cout << "\tWeightDistribution-1-a3=a3 " << endl;
01441     cout << "\tWeightDistribution-1-b1=b1 " << endl;
01442     cout << "\tWeightDistribution-1-b2=b2 " << endl;
01443     cout << "\tWeightDistribution-1-b3=b3 " << endl;
01444     cout << "\tWeightDistribution-1-x1=x1 " << endl;
01445     cout << "\tWeightDistribution-1-y1=y1 " << endl;
01446     cout << "\tWeightDistribution-1-z1=z1 " << endl;
01447     cout << "\tWeightDistribution-1-c=c " << endl;
01448     cout << "\tIt is possible to set step function to the result of weight equation." << endl;
01449     cout << "\tWeightDistribution-1-steps=step1,step2,step3:increasing order" << endl;
01450     cout << "\tWeightDistribution-1-values=val1,val2,val3,val4." << endl;
01451     cout << "\t\tIf w < step1 -> w = val1" << endl;
01452     cout << "\t\tElse if w < step2 -> w = val2" << endl;
01453     cout << "\t\tElse if w < step3 -> w = val3" << endl;
01454     cout << "\t\tElse  -> w = val4" << endl;
01455     cout <<"- Holes:" << endl;
01456     cout << "\thole-1:SPHERE,XCENTER,YCENTER,ZCENTER,RADIUS (only for dim=3)" << endl;
01457     cout << "\thole-2:CUBE,XCENTER,YCENTER,ZCENTER,EDGE (only for dim=3)" << endl;
01458     cout << "\thole-3:RECTANGULAR_PRISM,XCENTER,YCENTER,ZCENTER,XEDGE,YEDGE,ZEDGE (only for dim=3)" << endl;
01459     cout << "\thole-4:SQUARE,XCENTER,YCENTER,EDGE (only for dim=2)" << endl;
01460     cout << "\thole-5:RECTANGLE,XCENTER,YCENTER,XEDGE,YEDGE (only for dim=2)" << endl;
01461     cout << "\thole-6:CIRCLE,XCENTER,YCENTER,RADIUS (only for dim=2)" << endl;
01462     cout << "- out_file:out_file_path : if provided output will be written to files." << endl;
01463     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;
01464 
01465   }
01466 
01467   GeometricGenerator(Teuchos::ParameterList &params, const RCP<const Teuchos::Comm<int> > & comm_){
01468     this->wd = NULL;
01469     this->holes = NULL; //to represent if there is any hole in the input
01470     this->coordinate_dimension = 0;  //dimension of the geometry
01471     this->weight_dimension = 0;  //dimension of the geometry
01472     this->worldSize = comm_->getSize(); //comminication world object.
01473     this->numGlobalCoords = 0;  //global number of coordinates requested to be created.
01474     this->loadDistributions = NULL; //sized as the number of processors, the load of each processor.
01475     //this->distinctCoordinates = NULL; // if processors have different or same range for coordinates to be created.
01476     this->coordinateDistributions = NULL;
01477     this->holeCount = 0;
01478     this->distributionCount = 0;
01479     this->outfile = "";
01480     //this->points =  NULL;
01481 
01482     /*
01483     this->minx = 0;
01484     this->maxx = 0;
01485     this->miny = 0;
01486     this->maxy = 0;
01487     this->minz = 0;
01488     this->maxz = 0;
01489      */
01490     this->loadDistSet = false;
01491     this->distinctCoordSet = false;
01492     this->myRank = comm_->getRank();
01493 
01494     try {
01495       this->parseParams(params);
01496     }
01497     catch(std::string s){
01498       if(myRank == 0){
01499         print_description();
01500       }
01501       throw s;
01502     }
01503 
01504 
01505 
01506 
01507     lno_t myPointCount = 0;
01508     this->numGlobalCoords = 0;
01509 
01510     gno_t prefixSum = 0;
01511     for(int i = 0; i < this->distributionCount; ++i){
01512       for(int ii = 0; ii < this->worldSize; ++ii){
01513         lno_t increment  = lno_t (this->coordinateDistributions[i]->numPoints * this->loadDistributions[ii]);
01514         if (ii < this->coordinateDistributions[i]->numPoints % this->worldSize){
01515           increment += 1;
01516         }
01517         this->numGlobalCoords += increment;
01518         if(ii < myRank){
01519           prefixSum += increment;
01520         }
01521       }
01522       myPointCount += lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
01523       if (myRank < this->coordinateDistributions[i]->numPoints % this->worldSize){
01524         myPointCount += 1;
01525       }
01526     }
01527 
01528     this->coords = new T *[this->coordinate_dimension];
01529     for(int i = 0; i < this->coordinate_dimension; ++i){
01530       this->coords[i] = new T[myPointCount];
01531     }
01532 
01533     for (int ii = 0; ii < this->coordinate_dimension; ++ii){
01534 #ifdef HAVE_ZOLTAN2_OMP
01535 #pragma omp parallel for
01536 #endif
01537       for(int i = 0; i < myPointCount; ++i){
01538         this->coords[ii][i] = 0;
01539       }
01540     }
01541 
01542     this->numLocalCoords = 0;
01543     srand (myRank + 1);
01544     for (int i = 0; i < distributionCount; ++i){
01545 
01546       lno_t requestedPointCount = lno_t(this->coordinateDistributions[i]->numPoints *  this->loadDistributions[myRank]);
01547       if (myRank < this->coordinateDistributions[i]->numPoints % this->worldSize){
01548         requestedPointCount += 1;
01549       }
01550       //cout << "req:" << requestedPointCount << endl;
01551       //this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->points + this->numLocalCoords, this->holes, this->holeCount,  this->loadDistributions, myRank);
01552       this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->coords, this->numLocalCoords, this->holes, this->holeCount,  this->loadDistributions, myRank);
01553       this->numLocalCoords += requestedPointCount;
01554     }
01555 
01556 
01557 
01558     if (this->distinctCoordSet){
01559       //TODO: Partition and migration.
01560     }
01561 
01562     if(this->outfile != ""){
01563 
01564       std::ofstream myfile;
01565       myfile.open ((this->outfile + toString<int>(myRank)).c_str());
01566       for(lno_t i = 0; i < this->numLocalCoords; ++i){
01567 
01568         myfile << this->coords[0][i];
01569         if(this->coordinate_dimension > 1){
01570           myfile << " " << this->coords[1][i];
01571         }
01572         if(this->coordinate_dimension > 2){
01573           myfile << " " << this->coords[2][i];
01574         }
01575         myfile << std::endl;
01576       }
01577       myfile.close();
01578     }
01579 
01580 
01581 
01582     /*
01583     Zoltan2::XpetraMultiVectorInput < Tpetra::MultiVector<T, lno_t, gno_t, node_t> > xmv (RCP < Tpetra::MultiVector<T, lno_t, gno_t, node_t> > (tmVector));
01584 
01585     RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmVector2;
01586     Zoltan2::PartitioningSolution< Tpetra::MultiVector<T, lno_t, gno_t, node_t> > solution;
01587     xmv.applyPartitioningSolution<Tpetra::MultiVector<T, lno_t, gno_t, node_t> >(this->tmVector, &tmVector2, solution);
01588      */
01589 
01590     this->wghts = new T *[this->weight_dimension];
01591     for(int i = 0; i < this->weight_dimension; ++i){
01592       this->wghts[i] = new T[this->numLocalCoords];
01593     }
01594 
01595     for(int ii = 0; ii < this->weight_dimension; ++ii){
01596       switch(this->coordinate_dimension){
01597       case 1:
01598 #ifdef HAVE_ZOLTAN2_OMP
01599 #pragma omp parallel for
01600 #endif
01601         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01602           this->wghts[ii][i] = this->wd[ii]->get1DWeight(this->coords[0][i]);
01603         }
01604         break;
01605       case 2:
01606 #ifdef HAVE_ZOLTAN2_OMP
01607 #pragma omp parallel for
01608 #endif
01609         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01610           this->wghts[ii][i] = this->wd[ii]->get2DWeight(this->coords[0][i], this->coords[1][i]);
01611         }
01612         break;
01613       case 3:
01614 #ifdef HAVE_ZOLTAN2_OMP
01615 #pragma omp parallel for
01616 #endif
01617         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01618           this->wghts[ii][i] = this->wd[ii]->get3DWeight(this->coords[0][i], this->coords[1][i], this->coords[2][i]);
01619         }
01620         break;
01621       }
01622     }
01623   }
01624 
01625   int getWeightDimension(){
01626     return this->weight_dimension;
01627   }
01628   int getCoordinateDimension(){
01629     return this->coordinate_dimension;
01630   }
01631   lno_t getNumLocalCoords(){
01632     return this->numLocalCoords;
01633   }
01634   gno_t getNumGlobalCoords(){
01635     return this->numGlobalCoords;
01636   }
01637 
01638   T **getLocalCoordinatesView(){
01639     return this->coords;
01640   }
01641 
01642   T **getLocalWeightsView(){
01643     return this->wghts;
01644   }
01645 
01646   void getLocalCoordinatesCopy( T ** c){
01647     for(int ii = 0; ii < this->coordinate_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         c[ii][i] = this->coords[ii][i];
01653       }
01654     }
01655   }
01656 
01657   void getLocalWeightsCopy(T **w){
01658     for(int ii = 0; ii < this->weight_dimension; ++ii){
01659 #ifdef HAVE_ZOLTAN2_OMP
01660 #pragma omp parallel for
01661 #endif
01662       for (lno_t i = 0; i < this->numLocalCoords; ++i){
01663         w[ii][i] = this->wghts[ii][i];
01664       }
01665     }
01666   }
01667 };