Zoltan2
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_XpetraMultiVectorAdapter.hpp>
00059 #include <Zoltan2_PartitioningSolution.hpp>
00060 #include <Teuchos_ArrayViewDecl.hpp>
00061 #include <Teuchos_RCP.hpp>
00062 #include <Tpetra_Distributor.hpp>
00063 #include <Zoltan2_PartitioningProblem.hpp>
00064 
00065 
00066 //#define HAVE_ZOLTAN2_ZOLTAN
00067 #ifdef HAVE_ZOLTAN2_ZOLTAN
00068 #include <zoltan.h>
00069 #endif
00070 
00071 using Teuchos::CommandLineProcessor;
00072 
00073 namespace GeometricGen{
00074 #define CATCH_EXCEPTIONS(pp) \
00075         catch (std::runtime_error &e) { \
00076             cout << "Runtime exception returned from " << pp << ": " \
00077             << e.what() << " FAIL" << endl; \
00078             return -1; \
00079         } \
00080         catch (std::logic_error &e) { \
00081             cout << "Logic exception returned from " << pp << ": " \
00082             << e.what() << " FAIL" << endl; \
00083             return -1; \
00084         } \
00085         catch (std::bad_alloc &e) { \
00086             cout << "Bad_alloc exception returned from " << pp << ": " \
00087             << e.what() << " FAIL" << endl; \
00088             return -1; \
00089         } \
00090         catch (std::exception &e) { \
00091             cout << "Unknown exception returned from " << pp << ": " \
00092             << e.what() << " FAIL" << endl; \
00093             return -1; \
00094         }
00095 
00096 
00097 
00098 
00099 #ifdef HAVE_ZOLTAN2_ZOLTAN
00100 
00101 
00102 template <typename tMVector_t>
00103 class DOTS{
00104 public:
00105   vector<vector<float> > weights;
00106   tMVector_t *coordinates;
00107 };
00108 
00109 template <typename tMVector_t>
00110 int getNumObj(void *data, int *ierr)
00111 {
00112   *ierr = 0;
00113   DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
00114   return dots_->coordinates->getLocalLength();
00115 }
00117 template <typename tMVector_t>
00118 void getCoords(void *data, int numGid, int numLid,
00119   int numObj, ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids,
00120   int dim, double *coords_, int *ierr)
00121 {
00122     // I know that Zoltan asks for coordinates in gid order.
00123     if (dim == 3){
00124       *ierr = 0;
00125       DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
00126       double *val = coords_;
00127       const scalar_t *x = dots_->coordinates->getData(0).getRawPtr();
00128       const scalar_t *y = dots_->coordinates->getData(1).getRawPtr();
00129       const scalar_t *z = dots_->coordinates->getData(2).getRawPtr();
00130       for (int i=0; i < numObj; i++){
00131         *val++ = static_cast<double>(x[i]);
00132         *val++ = static_cast<double>(y[i]);
00133         *val++ = static_cast<double>(z[i]);
00134       }
00135     }
00136     else {
00137       *ierr = 0;
00138       DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
00139       double *val = coords_;
00140       const scalar_t *x = dots_->coordinates->getData(0).getRawPtr();
00141       const scalar_t *y = dots_->coordinates->getData(1).getRawPtr();
00142       for (int i=0; i < numObj; i++){
00143         *val++ = static_cast<double>(x[i]);
00144         *val++ = static_cast<double>(y[i]);
00145       }
00146 
00147 
00148     }
00149 }
00150 template <typename tMVector_t>
00151 int getDim(void *data, int *ierr)
00152 {
00153   *ierr = 0;
00154   DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
00155   int dim =  dots_->coordinates->getNumVectors();
00156 
00157   return dim;
00158 }
00159 
00161 template <typename tMVector_t>
00162 void getObjList(void *data, int numGid, int numLid,
00163   ZOLTAN_ID_PTR gids, ZOLTAN_ID_PTR lids,
00164   int num_wgts, float *obj_wgts, int *ierr)
00165 {
00166   *ierr = 0;
00167   DOTS<tMVector_t> *dots_ = (DOTS<tMVector_t> *) data;
00168 
00169   size_t localLen = dots_->coordinates->getLocalLength();
00170   const gno_t *ids =
00171                dots_->coordinates->getMap()->getNodeElementList().getRawPtr();
00172 
00173   if (sizeof(ZOLTAN_ID_TYPE) == sizeof(gno_t))
00174     memcpy(gids, ids, sizeof(ZOLTAN_ID_TYPE) * localLen);
00175   else
00176     for (size_t i=0; i < localLen; i++)
00177       gids[i] = static_cast<ZOLTAN_ID_TYPE>(ids[i]);
00178 
00179   if (num_wgts > 0){
00180     float *wgts = obj_wgts;
00181     for (size_t i=0; i < localLen; i++)
00182       for (int w=0; w < num_wgts; w++)
00183         *wgts++ = dots_->weights[w][i];
00184   }
00185 }
00186 #endif
00187 
00188 
00189 enum shape {SQUARE, RECTANGLE, CIRCLE, CUBE, RECTANGULAR_PRISM, SPHERE};
00190 const std::string shapes[] = {"SQUARE", "RECTANGLE", "CIRCLE", "CUBE", "RECTANGULAR_PRISM", "SPHERE"};
00191 #define SHAPE_COUNT 6
00192 
00193 enum distribution {normal, uniform};
00194 const std::string distribution[] = {"distribution", "uniform"};
00195 #define DISTRIBUTION_COUNT 2
00196 
00197 #define HOLE_ALLOC_STEP  10
00198 #define MAX_WEIGHT_DIM  10
00199 #define INVALID(STR) "Invalid argument at " + STR
00200 #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"
00201 
00202 #define INVALID_SHAPE_ARG(SHAPE, REQUIRED) "Invalid argument count for shape " + SHAPE + ". Requires " + REQUIRED + " argument(s)."
00203 #define MAX_ITER_ALLOWED 500
00204 
00205 // KDD const std::string weight_distribution = "WeightDistribution";
00206 const std::string weight_distribution_string = "WeightDistribution-";
00207 template <typename T>
00208 struct CoordinatePoint {
00209   T x;
00210   T y;
00211   T z;
00212   /*
00213   bool isInArea(int dim, T *minCoords, T *maxCoords){
00214     dim = min(3, dim);
00215     for (int i = 0; i < dim; ++i){
00216       bool maybe = true;
00217       switch(i){
00218       case 0:
00219         maybe = x < maxCoords[i] && x > maxCoords[i];
00220         break;
00221       case 1:
00222         maybe = y < maxCoords[i] && y > maxCoords[i];
00223         break;
00224       case 2:
00225         maybe = z < maxCoords[i] && z > maxCoords[i];
00226         break;
00227       }
00228       if (!maybe){
00229         return false;
00230       }
00231     }
00232     return true;
00233   }
00234    */
00235   CoordinatePoint(){
00236     x = 0;y=0;z=0;
00237   }
00238 };
00239 
00240 template <typename T>
00241 class Hole{
00242 public:
00243   CoordinatePoint<T> center;
00244   T edge1, edge2, edge3;
00245   Hole(CoordinatePoint<T> center_, T edge1_, T edge2_, T edge3_){
00246     this->center.x = center_.x;
00247     this->center.y = center_.y;
00248     this->center.z = center_.z;
00249     this->edge1 = edge1_;
00250     this->edge2 = edge2_;
00251     this->edge3 = edge3_;
00252   }
00253   virtual bool isInArea(CoordinatePoint<T> dot) = 0;
00254   virtual ~Hole(){}
00255 };
00256 
00257 template <typename T>
00258 class SquareHole:public Hole<T>{
00259 public:
00260   SquareHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00261   virtual ~SquareHole(){}
00262   virtual bool isInArea(CoordinatePoint<T> dot){
00263     return fabs(dot.x - this->center.x) < this->edge1 / 2 && fabs(dot.y - this->center.y) < this->edge1 / 2;
00264   }
00265 };
00266 
00267 template <typename T>
00268 class RectangleHole:public Hole<T>{
00269 public:
00270   RectangleHole(CoordinatePoint<T> center_  , T edge_x_, T edge_y_): Hole<T>(center_, edge_x_,  edge_y_, 0){}
00271   virtual bool isInArea(CoordinatePoint<T> dot){
00272     return fabs(dot.x - this->center.x) < this->edge1 / 2 && fabs(dot.y - this->center.y) < this->edge2 / 2;
00273   }
00274   virtual ~RectangleHole(){}
00275 };
00276 
00277 template <typename T>
00278 class CircleHole:public Hole<T>{
00279 public:
00280   CircleHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00281   virtual ~CircleHole(){}
00282   virtual bool isInArea(CoordinatePoint<T> dot){
00283     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;
00284   }
00285 };
00286 
00287 template <typename T>
00288 class CubeHole:public Hole<T>{
00289 public:
00290   CubeHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00291   virtual ~CubeHole(){}
00292   virtual bool isInArea(CoordinatePoint<T> dot){
00293     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;
00294   }
00295 };
00296 
00297 template <typename T>
00298 class RectangularPrismHole:public Hole<T>{
00299 public:
00300   RectangularPrismHole(CoordinatePoint<T> center_  , T edge_x_, T edge_y_, T edge_z_): Hole<T>(center_, edge_x_,  edge_y_, edge_z_){}
00301   virtual ~RectangularPrismHole(){}
00302   virtual bool isInArea(CoordinatePoint<T> dot){
00303     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;
00304   }
00305 };
00306 
00307 template <typename T>
00308 class SphereHole:public Hole<T>{
00309 public:
00310   SphereHole(CoordinatePoint<T> center_ , T edge_): Hole<T>(center_, edge_, 0 , 0){}
00311   virtual ~SphereHole(){}
00312   virtual bool isInArea(CoordinatePoint<T> dot){
00313     return  fabs((dot.x - this->center.x) * (dot.x - this->center.x) * (dot.x - this->center.x)) +
00314         fabs((dot.y - this->center.y) * (dot.y - this->center.y) * (dot.y - this->center.y)) +
00315         fabs((dot.z - this->center.z) * (dot.z - this->center.z) * (dot.z - this->center.z))
00316         <
00317         this->edge1 * this->edge1 * this->edge1;
00318   }
00319 };
00320 
00321 template <typename T, typename weighttype>
00322 class WeightDistribution{
00323 public:
00324   virtual weighttype getWeight(CoordinatePoint<T> P)=0;
00325   virtual weighttype get1DWeight(T x)=0;
00326   virtual weighttype get2DWeight(T x, T y)=0;
00327   virtual weighttype get3DWeight(T x, T y, T z)=0;
00328   WeightDistribution(){};
00329   virtual ~WeightDistribution(){};
00330 };
00331 
00350 template <typename T, typename weighttype>
00351 class SteppedEquation:public WeightDistribution<T, weighttype>{
00352   T a1,a2,a3;
00353   T b1,b2,b3;
00354   T c;
00355   T x1,y1,z1;
00356 
00357   T *steps;
00358   weighttype *values;
00359   int stepCount;
00360 public:
00361   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>(){
00362     this->a1 = a1_;
00363     this->a2 = a2_;
00364     this->a3 = a3_;
00365     this->b1 = b1_;
00366     this->b2 = b2_;
00367     this->b3 = b3_;
00368     this->c = c_;
00369     this->x1 = x1_;
00370     this->y1 = y1_;
00371     this->z1 = z1_;
00372 
00373 
00374     this->stepCount = stepCount_;
00375     if(this->stepCount > 0){
00376       this->steps = new T[this->stepCount];
00377       this->values = new T[this->stepCount + 1];
00378 
00379       for (int i = 0; i < this->stepCount; ++i){
00380         this->steps[i] = steps_[i];
00381         this->values[i] = values_[i];
00382       }
00383       this->values[this->stepCount] = values_[this->stepCount];
00384     }
00385 
00386   }
00387 
00388   virtual ~SteppedEquation(){
00389     if(this->stepCount > 0){
00390       delete [] this->steps;
00391       delete [] this->values;
00392     }
00393   }
00394 
00395 
00396   virtual weighttype get1DWeight(T x){
00397     T expressionRes = this->a1 * pow( (x - this->x1), b1) + c;
00398     if(this->stepCount > 0){
00399       for (int i = 0; i < this->stepCount; ++i){
00400         if (expressionRes < this->steps[i]) return this->values[i];
00401       }
00402       return this->values[this->stepCount];
00403     }
00404     else {
00405       return weighttype(expressionRes);
00406     }
00407   }
00408 
00409   virtual weighttype get2DWeight(T x, T y){
00410     T expressionRes = this->a1 * pow( (x - this->x1), b1) + this->a2 * pow( (y - this->y1), b2) + c;
00411     if(this->stepCount > 0){
00412       for (int i = 0; i < this->stepCount; ++i){
00413         if (expressionRes < this->steps[i]) return this->values[i];
00414       }
00415       return this->values[this->stepCount];
00416     }
00417     else {
00418       return weighttype(expressionRes);
00419     }
00420   }
00421 
00422   void print (T x, T y, T z){
00423     cout << this->a1 << "*" <<  "math.pow( (" << x  << "- " <<  this->x1 << " ), " << b1 <<")" << "+" <<  this->a2<< "*"<<  "math.pow( (" << y << "-" <<  this->y1 << "), " <<
00424         b2 << " ) + " << this->a3 << " * math.pow( (" << z << "-" <<  this->z1 << "), " << b3 << ")+ " << c << " == " <<
00425         this->a1 * pow( (x - this->x1), b1) + this->a2 * pow( (y - this->y1), b2) + this->a3 * pow( (z - this->z1), b3) + c << endl;
00426 
00427   }
00428 
00429   virtual weighttype get3DWeight(T x, T y, T z){
00430     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;
00431 
00432     //this->print(x,y,z);
00433     if(this->stepCount > 0){
00434       for (int i = 0; i < this->stepCount; ++i){
00435         if (expressionRes < this->steps[i]) {
00436           //cout << "0exp:" << expressionRes << " step:" << steps[i] << " value:" << values[i] << endl;
00437           return this->values[i];
00438         }
00439       }
00440 
00441       //cout << "1exp:" << expressionRes << " step:" << steps[stepCount] << " value:" << values[stepCount] << endl;
00442       return this->values[this->stepCount];
00443     }
00444     else {
00445       return weighttype(expressionRes);
00446     }
00447   }
00448 
00449   virtual weighttype getWeight(CoordinatePoint<T> p){
00450     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);
00451     if(this->stepCount > 0){
00452       for (int i = 0; i < this->stepCount; ++i){
00453         if (expressionRes < this->steps[i]) return this->values[i];
00454       }
00455       return this->values[this->stepCount];
00456     }
00457     else {
00458       return weighttype(expressionRes);
00459     }
00460   }
00461 };
00462 
00463 template <typename T, typename lno_t, typename gno_t>
00464 class CoordinateDistribution{
00465 public:
00466   gno_t numPoints;
00467   int dimension;
00468   lno_t requested;
00469   gno_t assignedPrevious;
00470   int worldSize;
00471   virtual ~CoordinateDistribution(){}
00472 
00473   CoordinateDistribution(gno_t np_, int dim, int wSize) :
00474     numPoints(np_), dimension(dim), requested(0), assignedPrevious(0),
00475     worldSize(wSize){}
00476 
00477   virtual CoordinatePoint<T> getPoint(gno_t point_index, unsigned int &state) = 0;
00478   virtual T getXCenter() = 0;
00479   virtual T getXRadius() =0;
00480 
00481   void GetPoints(lno_t requestedPointcount, CoordinatePoint<T> *points /*preallocated sized numPoints*/,
00482       Hole<T> **holes, lno_t holeCount,
00483       float *sharedRatios_, int myRank){
00484 
00485     for (int i = 0; i < myRank; ++i){
00486       //cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< endl;
00487       this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
00488       if (i < this->numPoints % this->worldSize){
00489         this->assignedPrevious += 1;
00490       }
00491     }
00492 
00493     this->requested = requestedPointcount;
00494 
00495     unsigned int slice =  UINT_MAX/(this->worldSize);
00496     unsigned int stateBegin = myRank * slice;
00497 
00498 //#ifdef HAVE_ZOLTAN2_OMP
00499 //#pragma omp parallel for
00500 //#endif
00501 #ifdef HAVE_ZOLTAN2_OMP
00502 #pragma omp parallel
00503 #endif
00504     {
00505       int me = 0;
00506       int tsize = 1;
00507 #ifdef HAVE_ZOLTAN2_OMP
00508       me = omp_get_thread_num();
00509       tsize = omp_get_num_threads();
00510 #endif
00511       unsigned int state = stateBegin + me * slice/(tsize);
00512 
00513 #ifdef HAVE_ZOLTAN2_OMP
00514 #pragma omp for
00515 #endif
00516       for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
00517         lno_t iteration = 0;
00518         while(1){
00519           if(++iteration > MAX_ITER_ALLOWED) {
00520             throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00521           }
00522           CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, &state);
00523 
00524           bool isInHole = false;
00525           for(lno_t i = 0; i < holeCount; ++i){
00526             if(holes[i][0].isInArea(p)){
00527               isInHole = true;
00528               break;
00529             }
00530           }
00531           if(isInHole) continue;
00532           points[cnt].x = p.x;
00533 
00534           points[cnt].y = p.y;
00535           points[cnt].z = p.z;
00536           break;
00537         }
00538       }
00539     }
00540 //#pragma omp parallel
00541       /*
00542     {
00543 
00544       lno_t cnt = 0;
00545       lno_t iteration = 0;
00546       while (cnt < requestedPointcount){
00547         if(++iteration > MAX_ITER_ALLOWED) {
00548           throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00549         }
00550         CoordinatePoint <T> p = this->getPoint();
00551 
00552         bool isInHole = false;
00553         for(lno_t i = 0; i < holeCount; ++i){
00554           if(holes[i][0].isInArea(p)){
00555             isInHole = true;
00556             break;
00557           }
00558         }
00559         if(isInHole) continue;
00560         iteration = 0;
00561         points[cnt].x = p.x;
00562         points[cnt].y = p.y;
00563         points[cnt].z = p.z;
00564         ++cnt;
00565       }
00566     }
00567     */
00568   }
00569 
00570   void GetPoints(lno_t requestedPointcount, T **coords/*preallocated sized numPoints*/, lno_t tindex,
00571       Hole<T> **holes, lno_t holeCount,
00572       float *sharedRatios_, int myRank){
00573 
00574     for (int i = 0; i < myRank; ++i){
00575       //cout << "me:" << myRank << " i:" << i << " s:" << sharedRatios_[i]<< endl;
00576       this->assignedPrevious += lno_t(sharedRatios_[i] * this->numPoints);
00577       if (i < this->numPoints % this->worldSize){
00578         this->assignedPrevious += 1;
00579       }
00580     }
00581 
00582     this->requested = requestedPointcount;
00583 
00584     unsigned int slice =  UINT_MAX/(this->worldSize);
00585     unsigned int stateBegin = myRank * slice;
00586 #ifdef HAVE_ZOLTAN2_OMP
00587 #pragma omp parallel
00588 #endif
00589     {
00590       int me = 0;
00591       int tsize = 1;
00592 #ifdef HAVE_ZOLTAN2_OMP
00593       me = omp_get_thread_num();
00594       tsize = omp_get_num_threads();
00595 #endif
00596       unsigned int state = stateBegin + me * (slice/(tsize));
00597       /*
00598 #pragma omp critical
00599       {
00600 
00601         cout << "myRank:" << me << " stateBeg:" << stateBegin << " tsize:" << tsize << " state:" << state <<  " slice: " << slice / tsize <<  endl;
00602       }
00603       */
00604 #ifdef HAVE_ZOLTAN2_OMP
00605 #pragma omp for
00606 #endif
00607       for(lno_t cnt = 0; cnt < requestedPointcount; ++cnt){
00608         lno_t iteration = 0;
00609         while(1){
00610           if(++iteration > MAX_ITER_ALLOWED) {
00611             throw "Max number of Iteration is reached for point creation. Check the area criteria or hole coordinates.";
00612           }
00613           CoordinatePoint <T> p = this->getPoint( this->assignedPrevious + cnt, state);
00614 
00615           bool isInHole = false;
00616           for(lno_t i = 0; i < holeCount; ++i){
00617             if(holes[i][0].isInArea(p)){
00618               isInHole = true;
00619               break;
00620             }
00621           }
00622           if(isInHole) continue;
00623           coords[0][cnt + tindex] = p.x;
00624           if(this->dimension > 1){
00625             coords[1][cnt + tindex] = p.y;
00626             if(this->dimension > 2){
00627               coords[2][cnt + tindex] = p.z;
00628             }
00629           }
00630           break;
00631         }
00632       }
00633     }
00634   }
00635 };
00636 
00637 template <typename T, typename lno_t, typename gno_t>
00638 class CoordinateNormalDistribution:public CoordinateDistribution<T,lno_t,gno_t>{
00639 public:
00640   CoordinatePoint<T> center;
00641   T standartDevx;
00642   T standartDevy;
00643   T standartDevz;
00644 
00645 
00646   virtual T getXCenter(){
00647     return center.x;
00648   }
00649   virtual T getXRadius(){
00650     return standartDevx;
00651   }
00652 
00653   CoordinateNormalDistribution(gno_t np_, int dim, CoordinatePoint<T> center_ ,
00654                                T sd_x, T sd_y, T sd_z, int wSize) :
00655     CoordinateDistribution<T,lno_t,gno_t>(np_,dim,wSize),
00656     standartDevx(sd_x), standartDevy(sd_y), standartDevz(sd_z)
00657   {
00658     this->center.x = center_.x;
00659     this->center.y = center_.y;
00660     this->center.z = center_.z;
00661   }
00662 
00663   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00664 
00665     pindex = 0; // not used in normal distribution.
00666     CoordinatePoint <T> p;
00667 
00668     for(int i = 0; i < this->dimension; ++i){
00669       switch(i){
00670       case 0:
00671         p.x = normalDist(this->center.x, this->standartDevx, state);
00672         break;
00673       case 1:
00674         p.y = normalDist(this->center.y, this->standartDevy, state);
00675         break;
00676       case 2:
00677         p.z = normalDist(this->center.z, this->standartDevz, state);
00678         break;
00679       default:
00680         throw "unsupported dimension";
00681       }
00682     }
00683     return p;
00684   }
00685 
00686   virtual ~CoordinateNormalDistribution(){};
00687 private:
00688   T normalDist(T center_, T sd, unsigned int &state) {
00689     static bool derived=false;
00690     static T storedDerivation;
00691     T polarsqrt, normalsquared, normal1, normal2;
00692     if (!derived) {
00693       do {
00694         normal1=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
00695         normal2=2.0*( T(rand_r(&state))/T(RAND_MAX) ) - 1.0;
00696         normalsquared=normal1*normal1+normal2*normal2;
00697       } while ( normalsquared>=1.0 || normalsquared == 0.0);
00698 
00699       polarsqrt=sqrt(-2.0*log(normalsquared)/normalsquared);
00700       storedDerivation=normal1*polarsqrt;
00701       derived=true;
00702       return normal2*polarsqrt*sd + center_;
00703     }
00704     else {
00705       derived=false;
00706       return storedDerivation*sd + center_;
00707     }
00708   }
00709 };
00710 
00711 template <typename T, typename lno_t, typename gno_t>
00712 class CoordinateUniformDistribution:public CoordinateDistribution<T,lno_t, gno_t>{
00713 public:
00714   T leftMostx;
00715   T rightMostx;
00716   T leftMosty;
00717   T rightMosty;
00718   T leftMostz;
00719   T rightMostz;
00720 
00721 
00722   virtual T getXCenter(){
00723     return (rightMostx - leftMostx)/2  + leftMostx;
00724   }
00725   virtual T getXRadius(){
00726     return (rightMostx - leftMostx)/2;
00727   }
00728 
00729 
00730   CoordinateUniformDistribution(gno_t np_, int dim, T l_x, T r_x, T l_y, T r_y,
00731                                 T l_z, T r_z, int wSize ) :
00732       CoordinateDistribution<T,lno_t,gno_t>(np_,dim,wSize),
00733       leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y),
00734       leftMostz(l_z), rightMostz(r_z){}
00735 
00736   virtual ~CoordinateUniformDistribution(){};
00737   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00738 
00739 
00740     pindex = 0; //not used in uniform dist.
00741     CoordinatePoint <T> p;
00742     for(int i = 0; i < this->dimension; ++i){
00743       switch(i){
00744       case 0:
00745         p.x = uniformDist(this->leftMostx, this->rightMostx, state);
00746         break;
00747       case 1:
00748         p.y = uniformDist(this->leftMosty, this->rightMosty, state);
00749         break;
00750       case 2:
00751         p.z = uniformDist(this->leftMostz, this->rightMostz, state);
00752         break;
00753       default:
00754         throw "unsupported dimension";
00755       }
00756     }
00757     return p;
00758   }
00759 
00760 private:
00761 
00762   T uniformDist(T a, T b, unsigned int &state)
00763   {
00764     return (b-a)*(rand_r(&state) / double(RAND_MAX)) + a;
00765   }
00766 };
00767 
00768 template <typename T, typename lno_t, typename gno_t>
00769 class CoordinateGridDistribution:public CoordinateDistribution<T,lno_t,gno_t>{
00770 public:
00771   T leftMostx;
00772   T rightMostx;
00773   T leftMosty;
00774   T rightMosty;
00775   T leftMostz;
00776   T rightMostz;
00777   gno_t along_X, along_Y, along_Z;
00778   //T currentX, currentY, currentZ;
00779   T processCnt;
00780   int myRank;
00781   T xstep, ystep, zstep;
00782   gno_t xshift, yshift, zshift;
00783 
00784   virtual T getXCenter(){
00785     return (rightMostx - leftMostx)/2  + leftMostx;
00786   }
00787   virtual T getXRadius(){
00788     return (rightMostx - leftMostx)/2;
00789   }
00790 
00791 
00792   CoordinateGridDistribution(gno_t alongX, gno_t alongY, gno_t alongZ, int dim,
00793                              T l_x, T r_x, T l_y, T r_y, T l_z, T r_z ,
00794                              int myRank_, int wSize) :
00795       CoordinateDistribution<T,lno_t,gno_t>(alongX * alongY * alongZ,dim,wSize),
00796       leftMostx(l_x), rightMostx(r_x), leftMosty(l_y), rightMosty(r_y), leftMostz(l_z), rightMostz(r_z), myRank(myRank_){
00797     //currentX = leftMostx, currentY = leftMosty, currentZ = leftMostz;
00798     this->processCnt = 0;
00799     this->along_X = alongX; this->along_Y = alongY; this->along_Z = alongZ;
00800 
00801     if(this->along_X > 1)
00802       this->xstep = (rightMostx - leftMostx) / (alongX - 1);
00803     else
00804       this->xstep = 1;
00805     if(this->along_Y > 1)
00806       this->ystep = (rightMosty - leftMosty) / (alongY - 1);
00807     else
00808       this->ystep = 1;
00809     if(this->along_Z > 1)
00810       this->zstep = (rightMostz - leftMostz) / (alongZ - 1);
00811     else
00812       this->zstep = 1;
00813     xshift = 0; yshift = 0; zshift = 0;
00814   }
00815 
00816   virtual ~CoordinateGridDistribution(){};
00817   virtual CoordinatePoint<T> getPoint(gno_t pindex, unsigned int &state){
00818     //lno_t before = processCnt + this->assignedPrevious;
00819     //cout << "before:" << processCnt << " " << this->assignedPrevious << endl;
00820     //lno_t xshift = 0, yshift = 0, zshift = 0;
00821 
00822     //lno_t tmp = before % (this->along_X * this->along_Y);
00823     //xshift  = tmp % this->along_X;
00824     //yshift = tmp / this->along_X;
00825     //zshift = before / (this->along_X * this->along_Y);
00826 
00827     state = 0; //not used here
00828     this->zshift = pindex / (along_X * along_Y);
00829     this->yshift = (pindex % (along_X * along_Y)) / along_X;
00830     this->xshift = (pindex % (along_X * along_Y)) % along_X;
00831 
00832     //this->xshift = pindex / (along_Z * along_Y);
00833    // this->zshift = (pindex % (along_Z * along_Y)) / along_Y;
00834     //this->yshift = (pindex % (along_Z * along_Y)) % along_Y;
00835 
00836     CoordinatePoint <T> p;
00837     p.x = xshift * this->xstep + leftMostx;
00838     p.y = yshift * this->ystep + leftMosty;
00839     p.z = zshift * this->zstep + leftMostz;
00840 /*
00841     ++xshift;
00842     if(xshift == this->along_X){
00843       ++yshift;
00844       xshift = 0;
00845       if(yshift == this->along_Y){
00846         ++zshift;
00847         yshift = 0;
00848       }
00849     }
00850 */
00851     /*
00852     if(this->processCnt == 0){
00853       this->xshift = this->assignedPrevious / (along_Z * along_Y);
00854       //this->yshift = (this->assignedPrevious % (along_X * along_Y)) / along_X;
00855       this->zshift = (this->assignedPrevious % (along_Z * along_Y)) / along_Y;
00856       //this->xshift = (this->assignedPrevious % (along_X * along_Y)) % along_X;
00857       this->yshift = (this->assignedPrevious % (along_Z * along_Y)) % along_Y;
00858       ++this->processCnt;
00859     }
00860 
00861     CoordinatePoint <T> p;
00862     p.x = xshift * this->xstep + leftMostx;
00863     p.y = yshift * this->ystep + leftMosty;
00864     p.z = zshift * this->zstep + leftMostz;
00865 
00866     ++yshift;
00867     if(yshift == this->along_Y){
00868       ++zshift;
00869       yshift = 0;
00870       if(zshift == this->along_Z){
00871         ++xshift;
00872         zshift = 0;
00873       }
00874 
00875     }
00876     */
00877     /*
00878     if(this->requested - 1 > this->processCnt){
00879       this->processCnt++;
00880     } else {
00881       cout << "req:" << this->requested << " pr:" <<this->processCnt << endl;
00882       cout << "p:" << p.x << " " << p.y << " " << p.z <<  endl ;
00883       cout << "s:" << xshift << " " << yshift << " " << zshift <<  endl ;
00884       cout << "st:" << this->xstep << " " << this->ystep << " " << this->zstep <<  endl ;
00885     }
00886     */
00887     return p;
00888   }
00889 
00890 private:
00891 
00892 };
00893 
00894 template <typename T, typename lno_t, typename gno_t, typename node_t>
00895 class GeometricGenerator {
00896 private:
00897   Hole<T> **holes; //to represent if there is any hole in the input
00898   int holeCount;
00899   int coordinate_dimension;  //dimension of the geometry
00900   gno_t numGlobalCoords;  //global number of coordinates requested to be created.
00901   lno_t numLocalCoords;
00902   float *loadDistributions; //sized as the number of processors, the load of each processor.
00903   bool loadDistSet;
00904   bool distinctCoordSet;
00905   CoordinateDistribution<T, lno_t,gno_t> **coordinateDistributions;
00906   int distributionCount;
00907   //CoordinatePoint<T> *points;
00908   T **coords;
00909   T **wghts;
00910   WeightDistribution<T,T> **wd;
00911   int numWeightsPerCoord;
00912   int predistribution;
00913   RCP<const Teuchos::Comm<int> > comm;
00914   //RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmVector;
00915   //RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmwVector;
00916   int worldSize;
00917   int myRank;
00918   T minx;
00919   T maxx;
00920   T miny;
00921   T maxy;
00922   T minz;
00923   T maxz;
00924   std::string outfile;
00925   float perturbation_ratio;
00926 
00927   typedef Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> tMVector_t;
00928   typedef Tpetra::Map<lno_t, gno_t, node_t> tMap_t;
00929 
00930 
00931   template <typename tt>
00932   tt getParamVal( const Teuchos::ParameterEntry& pe, const std::string &paramname){
00933     tt returnVal;
00934     try {
00935       returnVal = Teuchos::getValue<tt>(pe);
00936     }
00937     catch (...){
00938       throw INVALID(paramname);
00939     }
00940     return returnVal;
00941   }
00942 
00943   int countChar (std::string inStr, char cntChar){
00944     int cnt = 0;
00945     for (unsigned int i = 0; i < inStr.size(); ++i){
00946       if (inStr[i] == cntChar) {
00947         cnt++;
00948       }
00949     }
00950     return cnt;
00951   }
00952 
00953   template <typename tt>
00954   std::string toString(tt obj){
00955     std::stringstream ss (std::stringstream::in |std::stringstream::out);
00956     ss << obj;
00957     std::string tmp = "";
00958     ss >> tmp;
00959     return tmp;
00960   }
00961 
00962   template <typename tt>
00963   tt fromString(std::string obj){
00964     std::stringstream ss (std::stringstream::in | std::stringstream::out);
00965     ss << obj;
00966     tt tmp;
00967     ss >> tmp;
00968     if (ss.fail()){
00969       throw "Cannot convert string " + obj;
00970     }
00971     return tmp;
00972   }
00973 
00974 
00975   void splitString(std::string inStr, char splitChar, std::string *outSplittedStr){
00976     std::stringstream ss (std::stringstream::in | std::stringstream::out);
00977     ss << inStr;
00978     int cnt = 0;
00979     while (!ss.eof()){
00980       std::string tmp = "";
00981       std::getline(ss, tmp ,splitChar);
00982       outSplittedStr[cnt++] = tmp;
00983     }
00984   }
00985 
00986   /*
00987   void getDistinctCoordinateDescription(std::string distinctDescription){
00988 
00989     this->distinctCoordinates = new bool[this->dimension];
00990     if(distinctDescription != ""){
00991       int argCnt = this->countChar(distinctDescription, ',') + 1;
00992       if(argCnt != this->dimension) {
00993         throw "Invalid parameter count for distinct_coordinates. Size should be equal to dimension.";
00994       }
00995 
00996       std::string *splittedStr = new std::string[argCnt];
00997       splitString(holeDescription, ',', splittedStr);
00998 
00999       for(int i = 0; i < argCnt; ++i){
01000         if(splittedStr[i] == "T"){
01001           distinctCoordinates[i] = true;
01002         } else if(splittedStr[i] == "F"){
01003           distinctCoordinates[i] = false;
01004         } else {
01005           throw "Invalid parameter " + splittedStr[i] + " for distinct_coordinates.";
01006         }
01007       }
01008       delete []splittedStr;
01009     }
01010     else {
01011       for(int i = 0; i < this->dimension; ++i){
01012         distinctCoordinates[i] = false;
01013       }
01014     }
01015   }
01016    */
01017 
01018 
01019   void getCoordinateDistributions(std::string coordinate_distributions){
01020     if(coordinate_distributions == ""){
01021       throw "At least one distribution function must be provided to coordinate generator.";
01022     }
01023 
01024     int argCnt = this->countChar(coordinate_distributions, ',') + 1;
01025     std::string *splittedStr = new std::string[argCnt];
01026     splitString(coordinate_distributions, ',', splittedStr);
01027     coordinateDistributions = (CoordinateDistribution<T, lno_t,gno_t> **) malloc(sizeof (CoordinateDistribution<T, lno_t,gno_t> *) * 1);
01028     for(int i = 0; i < argCnt; ){
01029       coordinateDistributions = (CoordinateDistribution<T, lno_t,gno_t> **)realloc((void *)coordinateDistributions, (this->distributionCount + 1)* sizeof(CoordinateDistribution<T, lno_t,gno_t> *));
01030 
01031       std::string distName = splittedStr[i++];
01032       gno_t np_ = 0;
01033       if(distName == "NORMAL"){
01034         int reqArg = 5;
01035         if (this->coordinate_dimension == 3){
01036           reqArg = 7;
01037         }
01038         if(i + reqArg > argCnt) {
01039           std::string tmp = toString<int>(reqArg);
01040           throw INVALID_SHAPE_ARG(distName, tmp);
01041         }
01042         np_ = fromString<gno_t>(splittedStr[i++]);
01043         CoordinatePoint<T> pp;
01044 
01045         pp.x = fromString<T>(splittedStr[i++]);
01046         pp.y = fromString<T>(splittedStr[i++]);
01047         pp.z = 0;
01048         if(this->coordinate_dimension == 3){
01049           pp.z = fromString<T>(splittedStr[i++]);
01050         }
01051 
01052         T sd_x = fromString<T>(splittedStr[i++]);
01053         T sd_y = fromString<T>(splittedStr[i++]);
01054         T sd_z = 0;
01055         if(this->coordinate_dimension == 3){
01056           sd_z = fromString<T>(splittedStr[i++]);
01057         }
01058         this->coordinateDistributions[this->distributionCount++] = new CoordinateNormalDistribution<T, lno_t,gno_t>(np_, this->coordinate_dimension, pp , sd_x, sd_y, sd_z, this->worldSize );
01059 
01060       } else if(distName == "UNIFORM" ){
01061         int reqArg = 5;
01062         if (this->coordinate_dimension == 3){
01063           reqArg = 7;
01064         }
01065         if(i + reqArg > argCnt) {
01066           std::string tmp = toString<int>(reqArg);
01067           throw INVALID_SHAPE_ARG(distName, tmp);
01068         }
01069         np_ = fromString<gno_t>(splittedStr[i++]);
01070         T l_x = fromString<T>(splittedStr[i++]);
01071         T r_x = fromString<T>(splittedStr[i++]);
01072         T l_y = fromString<T>(splittedStr[i++]);
01073         T r_y = fromString<T>(splittedStr[i++]);
01074 
01075         T l_z = 0, r_z = 0;
01076 
01077         if(this->coordinate_dimension == 3){
01078           l_z = fromString<T>(splittedStr[i++]);
01079           r_z = fromString<T>(splittedStr[i++]);
01080         }
01081 
01082         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 );
01083       } else if (distName == "GRID"){
01084         int reqArg = 6;
01085         if(this->coordinate_dimension == 3){
01086           reqArg = 9;
01087         }
01088         if(i + reqArg > argCnt) {
01089           std::string tmp = toString<int>(reqArg);
01090           throw INVALID_SHAPE_ARG(distName, tmp);
01091         }
01092 
01093         gno_t np_x = fromString<gno_t>(splittedStr[i++]);
01094         gno_t np_y = fromString<gno_t>(splittedStr[i++]);
01095         gno_t np_z = 1;
01096 
01097 
01098         if(this->coordinate_dimension == 3){
01099           np_z = fromString<gno_t>(splittedStr[i++]);
01100         }
01101 
01102         np_ = np_x * np_y * np_z;
01103         T l_x = fromString<T>(splittedStr[i++]);
01104         T r_x = fromString<T>(splittedStr[i++]);
01105         T l_y = fromString<T>(splittedStr[i++]);
01106         T r_y = fromString<T>(splittedStr[i++]);
01107 
01108         T l_z = 0, r_z = 0;
01109 
01110         if(this->coordinate_dimension == 3){
01111           l_z = fromString<T>(splittedStr[i++]);
01112           r_z = fromString<T>(splittedStr[i++]);
01113         }
01114 
01115         if(np_x < 1 || np_z < 1 || np_y < 1 ){
01116           throw "Provide at least 1 point along each dimension for grid test.";
01117         }
01118         //cout << "ly:" << l_y << " ry:" << r_y << endl;
01119         this->coordinateDistributions[this->distributionCount++] = new CoordinateGridDistribution<T, lno_t,gno_t>
01120         (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);
01121 
01122       }
01123       else {
01124         std::string tmp = toString<int>(this->coordinate_dimension);
01125         throw INVALIDSHAPE(distName, tmp);
01126       }
01127       this->numGlobalCoords += (gno_t) np_;
01128     }
01129     delete []splittedStr;
01130   }
01131 
01132   void getProcLoadDistributions(std::string proc_load_distributions){
01133 
01134     this->loadDistributions = new float[this->worldSize];
01135     if(proc_load_distributions == ""){
01136       float r = 1.0 / this->worldSize;
01137       for(int i = 0; i < this->worldSize; ++i){
01138         this->loadDistributions[i] = r;
01139       }
01140     } else{
01141 
01142 
01143       int argCnt = this->countChar(proc_load_distributions, ',') + 1;
01144       if(argCnt != this->worldSize) {
01145         throw "Invalid parameter count load distributions. Given " + toString<int>(argCnt) + " processor size is " + toString<int>(this->worldSize);
01146       }
01147       std::string *splittedStr = new std::string[argCnt];
01148       splitString(proc_load_distributions, ',', splittedStr);
01149       for(int i = 0; i < argCnt; ++i){
01150         this->loadDistributions[i] = fromString<float>(splittedStr[i]);
01151       }
01152       delete []splittedStr;
01153 
01154 
01155       float sum = 0;
01156       for(int i = 0; i < this->worldSize; ++i){
01157         sum += this->loadDistributions[i];
01158       }
01159       if (fabs(sum - 1.0) > 10*std::numeric_limits<float>::epsilon()){
01160         throw "Processor load ratios do not sum to 1.0.";
01161       }
01162     }
01163 
01164   }
01165 
01166   void getHoles(std::string holeDescription){
01167     if(holeDescription == ""){
01168       return;
01169     }
01170     this->holes =  (Hole<T> **) malloc(sizeof (Hole <T>*));
01171     int argCnt = this->countChar(holeDescription, ',') + 1;
01172     std::string *splittedStr = new std::string[argCnt];
01173     splitString(holeDescription, ',', splittedStr);
01174 
01175     for(int i = 0; i < argCnt; ){
01176       holes = (Hole<T> **)realloc((void *)holes, (this->holeCount + 1)* sizeof(Hole<T> *));
01177 
01178       std::string shapeName = splittedStr[i++];
01179       if(shapeName == "SQUARE" && this->coordinate_dimension == 2){
01180         if(i + 3 > argCnt) {
01181           throw INVALID_SHAPE_ARG(shapeName, "3");
01182         }
01183         CoordinatePoint<T> pp;
01184         pp.x = fromString<T>(splittedStr[i++]);
01185         pp.y = fromString<T>(splittedStr[i++]);
01186         T edge = fromString<T>(splittedStr[i++]);
01187         this->holes[this->holeCount++] = new SquareHole<T>(pp, edge);
01188       } else if(shapeName == "RECTANGLE" && this->coordinate_dimension == 2){
01189         if(i + 4 > argCnt) {
01190           throw INVALID_SHAPE_ARG(shapeName, "4");
01191         }
01192         CoordinatePoint<T> pp;
01193         pp.x = fromString<T>(splittedStr[i++]);
01194         pp.y = fromString<T>(splittedStr[i++]);
01195         T edgex = fromString<T>(splittedStr[i++]);
01196         T edgey = fromString<T>(splittedStr[i++]);
01197 
01198         this->holes[this->holeCount++] = new RectangleHole<T>(pp, edgex,edgey);
01199       } else if(shapeName == "CIRCLE" && this->coordinate_dimension == 2){
01200         if(i + 3 > argCnt) {
01201           throw INVALID_SHAPE_ARG(shapeName, "3");
01202         }
01203         CoordinatePoint<T> pp;
01204         pp.x = fromString<T>(splittedStr[i++]);
01205         pp.y = fromString<T>(splittedStr[i++]);
01206         T r = fromString<T>(splittedStr[i++]);
01207         this->holes[this->holeCount++] = new CircleHole<T>(pp, r);
01208       }  else if(shapeName == "CUBE" && this->coordinate_dimension == 3){
01209         if(i + 4 > argCnt) {
01210           throw INVALID_SHAPE_ARG(shapeName, "4");
01211         }
01212         CoordinatePoint<T> pp;
01213         pp.x = fromString<T>(splittedStr[i++]);
01214         pp.y = fromString<T>(splittedStr[i++]);
01215         pp.z = fromString<T>(splittedStr[i++]);
01216         T edge = fromString<T>(splittedStr[i++]);
01217         this->holes[this->holeCount++] = new CubeHole<T>(pp, edge);
01218       }  else if(shapeName == "RECTANGULAR_PRISM" && this->coordinate_dimension == 3){
01219         if(i + 6 > argCnt) {
01220           throw INVALID_SHAPE_ARG(shapeName, "6");
01221         }
01222         CoordinatePoint<T> pp;
01223         pp.x = fromString<T>(splittedStr[i++]);
01224         pp.y = fromString<T>(splittedStr[i++]);
01225         pp.z = fromString<T>(splittedStr[i++]);
01226         T edgex = fromString<T>(splittedStr[i++]);
01227         T edgey = fromString<T>(splittedStr[i++]);
01228         T edgez = fromString<T>(splittedStr[i++]);
01229         this->holes[this->holeCount++] = new RectangularPrismHole<T>(pp, edgex, edgey, edgez);
01230 
01231       }  else if(shapeName == "SPHERE" && this->coordinate_dimension == 3){
01232         if(i + 4 > argCnt) {
01233           throw INVALID_SHAPE_ARG(shapeName, "4");
01234         }
01235         CoordinatePoint<T> pp;
01236         pp.x = fromString<T>(splittedStr[i++]);
01237         pp.y = fromString<T>(splittedStr[i++]);
01238         pp.z = fromString<T>(splittedStr[i++]);
01239         T r = fromString<T>(splittedStr[i++]);
01240         this->holes[this->holeCount++] = new SphereHole<T>(pp, r);
01241       }  else {
01242         std::string tmp = toString<int>(this->coordinate_dimension);
01243         throw INVALIDSHAPE(shapeName, tmp);
01244       }
01245     }
01246     delete [] splittedStr;
01247   }
01248 
01249   void getWeightDistribution(std::string *weight_distribution_arr, int wdimension){
01250     int wcount = 0;
01251 
01252     this->wd = new WeightDistribution<T,T> *[wdimension];
01253     for(int ii = 0; ii < MAX_WEIGHT_DIM; ++ii){
01254       std::string weight_distribution = weight_distribution_arr[ii];
01255       if(weight_distribution == "") continue;
01256       if(wcount == wdimension) {
01257         throw "Weight Dimension is provided as " + toString<int>(wdimension) + ". More weight distribution is provided.";
01258       }
01259 
01260       int count = this->countChar(weight_distribution, ' ');
01261       std::string *splittedStr = new string[count + 1];
01262       this->splitString(weight_distribution, ' ', splittedStr);
01263       //cout << count << endl;
01264       T c=1;
01265       T a1=0,a2=0,a3=0;
01266       T x1=0,y1=0,z1=0;
01267       T b1=1,b2=1,b3=1;
01268       T *steps = NULL;
01269       T *values= NULL;
01270       int stepCount = 0;
01271       int valueCount = 1;
01272 
01273       for (int i = 1; i < count + 1; ++i){
01274         int assignmentCount = this->countChar(splittedStr[i], '=');
01275         if (assignmentCount != 1){
01276           throw "Error at the argument " + splittedStr[i];
01277         }
01278         std::string parameter_vs_val[2];
01279         this->splitString(splittedStr[i], '=', parameter_vs_val);
01280         std::string parameter = parameter_vs_val[0];
01281         std::string value = parameter_vs_val[1];
01282         //cout << "parameter:" << parameter << " value:" << value << endl;
01283 
01284         if (parameter == "a1"){
01285           a1 = this->fromString<T>(value);
01286         }
01287         else if (parameter == "a2"){
01288           if(this->coordinate_dimension > 1){
01289             a2 = this->fromString<T>(value);
01290           }
01291           else {
01292             throw  parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01293           }
01294 
01295         }
01296         else if (parameter == "a3"){
01297           if(this->coordinate_dimension > 2){
01298             a3 = this->fromString<T>(value);
01299           }
01300           else {
01301             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01302           }
01303         }
01304         else if (parameter == "b1"){
01305           b1 = this->fromString<T>(value);
01306         }
01307         else if (parameter == "b2"){
01308           if(this->coordinate_dimension > 1){
01309             b2 = this->fromString<T>(value);
01310           }
01311           else {
01312             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01313           }
01314         }
01315         else if (parameter == "b3"){
01316 
01317           if(this->coordinate_dimension > 2){
01318             b3 = this->fromString<T>(value);
01319           }
01320           else {
01321             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01322           }
01323         }
01324         else if (parameter == "c"){
01325           c = this->fromString<T>(value);
01326         }
01327         else if (parameter == "x1"){
01328           x1 = this->fromString<T>(value);
01329         }
01330         else if (parameter == "y1"){
01331           if(this->coordinate_dimension > 1){
01332             y1 = this->fromString<T>(value);
01333           }
01334           else {
01335             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01336           }
01337         }
01338         else if (parameter == "z1"){
01339           if(this->coordinate_dimension > 2){
01340             z1 = this->fromString<T>(value);
01341           }
01342           else {
01343             throw parameter+ " argument is not valid when dimension is " + toString<int>(this->coordinate_dimension);
01344           }
01345         }
01346         else if (parameter == "steps"){
01347           stepCount = this->countChar(value, ',') + 1;
01348           std::string *stepstr = new std::string[stepCount];
01349           this->splitString(value, ',', stepstr);
01350           steps = new T[stepCount];
01351           for (int j = 0; j < stepCount; ++j){
01352             steps[j] = this->fromString<T>(stepstr[j]);
01353           }
01354           delete [] stepstr;
01355         }
01356         else if (parameter == "values"){
01357           valueCount = this->countChar(value, ',') + 1;
01358           std::string *stepstr = new std::string[valueCount];
01359           this->splitString(value, ',', stepstr);
01360           values = new T[valueCount];
01361           for (int j = 0; j < valueCount; ++j){
01362             values[j] = this->fromString<T>(stepstr[j]);
01363           }
01364           delete [] stepstr;
01365         }
01366         else {
01367           throw "Invalid parameter name at " + splittedStr[i];
01368         }
01369       }
01370 
01371       delete []splittedStr;
01372       if(stepCount + 1!= valueCount){
01373         throw "Step count: " + this->toString<int>(stepCount) + " must be 1 less than value count: " + this->toString<int>(valueCount);
01374       }
01375 
01376 
01377       this->wd[wcount++] =  new SteppedEquation<T,T>(a1, a2,  a3,  b1,  b2,  b3,  c,  x1,  y1,  z1, steps, values, stepCount);
01378 
01379       if(stepCount > 0){
01380         delete [] steps;
01381         delete [] values;
01382 
01383       }
01384     }
01385     if(wcount != this->numWeightsPerCoord){
01386       throw "Weight Dimension is provided as " + toString<int>(wdimension) + ". But " + toString<int>(wcount)+" weight distributions are provided.";
01387     }
01388   }
01389 
01390   void parseParams(const Teuchos::ParameterList & params){
01391     try {
01392       std::string holeDescription  = "";
01393       std::string proc_load_distributions = "";
01394       std::string distinctDescription = "";
01395       std::string coordinate_distributions = "";
01396       std::string numWeightsPerCoord_parameters[MAX_WEIGHT_DIM];
01397       for (int i = 0; i < MAX_WEIGHT_DIM; ++i){
01398         numWeightsPerCoord_parameters[i] = "";
01399       }
01400 
01401 
01402       for (Teuchos::ParameterList::ConstIterator pit = params.begin(); pit != params.end(); ++pit ){
01403         const std::string &paramName = params.name(pit);
01404 
01405         const Teuchos::ParameterEntry &pe = params.entry(pit);
01406 
01407         if(paramName.find("hole-") == 0){
01408           if(holeDescription == ""){
01409             holeDescription = getParamVal<std::string>(pe, paramName);
01410           }
01411           else {
01412             holeDescription +=","+getParamVal<std::string>(pe, paramName);
01413           }
01414         }
01415         else if(paramName.find("distribution-") == 0){
01416           if(coordinate_distributions == "")
01417             coordinate_distributions = getParamVal<std::string>(pe, paramName);
01418           else
01419             coordinate_distributions +=  ","+getParamVal<std::string>(pe, paramName);
01420           //cout << coordinate_distributions << endl;
01421           //TODO coordinate distribution description
01422         }
01423 
01424         else if (paramName.find(weight_distribution_string) == 0){
01425           std::string weight_dist_param = paramName.substr(weight_distribution_string.size());
01426           int dash_pos = weight_dist_param.find("-");
01427           std::string distribution_index_string = weight_dist_param.substr(0, dash_pos);
01428           int distribution_index = fromString<int>(distribution_index_string);
01429 
01430           if(distribution_index >= MAX_WEIGHT_DIM){
01431             throw "Given distribution index:" + distribution_index_string + " larger than maximum allowed number of weights:" + toString<int>(MAX_WEIGHT_DIM);
01432           }
01433           numWeightsPerCoord_parameters[distribution_index] +=  " " + weight_dist_param.substr(dash_pos + 1)+ "="+ getParamVal<std::string>(pe, paramName);
01434         }
01435         else if(paramName == "dim"){
01436           int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
01437           if(dim < 2 && dim > 3){
01438             throw INVALID(paramName);
01439           } else {
01440             this->coordinate_dimension = dim;
01441           }
01442         }
01443         else if(paramName == "wdim"){
01444           int dim = fromString<int>(getParamVal<std::string>(pe, paramName));
01445           if(dim < 1 && dim > MAX_WEIGHT_DIM){
01446             throw INVALID(paramName);
01447           } else {
01448             this->numWeightsPerCoord = dim;
01449           }
01450         }
01451         else if(paramName == "predistribution"){
01452           int pre_distribution = fromString<int>(getParamVal<std::string>(pe, paramName));
01453           if(pre_distribution < 0 && pre_distribution > 3){
01454             throw INVALID(paramName);
01455           } else {
01456             this->predistribution = pre_distribution;
01457           }
01458         }
01459         else if(paramName == "perturbation_ratio"){
01460           float perturbation = fromString<float>(getParamVal<std::string>(pe, paramName));
01461           if(perturbation < 0 && perturbation > 1){
01462             throw INVALID(paramName);
01463           } else {
01464             this->perturbation_ratio = perturbation;
01465           }
01466         }
01467 
01468 
01469         else if(paramName == "proc_load_distributions"){
01470           proc_load_distributions = getParamVal<std::string>(pe, paramName);
01471           this->loadDistSet = true;
01472         }
01473 
01474         else if(paramName == "distinct_coordinates"){
01475           distinctDescription = getParamVal<std::string>(pe, paramName);
01476           if(distinctDescription == "T"){
01477             this->distinctCoordSet = true;
01478           } else if(distinctDescription == "F"){
01479             this->distinctCoordSet = true;
01480           } else {
01481             throw "Invalid parameter for distinct_coordinates: " + distinctDescription + ". Candiates: T or F." ;
01482           }
01483         }
01484 
01485         else if(paramName == "out_file"){
01486           this->outfile = getParamVal<std::string>(pe, paramName);
01487         }
01488 
01489         else {
01490           throw INVALID(paramName);
01491         }
01492       }
01493 
01494 
01495       if(this->coordinate_dimension == 0){
01496         throw "Dimension must be provided to coordinate generator.";
01497       }
01498       /*
01499       if(this->globalNumCoords == 0){
01500         throw "Number of coordinates must be provided to coordinate generator.";
01501       }
01502        */
01503       /*
01504       if(maxx <= minx ){
01505         throw "Error: maxx= "+ toString<T>(maxx)+ " and minx=" + toString<T>(minx);
01506       }
01507       if(maxy <= miny ){
01508         throw "Error: maxy= "+ toString<T>(maxy)+ " and miny=" + toString<T>(miny);
01509 
01510       }
01511       if(this->dimension == 3 && maxz <= minz ){
01512         throw "Error: maxz= "+ toString<T>(maxz)+ " and minz=" + toString<T>(minz);
01513       }
01514        */
01515       if (this->loadDistSet && this->distinctCoordSet){
01516         throw "distinct_coordinates and proc_load_distributions parameters cannot be satisfied together.";
01517       }
01518       this->getHoles(holeDescription);
01519       //this->getDistinctCoordinateDescription(distinctDescription);
01520       this->getProcLoadDistributions(proc_load_distributions);
01521       this->getCoordinateDistributions(coordinate_distributions);
01522       this->getWeightDistribution(numWeightsPerCoord_parameters, this->numWeightsPerCoord);
01523       /*
01524       if(this->numGlobalCoords <= 0){
01525         throw "Must have at least 1 point";
01526       }
01527        */
01528     }
01529     catch(std::string s){
01530       throw s;
01531     }
01532 
01533     catch(char * s){
01534       throw s;
01535     }
01536 
01537     catch(char const* s){
01538       throw s;
01539     }
01540 
01541   }
01542 public:
01543 
01544   ~GeometricGenerator(){
01545     if(this->holes){
01546       for (int i = 0; i < this->holeCount; ++i){
01547         delete this->holes[i];
01548       }
01549       free (this->holes);
01550     }
01551 
01552 
01553     delete []loadDistributions; //sized as the number of processors, the load of each processor.
01554     //delete []distinctCoordinates; // if processors have different or same range for coordinates to be created.
01555     if(coordinateDistributions){
01556 
01557       for (int i = 0; i < this->distributionCount; ++i){
01558         delete this->coordinateDistributions[i];
01559       }
01560       free (this->coordinateDistributions);
01561     }
01562     if (this->wd){
01563       for (int i = 0; i < this->numWeightsPerCoord; ++i){
01564         delete this->wd[i];
01565       }
01566       delete []this->wd;
01567     }
01568 
01569     if(this->numWeightsPerCoord){
01570       for(int i = 0; i < this->numWeightsPerCoord; ++i)
01571       delete [] this->wghts[i];
01572       delete []this->wghts;
01573     }
01574     if(this->coordinate_dimension){
01575       for(int i = 0; i < this->coordinate_dimension; ++i)
01576       delete [] this->coords[i];
01577       delete [] this->coords;
01578     }
01579     //delete []this->points;
01580   }
01581 
01582   void print_description(){
01583     cout <<"\nGeometric Generator Parameter File Format:" << endl;
01584     cout <<"- dim=Coordinate Dimension: 2 or 3" << endl;
01585     cout <<"- Available distributions:" << endl;
01586     cout <<"\tUNIFORM: -> distribution-1=UNIFORM,NUMCOORDINATES,XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << endl;
01587     cout <<"\tGRID: -> distribution-2=GRID,XLENGTH,YLENGTH{,ZLENGTH},XMIN,XMAX,YMIN,YMAX{,ZMIN,ZMAX}" << endl;
01588     cout <<"\tNORMAL: -> distribution-3=NORMAL,XCENTER,YCENTER{,ZCENTER},XSD,YSD,{,ZSD}" << endl;
01589     cout <<"- wdim=numWeightsPerCoord:  There should be as many weight function as number of weights per coord." << endl;
01590     cout <<"- Weight Equation: w = (a1 * (x - x1)^b1) + (a2 * (y - y1)^b2) + (a3 * (z - z1)^b3) + c" << endl;
01591     cout << "Parameter settings:" << endl;
01592     cout << "\tWeightDistribution-1-a1=a1 " << endl;
01593     cout << "\tWeightDistribution-1-a2=a2 " << endl;
01594     cout << "\tWeightDistribution-1-a3=a3 " << endl;
01595     cout << "\tWeightDistribution-1-b1=b1 " << endl;
01596     cout << "\tWeightDistribution-1-b2=b2 " << endl;
01597     cout << "\tWeightDistribution-1-b3=b3 " << endl;
01598     cout << "\tWeightDistribution-1-x1=x1 " << endl;
01599     cout << "\tWeightDistribution-1-y1=y1 " << endl;
01600     cout << "\tWeightDistribution-1-z1=z1 " << endl;
01601     cout << "\tWeightDistribution-1-c=c " << endl;
01602     cout << "\tIt is possible to set step function to the result of weight equation." << endl;
01603     cout << "\tWeightDistribution-1-steps=step1,step2,step3:increasing order" << endl;
01604     cout << "\tWeightDistribution-1-values=val1,val2,val3,val4." << endl;
01605     cout << "\t\tIf w < step1 -> w = val1" << endl;
01606     cout << "\t\tElse if w < step2 -> w = val2" << endl;
01607     cout << "\t\tElse if w < step3 -> w = val3" << endl;
01608     cout << "\t\tElse  -> w = val4" << endl;
01609     cout <<"- Holes:" << endl;
01610     cout << "\thole-1:SPHERE,XCENTER,YCENTER,ZCENTER,RADIUS (only for dim=3)" << endl;
01611     cout << "\thole-2:CUBE,XCENTER,YCENTER,ZCENTER,EDGE (only for dim=3)" << endl;
01612     cout << "\thole-3:RECTANGULAR_PRISM,XCENTER,YCENTER,ZCENTER,XEDGE,YEDGE,ZEDGE (only for dim=3)" << endl;
01613     cout << "\thole-4:SQUARE,XCENTER,YCENTER,EDGE (only for dim=2)" << endl;
01614     cout << "\thole-5:RECTANGLE,XCENTER,YCENTER,XEDGE,YEDGE (only for dim=2)" << endl;
01615     cout << "\thole-6:CIRCLE,XCENTER,YCENTER,RADIUS (only for dim=2)" << endl;
01616     cout << "- out_file:out_file_path : if provided output will be written to files." << endl;
01617     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;
01618     cout << "- predistribution:distribution_option. Predistribution of the coordinates to the processors. 0 for NONE, 1 RCB, 2 MJ, 3 BLOCK." << endl;
01619     cout << "- perturbation_ratio:the percent of the local data which will be perturbed in order to simulate the changes in the dynamic partitioning. Float value between 0 and 1." << endl;
01620   }
01621 
01622   GeometricGenerator(Teuchos::ParameterList &params, const RCP<const Teuchos::Comm<int> > & comm_){
01623     this->wd = NULL;
01624     this->comm = comm_;
01625     this->holes = NULL; //to represent if there is any hole in the input
01626     this->coordinate_dimension = 0;  //dimension of the geometry
01627     this->numWeightsPerCoord = 0;
01628     this->worldSize = comm_->getSize(); //comminication world object.
01629     this->numGlobalCoords = 0;  //global number of coordinates requested to be created.
01630     this->loadDistributions = NULL; //sized as the number of processors, the load of each processor.
01631     //this->distinctCoordinates = NULL; // if processors have different or same range for coordinates to be created.
01632     this->coordinateDistributions = NULL;
01633     this->holeCount = 0;
01634     this->distributionCount = 0;
01635     this->outfile = "";
01636     this->predistribution = 0;
01637     this->perturbation_ratio = 0;
01638     //this->points =  NULL;
01639 
01640     /*
01641     this->minx = 0;
01642     this->maxx = 0;
01643     this->miny = 0;
01644     this->maxy = 0;
01645     this->minz = 0;
01646     this->maxz = 0;
01647      */
01648     this->loadDistSet = false;
01649     this->distinctCoordSet = false;
01650     this->myRank = comm_->getRank();
01651 
01652     try {
01653       this->parseParams(params);
01654     }
01655     catch(std::string s){
01656       if(myRank == 0){
01657         print_description();
01658       }
01659       throw s;
01660     }
01661 
01662 
01663 
01664 
01665     lno_t myPointCount = 0;
01666     this->numGlobalCoords = 0;
01667 
01668     gno_t prefixSum = 0;
01669     for(int i = 0; i < this->distributionCount; ++i){
01670       for(int ii = 0; ii < this->worldSize; ++ii){
01671         lno_t increment  = lno_t (this->coordinateDistributions[i]->numPoints * this->loadDistributions[ii]);
01672         if (ii < this->coordinateDistributions[i]->numPoints % this->worldSize){
01673           increment += 1;
01674         }
01675         this->numGlobalCoords += increment;
01676         if(ii < myRank){
01677           prefixSum += increment;
01678         }
01679       }
01680       myPointCount += lno_t(this->coordinateDistributions[i]->numPoints * this->loadDistributions[myRank]);
01681       if (myRank < this->coordinateDistributions[i]->numPoints % this->worldSize){
01682         myPointCount += 1;
01683       }
01684     }
01685 
01686     this->coords = new T *[this->coordinate_dimension];
01687     for(int i = 0; i < this->coordinate_dimension; ++i){
01688       this->coords[i] = new T[myPointCount];
01689     }
01690 
01691     for (int ii = 0; ii < this->coordinate_dimension; ++ii){
01692 #ifdef HAVE_ZOLTAN2_OMP
01693 #pragma omp parallel for
01694 #endif
01695       for(int i = 0; i < myPointCount; ++i){
01696         this->coords[ii][i] = 0;
01697       }
01698     }
01699 
01700     this->numLocalCoords = 0;
01701     srand ((myRank + 1) * this->numLocalCoords);
01702     for (int i = 0; i < distributionCount; ++i){
01703 
01704       lno_t requestedPointCount = lno_t(this->coordinateDistributions[i]->numPoints *  this->loadDistributions[myRank]);
01705       if (myRank < this->coordinateDistributions[i]->numPoints % this->worldSize){
01706         requestedPointCount += 1;
01707       }
01708       //cout << "req:" << requestedPointCount << endl;
01709       //this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->points + this->numLocalCoords, this->holes, this->holeCount,  this->loadDistributions, myRank);
01710       this->coordinateDistributions[i]->GetPoints(requestedPointCount,this->coords, this->numLocalCoords, this->holes, this->holeCount,  this->loadDistributions, myRank);
01711       this->numLocalCoords += requestedPointCount;
01712     }
01713 
01714     /*
01715 
01716     if (this->myRank >= 0){
01717         for(lno_t i = 0; i < this->numLocalCoords; ++i){
01718 
01719           cout <<"me:" << this->myRank << " "<< this->coords[0][i];
01720           if(this->coordinate_dimension > 1){
01721             cout << " " << this->coords[1][i];
01722           }
01723           if(this->coordinate_dimension > 2){
01724             cout  << " " << this->coords[2][i];
01725           }
01726           cout << std::endl;
01727         }
01728     }
01729   */
01730 
01731 
01732 
01733     if (this->predistribution){
01734       redistribute();
01735     }
01736 
01737 
01738 
01739     int scale = 3;
01740     if (this->perturbation_ratio > 0){
01741       this->perturb_data(this->perturbation_ratio, scale);
01742     }
01743     /*
01744     if (this->myRank >= 0){
01745       cout << "after distribution" << endl;
01746         for(lno_t i = 0; i < this->numLocalCoords; ++i){
01747 
01748           cout <<"me:" << this->myRank << " " << this->coords[0][i];
01749           if(this->coordinate_dimension > 1){
01750             cout << " " << this->coords[1][i];
01751           }
01752           if(this->coordinate_dimension > 2){
01753             cout  << " " << this->coords[2][i];
01754           }
01755           cout << std::endl;
01756         }
01757     }
01758 
01759     */
01760 
01761 
01762     if (this->distinctCoordSet){
01763       //TODO: Partition and migration.
01764     }
01765 
01766     if(this->outfile != ""){
01767 
01768       std::ofstream myfile;
01769       myfile.open ((this->outfile + toString<int>(myRank)).c_str());
01770       for(lno_t i = 0; i < this->numLocalCoords; ++i){
01771 
01772         myfile << this->coords[0][i];
01773         if(this->coordinate_dimension > 1){
01774           myfile << " " << this->coords[1][i];
01775         }
01776         if(this->coordinate_dimension > 2){
01777           myfile << " " << this->coords[2][i];
01778         }
01779         myfile << std::endl;
01780       }
01781       myfile.close();
01782 
01783       if (this->myRank == 0){
01784         std::ofstream gnuplotfile("gnu.gnuplot");
01785         for(int i = 0; i < this->worldSize; ++i){
01786           string s = "splot";
01787           if (this->coordinate_dimension == 2){
01788             s = "plot";
01789           }
01790           if (i > 0){
01791             s = "replot";
01792           }
01793           gnuplotfile << s << " \"" << (this->outfile + toString<int>(i)) << "\"" << endl;
01794         }
01795         gnuplotfile  << "pause -1" << endl;
01796         gnuplotfile.close();
01797       }
01798     }
01799 
01800 
01801 
01802     /*
01803     Zoltan2::XpetraMultiVectorAdapter < Tpetra::MultiVector<T, lno_t, gno_t, node_t> > xmv (RCP < Tpetra::MultiVector<T, lno_t, gno_t, node_t> > (tmVector));
01804 
01805     RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmVector2;
01806     Zoltan2::PartitioningSolution< Tpetra::MultiVector<T, lno_t, gno_t, node_t> > solution;
01807     xmv.applyPartitioningSolution<Tpetra::MultiVector<T, lno_t, gno_t, node_t> >(this->tmVector, &tmVector2, solution);
01808      */
01809     if (this->numWeightsPerCoord > 0){
01810       this->wghts = new T *[this->numWeightsPerCoord];
01811       for(int i = 0; i < this->numWeightsPerCoord; ++i){
01812         this->wghts[i] = new T[this->numLocalCoords];
01813       }
01814     }
01815 
01816     for(int ii = 0; ii < this->numWeightsPerCoord; ++ii){
01817       switch(this->coordinate_dimension){
01818       case 1:
01819 #ifdef HAVE_ZOLTAN2_OMP
01820 #pragma omp parallel for
01821 #endif
01822         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01823           this->wghts[ii][i] = this->wd[ii]->get1DWeight(this->coords[0][i]);
01824         }
01825         break;
01826       case 2:
01827 #ifdef HAVE_ZOLTAN2_OMP
01828 #pragma omp parallel for
01829 #endif
01830         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01831           this->wghts[ii][i] = this->wd[ii]->get2DWeight(this->coords[0][i], this->coords[1][i]);
01832         }
01833         break;
01834       case 3:
01835 #ifdef HAVE_ZOLTAN2_OMP
01836 #pragma omp parallel for
01837 #endif
01838         for (lno_t i = 0; i < this->numLocalCoords; ++i){
01839           this->wghts[ii][i] = this->wd[ii]->get3DWeight(this->coords[0][i], this->coords[1][i], this->coords[2][i]);
01840         }
01841         break;
01842       }
01843     }
01844   }
01845 
01846   //############################################################//
01848   //############################################################//
01849   void perturb_data(float used_perturbation_ratio, int scale){
01850     T *maxCoords= new T [this->coordinate_dimension];
01851     T *minCoords= new T [this->coordinate_dimension];
01852     for (int dim = 0; dim < this->coordinate_dimension; ++dim){
01853       minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
01854       for (lno_t i = 1; i < this->numLocalCoords; ++i){
01855         if (minCoords[dim] > this->coords[dim][i]){
01856           minCoords[dim] = this->coords[dim][i];
01857         }
01858 
01859         if (maxCoords[dim] < this->coords[dim][i]){
01860           maxCoords[dim] = this->coords[dim][i];
01861         }
01862       }
01863 
01864 
01865 
01866 
01867       T center = (maxCoords[dim] + minCoords[dim]) / 2;
01868 
01869       minCoords[dim] = center - (center - minCoords[dim]) * scale;
01870       maxCoords[dim] = (maxCoords[dim] - center) * scale + center;
01871 
01872     }
01873 
01874     gno_t numLocalToPerturb = this->numLocalCoords * used_perturbation_ratio;
01875     //cout << "numLocalToPerturb :" << numLocalToPerturb  << endl;
01876     for (int dim = 0; dim < this->coordinate_dimension; ++dim){
01877       T range = maxCoords[dim] - minCoords[dim];
01878       for (int i = 0; i < numLocalToPerturb; ++i){
01879         this->coords[dim][i] = (rand() / double (RAND_MAX)) * (range) +  minCoords[dim];
01880 
01881       }
01882     }
01883     delete []maxCoords;
01884     delete []minCoords;
01885   }
01886 
01887   //############################################################//
01889   //############################################################//
01890 
01891   //Returns the partitioning dimension as even as possible
01892   void getBestSurface (int remaining, int *dimProcs, int dim, int currentDim, int &bestSurface, int *bestDimProcs){
01893 
01894     if (currentDim < dim - 1){
01895       int ipx = 1;
01896       while(ipx <= remaining) {
01897         if(remaining % ipx == 0) {
01898           int nremain = remaining / ipx;
01899           dimProcs[currentDim] = ipx;
01900           getBestSurface (nremain, dimProcs, dim, currentDim + 1, bestSurface, bestDimProcs);
01901         }
01902         ipx++;
01903       }
01904     }
01905     else {
01906       dimProcs[currentDim] = remaining;
01907       int surface = 0;
01908       for (int i = 0; i < dim; ++i) surface += dimProcs[i];
01909       if (surface < bestSurface){
01910         bestSurface = surface;
01911         for (int i = 0; i < dim; ++i) bestDimProcs[i] = dimProcs[i];
01912       }
01913     }
01914 
01915   }
01916 
01917   //returns min and max coordinates along each dimension
01918   void getMinMaxCoords(T *globalMaxCoords, T *globalMinCoords){
01919     T *maxCoords= new T [this->coordinate_dimension];
01920     T *minCoords= new T [this->coordinate_dimension];
01921     for (int dim = 0; dim < this->coordinate_dimension; ++dim){
01922       minCoords[dim] = maxCoords[dim] = this->coords[dim][0];
01923       for (lno_t i = 1; i < this->numLocalCoords; ++i){
01924         if (minCoords[dim] > this->coords[dim][i]){
01925           minCoords[dim] = this->coords[dim][i];
01926         }
01927 
01928         if (maxCoords[dim] < this->coords[dim][i]){
01929           maxCoords[dim] = this->coords[dim][i];
01930         }
01931       }
01932     }
01933 
01934     reduceAll<int, T>( *(this->comm), Teuchos::REDUCE_MAX,
01935                             this->coordinate_dimension,
01936                             maxCoords,
01937                             globalMaxCoords);
01938 
01939 
01940     reduceAll<int, T>( *(this->comm), Teuchos::REDUCE_MIN,
01941                                 this->coordinate_dimension,
01942                                 minCoords,
01943                                 globalMinCoords);
01944 
01945     delete []maxCoords;
01946     delete []minCoords;
01947   }
01948 
01949 
01950   //performs a block partitioning.
01951   //then distributes the points of the overloaded parts to underloaded parts.
01952   void blockPartition(int *coordinate_grid_parts){
01953 
01954 
01955     //############################################################//
01957     //############################################################//
01958 
01959     T *maxCoords= new T [this->coordinate_dimension];
01960     T *minCoords= new T [this->coordinate_dimension];
01961     //global min and max coordinates in each dimension.
01962     this->getMinMaxCoords(maxCoords, minCoords);
01963 
01964 
01965     //############################################################//
01967     //############################################################//
01968     int remaining = this->worldSize;
01969     int coord_dim = this->coordinate_dimension;
01970     int *dimProcs = new int[coord_dim];
01971     int *bestDimProcs = new int[coord_dim];
01972     int currentDim = 0;
01973 
01974     int bestSurface = 0;
01975     dimProcs[0] = remaining;
01976     for (int i = 1; i < coord_dim; ++i) dimProcs[i] = 1;
01977     for (int i = 0; i < coord_dim; ++i) bestSurface += dimProcs[i];
01978     for (int i = 0; i < coord_dim; ++i) bestDimProcs[i] = dimProcs[i];
01979     //divides the parts into dimensions as even as possible.
01980     getBestSurface ( remaining, dimProcs,  coord_dim,  currentDim, bestSurface, bestDimProcs);
01981 
01982 
01983     delete []dimProcs;
01984 
01985     //############################################################//
01987     //############################################################//
01988     int *shiftProcCount = new int[coord_dim];
01989     //how the consecutive parts along a dimension
01990     //differs in the index.
01991 
01992     int remainingProc = this->worldSize;
01993     for (int dim = 0; dim < coord_dim; ++dim){
01994       remainingProc = remainingProc /  bestDimProcs[dim];
01995       shiftProcCount[dim] = remainingProc;
01996     }
01997 
01998     T *dim_slices = new T[coord_dim];
01999     for (int dim = 0; dim < coord_dim; ++dim){
02000       T dim_range = maxCoords[dim] - minCoords[dim];
02001       T slice = dim_range / bestDimProcs[dim];
02002       dim_slices[dim] = slice;
02003     }
02004 
02005     //############################################################//
02007     //############################################################//
02008 
02009     gno_t *numPointsInParts = new gno_t[this->worldSize];
02010     gno_t *numGlobalPointsInParts = new gno_t[this->worldSize];
02011     gno_t *numPointsInPartsInclusiveUptoMyIndex = new gno_t[this->worldSize];
02012 
02013     gno_t *partBegins = new gno_t [this->worldSize];
02014     gno_t *partNext = new gno_t[this->numLocalCoords];
02015 
02016 
02017     for (int i = 0; i < this->numLocalCoords; ++i){
02018       partNext[i] = -1;
02019     }
02020     for (int i = 0; i < this->worldSize; ++i) {
02021       partBegins[i] = 1;
02022     }
02023 
02024     for (int i = 0; i < this->worldSize; ++i)
02025       numPointsInParts[i] = 0;
02026 
02027     for (int i = 0; i < this->numLocalCoords; ++i){
02028       int partIndex = 0;
02029       for (int dim = 0; dim < coord_dim; ++dim){
02030         int shift = int ((this->coords[dim][i] - minCoords[dim]) / dim_slices[dim]);
02031         if (shift >= bestDimProcs[dim]){
02032           shift = bestDimProcs[dim] - 1;
02033         }
02034         shift = shift * shiftProcCount[dim];
02035         partIndex += shift;
02036       }
02037       numPointsInParts[partIndex] += 1;
02038       coordinate_grid_parts[i] = partIndex;
02039 
02040       partNext[i] = partBegins[partIndex];
02041       partBegins[partIndex] = i;
02042 
02043     }
02044 
02045     //############################################################//
02047     //############################################################//
02048     reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
02049                                 this->worldSize,
02050                                 numPointsInParts,
02051                                 numGlobalPointsInParts);
02052 
02053 
02054       Teuchos::scan<int,gno_t>(
02055           *(this->comm), Teuchos::REDUCE_SUM,
02056           this->worldSize,
02057           numPointsInParts,
02058           numPointsInPartsInclusiveUptoMyIndex
02059       );
02060 
02061 
02062 
02063 
02064     /*
02065     gno_t totalSize = 0;
02066     for (int i = 0; i < this->worldSize; ++i){
02067       totalSize += numPointsInParts[i];
02068     }
02069     if (totalSize != this->numLocalCoords){
02070       cout << "me:" << this->myRank << " ts:" << totalSize << " nl:" << this->numLocalCoords << endl;
02071     }
02072     */
02073 
02074 
02075       //cout << "me:" << this->myRank << " ilk print" << endl;
02076 
02077     gno_t optimal_num = this->numGlobalCoords / double (this->worldSize) + 0.5;
02078 #ifdef printparts
02079     if (this->myRank == 0){
02080       gno_t totalSize = 0;
02081       for (int i = 0; i < this->worldSize; ++i){
02082         cout << "me:" << this->myRank << " NumPoints in part:" << i << " is: " << numGlobalPointsInParts[i] << endl;
02083         totalSize += numGlobalPointsInParts[i];
02084       }
02085       cout << "Total:" << totalSize << " ng:" << this->numGlobalCoords << endl;
02086       cout << "optimal_num:" << optimal_num << endl;
02087     }
02088 #endif
02089     gno_t *extraInPart = new gno_t [this->worldSize];
02090 
02091     gno_t extraExchange = 0;
02092     for (int i = 0; i < this->worldSize; ++i){
02093       extraInPart[i] = numGlobalPointsInParts[i] - optimal_num;
02094       extraExchange += extraInPart[i];
02095     }
02096     if (extraExchange != 0){
02097       int addition = -1;
02098       if (extraExchange < 0) addition = 1;
02099       for (int i = 0; i < extraExchange; ++i){
02100         extraInPart[i % this->worldSize] += addition;
02101       }
02102     }
02103 
02104     //############################################################//
02106     //############################################################//
02107 
02108       int overloadedPartCount = 0;
02109       int *overloadedPartIndices = new int [this->worldSize];
02110 
02111 
02112       int underloadedPartCount = 0;
02113       int *underloadedPartIndices = new int [this->worldSize];
02114 
02115       for (int i = 0; i < this->worldSize; ++i){
02116         if(extraInPart[i] > 0){
02117           overloadedPartIndices[overloadedPartCount++] = i;
02118         }
02119         else if(extraInPart[i] < 0){
02120           underloadedPartIndices[underloadedPartCount++] = i;
02121         }
02122       }
02123 
02124       int underloadpartindex = underloadedPartCount - 1;
02125 
02126 
02127       int numPartsISendFrom = 0;
02128       int *mySendFromParts = new int[this->worldSize * 2];
02129       gno_t *mySendFromPartsCounts = new gno_t[this->worldSize * 2];
02130 
02131       int numPartsISendTo = 0;
02132       int *mySendParts = new int[this->worldSize * 2];
02133       gno_t *mySendCountsToParts = new gno_t[this->worldSize * 2];
02134 
02135 
02136     //############################################################//
02141     //############################################################//
02142       for (int i = overloadedPartCount - 1; i >= 0; --i){
02143         //get the overloaded part
02144         //the overload
02145         int overloadPart = overloadedPartIndices[i];
02146         gno_t overload = extraInPart[overloadPart];
02147         gno_t myload = numPointsInParts[overloadPart];
02148 
02149 
02150         //the inclusive load of the processors up to me
02151         gno_t inclusiveLoadUpToMe = numPointsInPartsInclusiveUptoMyIndex[overloadPart];
02152 
02153         //the exclusive load of the processors up to me
02154         gno_t exclusiveLoadUptoMe = inclusiveLoadUpToMe - myload;
02155 
02156 
02157         if (exclusiveLoadUptoMe >= overload){
02158           //this processor does not have to convert anything.
02159           //for this overloaded part.
02160           //set the extra for this processor to zero.
02161           overloadedPartIndices[i] = -1;
02162           extraInPart[overloadPart] = 0;
02163           //then consume underloaded parts.
02164           while (overload > 0){
02165             int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
02166             gno_t underload = extraInPart[nextUnderLoadedPart];
02167             gno_t left = overload + underload;
02168 
02169             if(left >= 0){
02170               extraInPart[nextUnderLoadedPart] = 0;
02171               underloadedPartIndices[underloadpartindex--] = -1;
02172             }
02173             else {
02174               extraInPart[nextUnderLoadedPart] = left;
02175             }
02176             overload = left;
02177           }
02178         }
02179         else if (exclusiveLoadUptoMe < overload){
02180           //if the previous processors load is not enough
02181           //then this processor should convert some of its elements.
02182           gno_t mySendCount = overload - exclusiveLoadUptoMe;
02183           //how much more needed.
02184           gno_t sendAfterMe = 0;
02185           //if my load is not enough
02186           //then the next processor will continue to convert.
02187           if (mySendCount > myload){
02188             sendAfterMe = mySendCount - myload;
02189             mySendCount = myload;
02190           }
02191 
02192 
02193           //this processors will convert from overloaded part
02194           //as many as mySendCount items.
02195           mySendFromParts[numPartsISendFrom] = overloadPart;
02196             mySendFromPartsCounts[numPartsISendFrom++] = mySendCount;
02197 
02198             //first consume underloaded parts for the previous processors.
02199           while (exclusiveLoadUptoMe > 0){
02200             int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
02201             gno_t underload = extraInPart[nextUnderLoadedPart];
02202             gno_t left = exclusiveLoadUptoMe + underload;
02203 
02204             if(left >= 0){
02205               extraInPart[nextUnderLoadedPart] = 0;
02206               underloadedPartIndices[underloadpartindex--] = -1;
02207             }
02208             else {
02209               extraInPart[nextUnderLoadedPart] = left;
02210             }
02211             exclusiveLoadUptoMe = left;
02212           }
02213 
02214           //consume underloaded parts for my load.
02215           while (mySendCount > 0){
02216             int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
02217             gno_t underload = extraInPart[nextUnderLoadedPart];
02218             gno_t left = mySendCount + underload;
02219 
02220             if(left >= 0){
02221               mySendParts[numPartsISendTo] = nextUnderLoadedPart;
02222               mySendCountsToParts[numPartsISendTo++] = -underload;
02223 
02224               extraInPart[nextUnderLoadedPart] = 0;
02225               underloadedPartIndices[underloadpartindex--] = -1;
02226             }
02227             else {
02228               extraInPart[nextUnderLoadedPart] = left;
02229 
02230               mySendParts[numPartsISendTo] = nextUnderLoadedPart;
02231               mySendCountsToParts[numPartsISendTo++] = mySendCount;
02232 
02233             }
02234             mySendCount = left;
02235           }
02236           //consume underloaded parts for the load of the processors after my index.
02237           while (sendAfterMe > 0){
02238             int nextUnderLoadedPart = underloadedPartIndices[underloadpartindex];
02239             gno_t underload = extraInPart[nextUnderLoadedPart];
02240             gno_t left = sendAfterMe + underload;
02241 
02242             if(left >= 0){
02243               extraInPart[nextUnderLoadedPart] = 0;
02244               underloadedPartIndices[underloadpartindex--] = -1;
02245             }
02246             else {
02247               extraInPart[nextUnderLoadedPart] = left;
02248             }
02249             sendAfterMe = left;
02250           }
02251         }
02252       }
02253 
02254 
02255     //############################################################//
02257     //############################################################//
02258       for (int i = 0 ; i < numPartsISendFrom; ++i){
02259 
02260         //get the part from which the elements will be converted.
02261         int sendFromPart = mySendFromParts[i];
02262         gno_t sendCount = mySendFromPartsCounts[i];
02263         while(sendCount > 0){
02264           int partToSendIndex = numPartsISendTo - 1;
02265           int partToSend = mySendParts[partToSendIndex];
02266 
02267           int sendCountToThisPart = mySendCountsToParts[partToSendIndex];
02268 
02269           //determine which part i should convert to
02270           //and how many to this part.
02271           if (sendCountToThisPart <= sendCount){
02272             mySendParts[partToSendIndex] = 0;
02273             mySendCountsToParts[partToSendIndex] = 0;
02274             --numPartsISendTo;
02275             sendCount -= sendCountToThisPart;
02276           }
02277           else {
02278             mySendCountsToParts[partToSendIndex] = sendCountToThisPart - sendCount;
02279             sendCountToThisPart = sendCount;
02280             sendCount = 0;
02281           }
02282 
02283 
02284             gno_t toChange = partBegins[sendFromPart];
02285           gno_t previous_begin = partBegins[partToSend];
02286 
02287           //do the conversion.
02288           for (int k = 0; k < sendCountToThisPart - 1; ++k){
02289             coordinate_grid_parts[toChange] = partToSend;
02290             toChange = partNext[toChange];
02291           }
02292           coordinate_grid_parts[toChange] = partToSend;
02293 
02294           gno_t newBegin = partNext[toChange];
02295           partNext[toChange] = previous_begin;
02296           partBegins[partToSend] = partBegins[sendFromPart];
02297           partBegins[sendFromPart] = newBegin;
02298         }
02299       }
02300 
02301       //if (this->myRank == 0) cout << "4" << endl;
02302 
02303 #ifdef printparts
02304 
02305 
02306       for (int i = 0; i < this->worldSize; ++i) numPointsInParts[i] = 0;
02307 
02308       for (int i = 0; i < this->numLocalCoords; ++i){
02309         numPointsInParts[coordinate_grid_parts[i]] += 1;
02310       }
02311 
02312     reduceAll<int, gno_t>( *(this->comm), Teuchos::REDUCE_SUM,
02313                                 this->worldSize,
02314                                 numPointsInParts,
02315                                 numGlobalPointsInParts);
02316     if (this->myRank == 0){
02317       cout << "reassigning" << endl;
02318       gno_t totalSize = 0;
02319       for (int i = 0; i < this->worldSize; ++i){
02320         cout << "me:" << this->myRank << " NumPoints in part:" << i << " is: " << numGlobalPointsInParts[i] << endl;
02321         totalSize += numGlobalPointsInParts[i];
02322       }
02323       cout << "Total:" << totalSize << " ng:" << this->numGlobalCoords << endl;
02324     }
02325 #endif
02326     delete []mySendCountsToParts;
02327     delete []mySendParts;
02328     delete []mySendFromPartsCounts;
02329     delete []mySendFromParts;
02330     delete []underloadedPartIndices;
02331     delete []overloadedPartIndices;
02332     delete []extraInPart;
02333     delete []partNext;
02334     delete []partBegins;
02335     delete []numPointsInPartsInclusiveUptoMyIndex;
02336     delete []numPointsInParts;
02337     delete []numGlobalPointsInParts;
02338 
02339     delete []shiftProcCount;
02340     delete []bestDimProcs;
02341     delete []dim_slices;
02342       delete []minCoords;
02343       delete []maxCoords;
02344   }
02345 
02346   //given the part numbers for each local coordinate,
02347   //distributes the coordinates to the corresponding processors.
02348   void distribute_points(int *coordinate_grid_parts){
02349 
02350     Tpetra::Distributor distributor(comm);
02351     ArrayView<const int> pIds( coordinate_grid_parts, this->numLocalCoords);
02352     /*
02353     for (int i = 0 ; i < this->numLocalCoords; ++i){
02354       cout << "me:" << this->myRank << " to part:" << coordinate_grid_parts[i] << endl;
02355     }
02356     */
02357     gno_t numMyNewGnos = distributor.createFromSends(pIds);
02358 
02359     //cout << "distribution step 1 me:" << this->myRank << " numLocal:"  <<numMyNewGnos << " old:" <<  numLocalCoords << endl;
02360 
02361     this->numLocalCoords = numMyNewGnos;
02362 
02363 
02364     ArrayRCP<T> recvBuf2(distributor.getTotalReceiveLength());
02365 
02366     for (int i = 0; i < this->coordinate_dimension; ++i){
02367       ArrayView<T> s(this->coords[i], this->numLocalCoords);
02368         distributor.doPostsAndWaits<T>(s, 1, recvBuf2());
02369       delete [] this->coords[i];
02370       this->coords[i] = new T[this->numLocalCoords];
02371         for (int j = 0; j < this->numLocalCoords; ++j){
02372           this->coords[i][j] = recvBuf2[j];
02373         }
02374 
02375     }
02376   }
02377 
02378   //calls MJ for p = numProcs
02379   int predistributeMJ(int *coordinate_grid_parts){
02380     int coord_dim = this->coordinate_dimension;
02381 
02382     lno_t numLocalPoints = this->numLocalCoords;
02383     gno_t numGlobalPoints = this->numGlobalCoords;
02384 
02385 
02386     //T **weight = NULL;
02387     //typedef Tpetra::MultiVector<T, lno_t, gno_t, node_t> tMVector_t;
02388     RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
02389         new Tpetra::Map<lno_t, gno_t, node_t> (numGlobalPoints, numLocalPoints, 0, comm));
02390 
02391     Teuchos::Array<Teuchos::ArrayView<const T> > coordView(coord_dim);
02392 
02393 
02394 
02395     for (int i=0; i < coord_dim; i++){
02396       if(numLocalPoints > 0){
02397         Teuchos::ArrayView<const T> a(coords[i], numLocalPoints);
02398         coordView[i] = a;
02399       } else{
02400         Teuchos::ArrayView<const T> a;
02401         coordView[i] = a;
02402       }
02403     }
02404 
02405     RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >tmVector = RCP< Tpetra::MultiVector<T, lno_t, gno_t, node_t> >(
02406         new Tpetra::MultiVector<T, lno_t, gno_t, node_t>( mp, coordView.view(0, coord_dim), coord_dim));
02407 
02408 
02409     RCP<const tMVector_t> coordsConst = Teuchos::rcp_const_cast<const tMVector_t>(tmVector);
02410     vector<const T *> weights;
02411     vector <int> stride;
02412 
02413     typedef Zoltan2::XpetraMultiVectorAdapter<tMVector_t> inputAdapter_t;
02414     //inputAdapter_t ia(coordsConst);
02415     inputAdapter_t ia(coordsConst,weights, stride);
02416 
02417     Teuchos::RCP <Teuchos::ParameterList> params ;
02418     params =RCP <Teuchos::ParameterList> (new Teuchos::ParameterList, true);
02419 
02420 
02421     params->set("algorithm", "multijagged");
02422     params->set("num_global_parts", this->worldSize);
02423 
02424     //TODO we need to fix the setting parts.
02425     //Although MJ sets the parts with
02426     //currently the part setting is not correct when migration is done.
02427     params->set("migration_check_option", 2);
02428 
02429 
02430     Zoltan2::PartitioningProblem<inputAdapter_t> *problem;
02431 
02432 
02433     try {
02434 #ifdef HAVE_ZOLTAN2_MPI
02435       problem = new Zoltan2::PartitioningProblem<inputAdapter_t>(&ia, params.getRawPtr(),
02436           MPI_COMM_WORLD);
02437 #else
02438       problem = new Zoltan2::PartitioningProblem<inputAdapter_t>(&ia, params.getRawPtr());
02439 #endif
02440     }
02441     CATCH_EXCEPTIONS("PartitioningProblem()")
02442 
02443     try {
02444       problem->solve();
02445     }
02446     CATCH_EXCEPTIONS("solve()")
02447 
02448     const typename inputAdapter_t::part_t *partIds = problem->getSolution().getPartList();
02449 
02450     for (int i = 0; i < this->numLocalCoords;++i){
02451       coordinate_grid_parts[i] = partIds[i];
02452       //cout << "me:" << this->myRank << " i:" << i << " goes to:" << partIds[i] << endl;
02453     }
02454     delete problem;
02455     return 0;
02456   }
02457 
02458 #ifdef HAVE_ZOLTAN2_ZOLTAN
02459   //calls RCP for p = numProcs
02460   int predistributeRCB(int *coordinate_grid_parts){
02461     int rank = this->myRank;
02462     int nprocs = this->worldSize;
02463     DOTS<tMVector_t> dots_;
02464 
02465     MEMORY_CHECK(rank==0 || rank==nprocs-1, "After initializing MPI");
02466 
02467 
02468     int nWeights = 0;
02469     int debugLevel=0;
02470     string memoryOn("memoryOn");
02471     string memoryOff("memoryOff");
02472     bool doMemory=false;
02473     int numGlobalParts = nprocs;
02474     int dummyTimer=0;
02475     bool remap=0;
02476 
02477     string balanceCount("balance_object_count");
02478     string balanceWeight("balance_object_weight");
02479     string mcnorm1("multicriteria_minimize_total_weight");
02480     string mcnorm2("multicriteria_balance_total_maximum");
02481     string mcnorm3("multicriteria_minimize_maximum_weight");
02482     string objective(balanceWeight);   // default
02483 
02484     // Process command line input
02485     CommandLineProcessor commandLine(false, true);
02486     //commandLine.setOption("size", &numGlobalCoords,
02487     //  "Approximate number of global coordinates.");
02488     int input_option = 0;
02489     commandLine.setOption("input_option", &input_option,
02490       "whether to use mesh creation, geometric generator, or file input");
02491     string inputFile = "";
02492 
02493     commandLine.setOption("input_file", &inputFile,
02494       "the input file for geometric generator or file input");
02495 
02496 
02497     commandLine.setOption("size", &numGlobalCoords,
02498       "Approximate number of global coordinates.");
02499     commandLine.setOption("numParts", &numGlobalParts,
02500       "Number of parts (default is one per proc).");
02501     commandLine.setOption("nWeights", &nWeights,
02502       "Number of weights per coordinate, zero implies uniform weights.");
02503     commandLine.setOption("debug", &debugLevel, "Zoltan1 debug level");
02504     commandLine.setOption("remap", "no-remap", &remap,
02505       "Zoltan1 REMAP parameter; disabled by default for scalability testing");
02506     commandLine.setOption("timers", &dummyTimer, "ignored");
02507     commandLine.setOption(memoryOn.c_str(), memoryOff.c_str(), &doMemory,
02508       "do memory profiling");
02509 
02510     string doc(balanceCount);
02511     doc.append(": ignore weights\n");
02512     doc.append(balanceWeight);
02513     doc.append(": balance on first weight\n");
02514     doc.append(mcnorm1);
02515     doc.append(": given multiple weights, balance their total.\n");
02516     doc.append(mcnorm3);
02517     doc.append(": given multiple weights, "
02518                "balance the maximum for each coordinate.\n");
02519     doc.append(mcnorm2);
02520     doc.append(": given multiple weights, balance the L2 norm of the weights.\n");
02521     commandLine.setOption("objective", &objective,  doc.c_str());
02522 
02523     CommandLineProcessor::EParseCommandLineReturn rc =
02524       commandLine.parse(0, NULL);
02525 
02526 
02527 
02528     if (rc != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
02529       if (rc == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) {
02530         if (rank==0) cout << "PASS" << endl;
02531         return 1;
02532       }
02533       else {
02534         if (rank==0) cout << "FAIL" << endl;
02535         return 0;
02536       }
02537     }
02538 
02539     //MEMORY_CHECK(doMemory && rank==0, "After processing parameters");
02540 
02541     // Create the data structure
02542 
02543     int coord_dim = this->coordinate_dimension;
02544 
02545 
02546     RCP<Tpetra::Map<lno_t, gno_t, node_t> > mp = rcp(
02547         new Tpetra::Map<lno_t, gno_t, node_t> (this->numGlobalCoords, this->numLocalCoords, 0, this->comm));
02548 
02549     Teuchos::Array<Teuchos::ArrayView<const T> > coordView(coord_dim);
02550     for (int i=0; i < coord_dim; i++){
02551       if(numLocalCoords > 0){
02552         Teuchos::ArrayView<const scalar_t> a(coords[i], numLocalCoords);
02553         coordView[i] = a;
02554       } else{
02555         Teuchos::ArrayView<const scalar_t> a;
02556         coordView[i] = a;
02557       }
02558     }
02559 
02560     tMVector_t *tmVector = new tMVector_t( mp, coordView.view(0, coord_dim), coord_dim);
02561 
02562     dots_.coordinates = tmVector;
02563     dots_.weights.resize(nWeights);
02564 
02565 
02566     MEMORY_CHECK(doMemory && rank==0, "After creating input");
02567 
02568     // Now call Zoltan to partition the problem.
02569 
02570     float ver;
02571     int aok = Zoltan_Initialize(0,NULL, &ver);
02572 
02573     if (aok != 0){
02574       printf("Zoltan_Initialize failed\n");
02575       exit(0);
02576     }
02577 
02578     struct Zoltan_Struct *zz;
02579     zz = Zoltan_Create(MPI_COMM_WORLD);
02580 
02581     Zoltan_Set_Param(zz, "LB_METHOD", "RCB");
02582     Zoltan_Set_Param(zz, "LB_APPROACH", "PARTITION");
02583     Zoltan_Set_Param(zz, "CHECK_GEOM", "0");
02584     Zoltan_Set_Param(zz, "NUM_GID_ENTRIES", "1");
02585     Zoltan_Set_Param(zz, "NUM_LID_ENTRIES", "0");
02586     Zoltan_Set_Param(zz, "RETURN_LISTS", "PART");
02587     std::ostringstream oss;
02588     oss << numGlobalParts;
02589     Zoltan_Set_Param(zz, "NUM_GLOBAL_PARTS", oss.str().c_str());
02590     oss.str("");
02591     oss << debugLevel;
02592     Zoltan_Set_Param(zz, "DEBUG_LEVEL", oss.str().c_str());
02593 
02594     if (remap)
02595       Zoltan_Set_Param(zz, "REMAP", "1");
02596     else
02597       Zoltan_Set_Param(zz, "REMAP", "0");
02598 
02599     if (objective != balanceCount){
02600       oss.str("");
02601       oss << nWeights;
02602       Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", oss.str().c_str());
02603 
02604       if (objective == mcnorm1)
02605         Zoltan_Set_Param(zz, "RCB_MULTICRITERIA_NORM", "1");
02606       else if (objective == mcnorm2)
02607         Zoltan_Set_Param(zz, "RCB_MULTICRITERIA_NORM", "2");
02608       else if (objective == mcnorm3)
02609         Zoltan_Set_Param(zz, "RCB_MULTICRITERIA_NORM", "3");
02610     }
02611     else{
02612       Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", "0");
02613     }
02614 
02615     Zoltan_Set_Num_Obj_Fn(zz, getNumObj<tMVector_t>, &dots_);
02616     Zoltan_Set_Obj_List_Fn(zz, getObjList<tMVector_t>, &dots_);
02617     Zoltan_Set_Num_Geom_Fn(zz,  getDim<tMVector_t>, &dots_);
02618     Zoltan_Set_Geom_Multi_Fn(zz, getCoords<tMVector_t>, &dots_);
02619 
02620     int changes, numGidEntries, numLidEntries, numImport, numExport;
02621     ZOLTAN_ID_PTR importGlobalGids, importLocalGids;
02622     ZOLTAN_ID_PTR exportGlobalGids, exportLocalGids;
02623     int *importProcs, *importToPart, *exportProcs, *exportToPart;
02624 
02625     MEMORY_CHECK(doMemory && rank==0, "Before Zoltan_LB_Partition");
02626 
02627 
02628     aok = Zoltan_LB_Partition(zz, &changes, &numGidEntries, &numLidEntries,
02629                               &numImport, &importGlobalGids, &importLocalGids,
02630                               &importProcs, &importToPart,
02631                               &numExport, &exportGlobalGids, &exportLocalGids,
02632                               &exportProcs, &exportToPart);
02633 
02634 
02635     MEMORY_CHECK(doMemory && rank==0, "After Zoltan_LB_Partition");
02636 
02637     for (lno_t i = 0; i < numLocalCoords; i++)
02638       coordinate_grid_parts[i] = exportToPart[i];
02639     Zoltan_Destroy(&zz);
02640     MEMORY_CHECK(doMemory && rank==0, "After Zoltan_Destroy");
02641 
02642     delete dots_.coordinates;
02643     return 0;
02644 }
02645 #endif
02646   void redistribute(){
02647     int *coordinate_grid_parts = new int[this->numLocalCoords];
02648     switch (this->predistribution){
02649     case 1:
02650 #ifdef HAVE_ZOLTAN2_ZOLTAN
02651       this->predistributeRCB(coordinate_grid_parts);
02652       break;
02653 #endif
02654     case 2:
02655 
02656       this->predistributeMJ(coordinate_grid_parts);
02657       break;
02658     case 3:
02659       //block
02660       blockPartition(coordinate_grid_parts);
02661       break;
02662     }
02663     this->distribute_points(coordinate_grid_parts);
02664 
02665     delete []coordinate_grid_parts;
02666 
02667 
02668   }
02669 
02670   //############################################################//
02672   //############################################################//
02673 
02674 
02675   int getNumWeights(){
02676     return this->numWeightsPerCoord;
02677   }
02678   int getCoordinateDimension(){
02679     return this->coordinate_dimension;
02680   }
02681   lno_t getNumLocalCoords(){
02682     return this->numLocalCoords;
02683   }
02684   gno_t getNumGlobalCoords(){
02685     return this->numGlobalCoords;
02686   }
02687 
02688   T **getLocalCoordinatesView(){
02689     return this->coords;
02690   }
02691 
02692   T **getLocalWeightsView(){
02693     return this->wghts;
02694   }
02695 
02696   void getLocalCoordinatesCopy( T ** c){
02697     for(int ii = 0; ii < this->coordinate_dimension; ++ii){
02698 #ifdef HAVE_ZOLTAN2_OMP
02699 #pragma omp parallel for
02700 #endif
02701       for (lno_t i = 0; i < this->numLocalCoords; ++i){
02702         c[ii][i] = this->coords[ii][i];
02703       }
02704     }
02705   }
02706 
02707   void getLocalWeightsCopy(T **w){
02708     for(int ii = 0; ii < this->numWeightsPerCoord; ++ii){
02709 #ifdef HAVE_ZOLTAN2_OMP
02710 #pragma omp parallel for
02711 #endif
02712       for (lno_t i = 0; i < this->numLocalCoords; ++i){
02713         w[ii][i] = this->wghts[ii][i];
02714       }
02715     }
02716   }
02717 };
02718 }