Zoltan2
Zoltan2_PartitioningProblem.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 
00050 #ifndef _ZOLTAN2_PARTITIONINGPROBLEM_HPP_
00051 #define _ZOLTAN2_PARTITIONINGPROBLEM_HPP_
00052 
00053 #include <Zoltan2_Problem.hpp>
00054 #include <Zoltan2_PartitioningAlgorithms.hpp>
00055 #include <Zoltan2_PartitioningSolution.hpp>
00056 #include <Zoltan2_PartitioningSolutionQuality.hpp>
00057 #include <Zoltan2_GraphModel.hpp>
00058 #include <Zoltan2_IdentifierModel.hpp>
00059 #include <Zoltan2_IntegerRangeList.hpp>
00060 #include <Zoltan2_MachineRepresentation.hpp>
00061 #include <Zoltan2_TaskMapping.hpp>
00062 
00063 #include <unistd.h>
00064 
00065 #ifdef HAVE_ZOLTAN2_OVIS
00066 #include <ovis.h>
00067 #endif
00068 
00069 namespace Zoltan2{
00070 
00094 template<typename Adapter>
00095 class PartitioningProblem : public Problem<Adapter>
00096 {
00097 public:
00098 
00099   typedef typename Adapter::scalar_t scalar_t;
00100   typedef typename Adapter::gid_t gid_t;
00101   typedef typename Adapter::gno_t gno_t;
00102   typedef typename Adapter::lno_t lno_t;
00103   typedef typename Adapter::user_t user_t;
00104   typedef typename Adapter::base_adapter_t base_adapter_t;
00105 
00106 #ifdef HAVE_ZOLTAN2_MPI
00107   typedef Teuchos::OpaqueWrapper<MPI_Comm> mpiWrapper_t;
00108 #endif
00109 
00112   ~PartitioningProblem() {};
00113 
00114 #ifdef HAVE_ZOLTAN2_MPI
00115 
00118   PartitioningProblem(Adapter *A, Teuchos::ParameterList *p, MPI_Comm comm); 
00119 
00120 #endif
00121 
00123   PartitioningProblem(Adapter *A, Teuchos::ParameterList *p) ;
00124 
00126   //
00127   //    \param updateInputData   If true this indicates that either
00128   //          this is the first attempt at solution, or that we
00129   //          are computing a new solution and the input data has
00130   //          changed since the previous solution was computed.  
00131   //          By input data we mean coordinates, topology, or weights.
00132   //          If false, this indicates that we are computing a
00133   //          new solution using the same input data was used for
00134   //          the previous solution, even though the parameters
00135   //          may have been changed.
00136   //
00137   //  For the sake of performance, we ask the caller to set \c updateInputData
00138   //  to false if he/she is computing a new solution using the same input data,
00139   //  but different problem parameters, than that which was used to compute 
00140   //  the most recent solution.
00141 
00142   void solve(bool updateInputData=true );
00143  
00145   //
00146   //   \return  a reference to the solution to the most recent solve().
00147 
00148   const PartitioningSolution<Adapter> &getSolution() {
00149     return *(solution_.getRawPtr());
00150   };
00151 
00160   const scalar_t getWeightImbalance(int idx=0) const {
00161     scalar_t imb = 0;
00162     if (!metrics_.is_null())
00163       metrics_->getWeightImbalance(imb, idx);
00164 
00165     return imb;
00166   }
00167 
00172   ArrayRCP<const MetricValues<scalar_t> > getMetrics() const {
00173    if (metrics_.is_null()){
00174       ArrayRCP<const MetricValues<scalar_t> > emptyMetrics;
00175       return emptyMetrics;
00176     }
00177     else
00178       return metrics_->getMetrics();
00179   }
00180 
00186   void printMetrics(std::ostream &os) const {
00187     if (metrics_.is_null())
00188       os << "No metrics available." << endl;
00189     else
00190       metrics_->printMetrics(os);
00191   };
00192 
00231   void setPartSizes(int len, partId_t *partIds, scalar_t *partSizes, 
00232     bool makeCopy=true) 
00233   { 
00234     setPartSizesForCriteria(0, len, partIds, partSizes, makeCopy);
00235   }
00236 
00271   void setPartSizesForCriteria(int criteria, int len, partId_t *partIds,
00272     scalar_t *partSizes, bool makeCopy=true) ;
00273 /*
00274   void setMachine(MachineRepresentation<typename Adapter::base_adapter_t::scalar_t> *machine);
00275 */
00278   void resetParameters(ParameterList *params)
00279   {
00280     Problem<Adapter>::resetParameters(params);  // creates new environment
00281     if (timer_.getRawPtr() != NULL)
00282       this->env_->setTimer(timer_);
00283   }
00284 
00289   const RCP<const Environment> & getEnvironment() const 
00290   {
00291     return this->envConst_;
00292   }
00293 
00294 private:
00295   void initializeProblem();
00296 
00297   void createPartitioningProblem(bool newData);
00298 
00299   RCP<PartitioningSolution<Adapter> > solution_;
00300 
00301   RCP<MachineRepresentation <typename Adapter::base_adapter_t::scalar_t>  > machine_;
00302 
00303   RCP<Comm<int> > problemComm_;
00304   RCP<const Comm<int> > problemCommConst_;
00305 
00306 #ifdef HAVE_ZOLTAN2_MPI
00307   MPI_Comm mpiComm_;
00308 #endif
00309 
00310   BaseAdapterType inputType_;
00311   ModelType modelType_;
00312   modelFlag_t graphFlags_;
00313   modelFlag_t idFlags_;
00314   modelFlag_t coordFlags_;
00315   std::string algorithm_;
00316 
00317   int numberOfWeights_;
00318 
00319   // Suppose Array<partId_t> partIds = partIds_[w].  If partIds.size() > 0
00320   // then the user supplied part sizes for weight index "w", and the sizes
00321   // corresponding to the Ids in partIds are partSizes[w].
00322   //
00323   // If numberOfWeights_ >= 0, then there is an Id and Sizes array for
00324   // for each weight.  Otherwise the user did not supply object weights,
00325   // but they can still specify part sizes. 
00326   // So numberOfCriteria_ is numberOfWeights_ or one, whichever is greater.
00327 
00328   ArrayRCP<ArrayRCP<partId_t> > partIds_;
00329   ArrayRCP<ArrayRCP<scalar_t> > partSizes_;
00330   int numberOfCriteria_;
00331 
00332   // Number of parts to be computed at each level in hierarchical partitioning.
00333   
00334   ArrayRCP<int> levelNumberParts_;
00335   bool hierarchical_;
00336 
00337   // Create a Timer if the user asked for timing stats.
00338 
00339   RCP<TimerManager> timer_;
00340 
00341   // Did the user request metrics?
00342 
00343   bool metricsRequested_;
00344   RCP<const PartitioningSolutionQuality<Adapter> > metrics_;
00345 };
00347 
00348 #ifdef HAVE_ZOLTAN2_MPI
00349 template <typename Adapter>
00350   PartitioningProblem<Adapter>::PartitioningProblem(Adapter *A, 
00351     ParameterList *p, MPI_Comm comm):
00352       Problem<Adapter>(A,p,comm), solution_(),
00353       problemComm_(), problemCommConst_(),
00354       inputType_(InvalidAdapterType), modelType_(InvalidModel), 
00355       graphFlags_(), idFlags_(), coordFlags_(), algorithm_(),
00356       numberOfWeights_(), partIds_(), partSizes_(), 
00357       numberOfCriteria_(), levelNumberParts_(), hierarchical_(false), 
00358       timer_(), metricsRequested_(false), metrics_()
00359 {
00360 
00361   initializeProblem();
00362 }
00363 #endif
00364 /*
00365 template <typename Adapter>
00366 void PartitioningProblem<Adapter>::setMachine(MachineRepresentation<typename Adapter::base_adapter_t::scalar_t> *machine){
00367   this->machine_ = RCP<MachineRepresentation<typename Adapter::base_adapter_t::scalar_t> > (machine, false);
00368 }
00369 */
00370 template <typename Adapter>
00371   PartitioningProblem<Adapter>::PartitioningProblem(Adapter *A, 
00372     ParameterList *p):
00373       Problem<Adapter>(A,p), solution_(),
00374       problemComm_(), problemCommConst_(),
00375       inputType_(InvalidAdapterType), modelType_(InvalidModel), 
00376       graphFlags_(), idFlags_(), coordFlags_(), algorithm_(),
00377       numberOfWeights_(), 
00378       partIds_(), partSizes_(), numberOfCriteria_(), 
00379       levelNumberParts_(), hierarchical_(false), timer_(),
00380       metricsRequested_(false), metrics_()
00381 {
00382   initializeProblem();
00383 }
00384 
00385 template <typename Adapter>
00386   void PartitioningProblem<Adapter>::initializeProblem()
00387 {
00388   HELLO;
00389 
00390   this->env_->debug(DETAILED_STATUS, "PartitioningProblem::initializeProblem");
00391 
00392   if (getenv("DEBUGME")){
00393     std::cout << getpid() << std::endl;
00394     sleep(15);
00395   }
00396 
00397 #ifdef HAVE_ZOLTAN2_OVIS
00398   ovis_enabled(this->comm_->getRank());
00399 #endif
00400 
00401   // Create a copy of the user's communicator.  
00402 
00403   problemComm_ = this->comm_->duplicate();
00404   problemCommConst_ = rcp_const_cast<const Comm<int> > (problemComm_);
00405 
00406   machine_ = RCP <Zoltan2::MachineRepresentation<typename Adapter::scalar_t> >(new Zoltan2::MachineRepresentation<typename Adapter::scalar_t>(problemComm_));
00407 #ifdef HAVE_ZOLTAN2_MPI
00408 
00409   // TPLs may want an MPI communicator
00410 
00411   mpiComm_ = Teuchos2MPI(problemComm_);
00412 
00413 #endif
00414 
00415   // Number of criteria is number of user supplied weights if non-zero.
00416   // Otherwise it is 1 and uniform weight is implied.
00417 
00418   numberOfWeights_ = this->inputAdapter_->getNumWeightsPerID();
00419 
00420   numberOfCriteria_ = (numberOfWeights_ > 1) ? numberOfWeights_ : 1;
00421 
00422   inputType_ = this->inputAdapter_->adapterType();
00423 
00424   // The Caller can specify part sizes in setPartSizes().  If he/she
00425   // does not, the part size arrays are empty.
00426 
00427   ArrayRCP<partId_t> *noIds = new ArrayRCP<partId_t> [numberOfCriteria_];
00428   ArrayRCP<scalar_t> *noSizes = new ArrayRCP<scalar_t> [numberOfCriteria_];
00429 
00430   partIds_ = arcp(noIds, 0, numberOfCriteria_, true);
00431   partSizes_ = arcp(noSizes, 0, numberOfCriteria_, true);
00432 
00433   if (this->env_->getDebugLevel() >= DETAILED_STATUS){
00434     std::ostringstream msg;
00435     msg << problemComm_->getSize() << " procs,"
00436       << numberOfWeights_ << " user-defined weights\n";
00437     this->env_->debug(DETAILED_STATUS, msg.str());
00438   }
00439 
00440   this->env_->memory("After initializeProblem");
00441 }
00442 
00443 template <typename Adapter>
00444   void PartitioningProblem<Adapter>::setPartSizesForCriteria(
00445     int criteria, int len, partId_t *partIds, scalar_t *partSizes, bool makeCopy) 
00446 {
00447   this->env_->localInputAssertion(__FILE__, __LINE__, "invalid length", 
00448     len>= 0, BASIC_ASSERTION);
00449 
00450   this->env_->localInputAssertion(__FILE__, __LINE__, "invalid criteria", 
00451     criteria >= 0 && criteria < numberOfCriteria_, BASIC_ASSERTION);
00452 
00453   if (len == 0){
00454     partIds_[criteria] = ArrayRCP<partId_t>();
00455     partSizes_[criteria] = ArrayRCP<scalar_t>();
00456     return;
00457   }
00458 
00459   this->env_->localInputAssertion(__FILE__, __LINE__, "invalid arrays", 
00460     partIds && partSizes, BASIC_ASSERTION);
00461 
00462   // The global validity of the partIds and partSizes arrays is performed
00463   // by the PartitioningSolution, which computes global part distribution and
00464   // part sizes.
00465 
00466   partId_t *z2_partIds = NULL;
00467   scalar_t *z2_partSizes = NULL;
00468   bool own_memory = false;
00469 
00470   if (makeCopy){
00471     z2_partIds = new partId_t [len];
00472     z2_partSizes = new scalar_t [len];
00473     this->env_->localMemoryAssertion(__FILE__, __LINE__, len, z2_partSizes);
00474     memcpy(z2_partIds, partIds, len * sizeof(partId_t));
00475     memcpy(z2_partSizes, partSizes, len * sizeof(scalar_t));
00476     own_memory=true;
00477   }
00478   else{
00479     z2_partIds = partIds;
00480     z2_partSizes = partSizes;
00481   }
00482 
00483   partIds_[criteria] = arcp(z2_partIds, 0, len, own_memory);
00484   partSizes_[criteria] = arcp(z2_partSizes, 0, len, own_memory);
00485 }
00486 
00487 template <typename Adapter>
00488 void PartitioningProblem<Adapter>::solve(bool updateInputData)
00489 {
00490   HELLO;
00491   this->env_->debug(DETAILED_STATUS, "Entering solve");
00492 
00493   // Create the computational model.
00494 
00495   this->env_->timerStart(MACRO_TIMERS, "create problem");
00496 
00497   createPartitioningProblem(updateInputData);
00498 
00499   this->env_->timerStop(MACRO_TIMERS, "create problem");
00500 
00501   // TODO: If hierarchical_
00502 
00503   // Create the solution. The algorithm will query the Solution
00504   //   for part and weight information. The algorithm will
00505   //   update the solution with part assignments and quality
00506   //   metrics.  The Solution object itself will convert our internal
00507   //   global numbers back to application global Ids if needed.
00508 
00509   RCP<const IdentifierMap<user_t> > idMap = 
00510     this->baseModel_->getIdentifierMap();
00511 
00512   PartitioningSolution<Adapter> *soln = NULL;
00513 
00514   this->env_->timerStart(MACRO_TIMERS, "create solution");
00515 
00516   try{
00517     soln = new PartitioningSolution<Adapter>( 
00518       this->envConst_, problemCommConst_, idMap, numberOfWeights_, 
00519       partIds_.view(0, numberOfCriteria_), 
00520       partSizes_.view(0, numberOfCriteria_));
00521   }
00522   Z2_FORWARD_EXCEPTIONS;
00523 
00524   solution_ = rcp(soln);
00525 
00526   this->env_->timerStop(MACRO_TIMERS, "create solution");
00527 
00528   this->env_->memory("After creating Solution");
00529 
00530   // Call the algorithm
00531 
00532   try {
00533     if (algorithm_ == std::string("scotch")){
00534       AlgPTScotch<Adapter>(this->envConst_, problemComm_,
00535 #ifdef HAVE_ZOLTAN2_MPI
00536         mpiComm_,
00537 #endif
00538         this->graphModel_, solution_);
00539     }
00540     else if (algorithm_ == std::string("block")){
00541       AlgBlock<Adapter> algblock(this->envConst_, problemComm_,
00542         this->identifierModel_, solution_);
00543       algblock.solve();
00544     }
00545     else if (algorithm_ == std::string("rcb")){
00546       AlgRCB<Adapter>(this->envConst_, problemComm_,
00547         this->coordinateModel_, solution_);
00548     }
00549     else if (algorithm_ == std::string("multijagged")){
00550       Zoltan2_AlgMJ<Adapter> *alg_mj = new Zoltan2_AlgMJ<Adapter>();
00551       alg_mj->multi_jagged_part( this->envConst_, problemComm_,
00552          this->coordinateModel_, solution_);
00553       delete alg_mj;
00554     }
00555     else{
00556       throw std::logic_error("partitioning algorithm not supported yet");
00557     }
00558   }
00559   Z2_FORWARD_EXCEPTIONS;
00560 
00561   //if mapping is requested
00562   const Teuchos::ParameterEntry *pe = this->envConst_->getParameters().getEntryPtr("mapping_type");
00563   int mapping_type = -1;
00564   if (pe){
00565     mapping_type = pe->getValue(&mapping_type);
00566   }
00567   //if mapping is 0 -- coordinate mapping
00568 
00569   if (mapping_type == 0){
00570 
00571     //partId_t *task_communication_xadj = NULL, *task_communication_adj = NULL;
00572     Zoltan2::CoordinateTaskMapper <Adapter, zoltan2_partId_t> *ctm=
00573                   new Zoltan2::CoordinateTaskMapper<Adapter,zoltan2_partId_t>(
00574                           problemComm_.getRawPtr(),
00575                           machine_.getRawPtr(),
00576                           this->coordinateModel_.getRawPtr(),
00577                           solution_.getRawPtr(),
00578                           this->envConst_.getRawPtr()
00579                           //,task_communication_xadj,
00580                           //task_communication_adj
00581                           );
00582     //for now just delete the object.
00583     delete ctm;
00584   }
00585   else if (mapping_type == 1){
00586     //if mapping is 1 -- graph mapping
00587   }
00588 
00589 #ifdef KDDKDD_SHOULD_NEVER_CHANGE_PROBLEMCOMM
00590 #ifdef HAVE_ZOLTAN2_MPI
00591 
00592   // The algorithm may have changed the communicator.  Change it back.
00593   // KDD:  Why would the algorithm change the communicator? TODO
00594   // KDD:  Should we allow such a side effect? TODO
00595 
00596   RCP<const mpiWrapper_t > wrappedComm = rcp(new mpiWrapper_t(mpiComm_));
00597   problemComm_ = rcp(new Teuchos::MpiComm<int>(wrappedComm));
00598   problemCommConst_ = rcp_const_cast<const Comm<int> > (problemComm_);
00599 
00600 #endif
00601 #endif
00602 
00603   if (metricsRequested_){
00604     typedef PartitioningSolution<Adapter> ps_t;
00605     typedef PartitioningSolutionQuality<Adapter> psq_t;
00606 
00607     psq_t *quality = NULL;
00608     RCP<const ps_t> solutionConst = rcp_const_cast<const ps_t>(solution_);
00609     RCP<const Adapter> adapter = rcp(this->inputAdapter_, false);
00610 
00611     try{
00612       quality = new psq_t(this->envConst_, problemCommConst_, adapter, 
00613         solutionConst);
00614     }
00615     Z2_FORWARD_EXCEPTIONS
00616 
00617     metrics_ = rcp(quality);
00618   }
00619 
00620   this->env_->debug(DETAILED_STATUS, "Exiting solve");
00621 }
00622 
00623 template <typename Adapter>
00624 void PartitioningProblem<Adapter>::createPartitioningProblem(bool newData)
00625 {
00626   HELLO;
00627   this->env_->debug(DETAILED_STATUS, 
00628     "PartitioningProblem::createPartitioningProblem");
00629 
00630   using std::string;
00631   using Teuchos::ParameterList;
00632 
00633   // A Problem object may be reused.  The input data may have changed and
00634   // new parameters or part sizes may have been set.
00635   //
00636   // Save these values in order to determine if we need to create a new model.
00637 
00638   ModelType previousModel = modelType_;
00639   modelFlag_t previousGraphModelFlags = graphFlags_;
00640   modelFlag_t previousIdentifierModelFlags = idFlags_;
00641   modelFlag_t previousCoordinateModelFlags = coordFlags_;
00642 
00643   modelType_ = InvalidModel;
00644   graphFlags_.reset();
00645   idFlags_.reset();
00646   coordFlags_.reset();
00647 
00649   // It's possible at this point that the Problem may want to
00650   // add problem parameters to the parameter list in the Environment. 
00651   //
00652   // Since the parameters in the Environment have already been
00653   // validated in its constructor, a new Environment must be created:
00655   // Teuchos::RCP<const Teuchos::Comm<int> > oldComm = this->env_->comm_;
00656   // const ParameterList &oldParams = this->env_->getUnvalidatedParameters();
00657   // 
00658   // ParameterList newParams = oldParams;
00659   // newParams.set("new_parameter", "new_value");
00660   // 
00661   // ParameterList &newPartParams = newParams.sublist("partitioning");
00662   // newPartParams.set("new_partitioning_parameter", "its_value");
00663   // 
00664   // this->env_ = rcp(new Environment(newParams, oldComm));
00666 
00667   this->env_->debug(DETAILED_STATUS, "    parameters");
00668   Environment &env = *(this->env_);
00669   ParameterList &pl = env.getParametersNonConst();
00670 
00671   std::string defString("default");
00672 
00673   // Did the user ask for computation of quality metrics?
00674 
00675   int yesNo=0;
00676   const Teuchos::ParameterEntry *pe = pl.getEntryPtr("compute_metrics");
00677   if (pe){
00678     yesNo = pe->getValue<int>(&yesNo);
00679     metricsRequested_ = true;
00680   }
00681 
00682   // Did the user specify a computational model?
00683 
00684   std::string model(defString);
00685   pe = pl.getEntryPtr("model");
00686   if (pe)
00687     model = pe->getValue<std::string>(&model);
00688 
00689   // Did the user specify an algorithm?
00690 
00691   std::string algorithm(defString);
00692   pe = pl.getEntryPtr("algorithm");
00693   if (pe)
00694     algorithm = pe->getValue<std::string>(&algorithm);
00695 
00696   // Possible algorithm requirements that must be conveyed to the model:
00697 
00698   bool needConsecutiveGlobalIds = false;
00699   bool removeSelfEdges= false;
00700 
00702   // Determine algorithm, model, and algorithm requirements.  This
00703   // is a first pass.  Feel free to change this and add to it.
00704   
00705   if (algorithm != defString){
00706     // Figure out the model required by the algorithm
00707     if (algorithm == std::string("block") ||
00708         algorithm == std::string("random") ||
00709         algorithm == std::string("cyclic") ){
00710 
00711       modelType_ = IdentifierModelType;
00712       algorithm_ = algorithm;
00713       needConsecutiveGlobalIds = true;
00714     }
00715     else if (algorithm == std::string("rcb") ||
00716              algorithm == std::string("rib") ||
00717              algorithm == std::string("multijagged") ||
00718              algorithm == std::string("hsfc")){
00719 
00720       modelType_ = CoordinateModelType;
00721       algorithm_ = algorithm;
00722     }
00723     else if (algorithm == std::string("metis") ||
00724              algorithm == std::string("parmetis") ||
00725              algorithm == std::string("scotch") ||
00726              algorithm == std::string("ptscotch")){
00727 
00728       modelType_ = GraphModelType;
00729       algorithm_ = algorithm;
00730       removeSelfEdges = true;
00731       needConsecutiveGlobalIds = true;
00732     }
00733     else if (algorithm == std::string("patoh") ||
00734              algorithm == std::string("phg")){
00735 
00736       if ((modelType_ != GraphModelType) &&
00737           (modelType_ != HypergraphModelType) ){
00738         modelType_ = HypergraphModelType;
00739       }
00740       algorithm_ = algorithm;
00741       needConsecutiveGlobalIds = true;
00742     }
00743     else{
00744       // Parameter list should ensure this does not happen.
00745       throw std::logic_error("parameter list algorithm is invalid");
00746     }
00747   }
00748   else if (model != defString){
00749     // Figure out the algorithm suggested by the model.
00750     if (model == std::string("hypergraph")){
00751       modelType_ = HypergraphModelType;
00752       if (problemComm_->getSize() > 1)
00753         algorithm_ = std::string("phg"); 
00754       else
00755         algorithm_ = std::string("patoh"); 
00756       needConsecutiveGlobalIds = true;
00757     }
00758     else if (model == std::string("graph")){
00759       modelType_ = GraphModelType;
00760 #ifdef HAVE_ZOLTAN2_SCOTCH
00761       if (problemComm_->getSize() > 1)
00762         algorithm_ = std::string("ptscotch"); 
00763       else
00764         algorithm_ = std::string("scotch"); 
00765       removeSelfEdges = true;
00766       needConsecutiveGlobalIds = true;
00767 #else
00768 #ifdef HAVE_ZOLTAN2_PARMETIS
00769       if (problemComm_->getSize() > 1)
00770         algorithm_ = std::string("parmetis"); 
00771       else
00772         algorithm_ = std::string("metis"); 
00773       removeSelfEdges = true;
00774       needConsecutiveGlobalIds = true;
00775 #else
00776       if (problemComm_->getSize() > 1)
00777         algorithm_ = std::string("phg"); 
00778       else
00779         algorithm_ = std::string("patoh"); 
00780       removeSelfEdges = true;
00781       needConsecutiveGlobalIds = true;
00782 #endif
00783 #endif
00784     }
00785     else if (model == std::string("geometry")){
00786       modelType_ = CoordinateModelType;
00787       algorithm_ = std::string("rcb");
00788     }
00789     else if (model == std::string("ids")){
00790       modelType_ = IdentifierModelType;
00791       algorithm_ = std::string("block");
00792       needConsecutiveGlobalIds = true;
00793     }
00794     else{
00795       // Parameter list should ensure this does not happen.
00796       env.localBugAssertion(__FILE__, __LINE__, 
00797         "parameter list model type is invalid", 1, BASIC_ASSERTION);
00798     }
00799   }
00800   else{   
00801     // Determine an algorithm and model suggested by the input type.
00802     //   TODO: this is a good time to use the time vs. quality parameter
00803     //     in choosing an algorithm, and setting some parameters
00804 
00805     if (inputType_ == MatrixAdapterType){
00806       modelType_ = HypergraphModelType;
00807       if (problemComm_->getSize() > 1)
00808         algorithm_ = std::string("phg"); 
00809       else
00810         algorithm_ = std::string("patoh"); 
00811     }
00812     else if (inputType_ == GraphAdapterType ||
00813         inputType_ == MeshAdapterType){
00814       modelType_ = GraphModelType;
00815       if (problemComm_->getSize() > 1)
00816         algorithm_ = std::string("phg"); 
00817       else
00818         algorithm_ = std::string("patoh"); 
00819     }
00820     else if (inputType_ == CoordinateAdapterType){
00821       modelType_ = CoordinateModelType;
00822       if(algorithm_ != std::string("multijagged"))
00823       algorithm_ = std::string("rcb");
00824     }
00825     else if (inputType_ == VectorAdapterType ||
00826              inputType_ == IdentifierAdapterType){
00827       modelType_ = IdentifierModelType;
00828       algorithm_ = std::string("block");
00829     }
00830     else{
00831       // This should never happen
00832       throw std::logic_error("input type is invalid");
00833     }
00834   }
00835 
00836   // Hierarchical partitioning?
00837 
00838   Array<int> valueList;
00839   pe = pl.getEntryPtr("topology");
00840 
00841   if (pe){
00842     valueList = pe->getValue<Array<int> >(&valueList);
00843 
00844     if (!Zoltan2::noValuesAreInRangeList<int>(valueList)){
00845       int *n = new int [valueList.size() + 1];
00846       levelNumberParts_ = arcp(n, 0, valueList.size() + 1, true);
00847       int procsPerNode = 1;
00848       for (int i=0; i < valueList.size(); i++){
00849         levelNumberParts_[i+1] = valueList[i];
00850         procsPerNode *= valueList[i];
00851       }
00852       // Number of parts in the first level
00853       levelNumberParts_[0] = env.numProcs_ / procsPerNode;
00854 
00855       if (env.numProcs_ % procsPerNode > 0)
00856         levelNumberParts_[0]++;
00857     }
00858   }
00859   else{
00860     levelNumberParts_.clear();
00861   }
00862 
00863   hierarchical_ = levelNumberParts_.size() > 0;
00864 
00865   // Object to be partitioned? (rows, columns, etc)
00866 
00867   std::string objectOfInterest(defString);
00868   pe = pl.getEntryPtr("objects_to_partition");
00869   if (pe)
00870     objectOfInterest = pe->getValue<std::string>(&objectOfInterest);
00871 
00873   // Set model creation flags, if any.
00874 
00875   this->env_->debug(DETAILED_STATUS, "    models");
00876   if (modelType_ == GraphModelType){
00877 
00878     // Any parameters in the graph sublist?
00879 
00880     std::string symParameter(defString);
00881     pe = pl.getEntryPtr("symmetrize_graph");
00882     if (pe){
00883       symParameter = pe->getValue<std::string>(&symParameter);
00884       if (symParameter == std::string("transpose"))
00885         graphFlags_.set(SYMMETRIZE_INPUT_TRANSPOSE);
00886       else if (symParameter == std::string("bipartite"))
00887         graphFlags_.set(SYMMETRIZE_INPUT_BIPARTITE);
00888     } 
00889 
00890     int sgParameter = 0;
00891     pe = pl.getEntryPtr("subset_graph");
00892     if (pe)
00893       sgParameter = pe->getValue<int>(&sgParameter);
00894 
00895     if (sgParameter == 1)
00896         graphFlags_.set(GRAPH_IS_A_SUBSET_GRAPH);
00897 
00898     // Any special behaviors required by the algorithm?
00899     
00900     if (removeSelfEdges)
00901       graphFlags_.set(SELF_EDGES_MUST_BE_REMOVED);
00902 
00903     if (needConsecutiveGlobalIds)
00904       graphFlags_.set(IDS_MUST_BE_GLOBALLY_CONSECUTIVE);
00905 
00906     // How does user input map to vertices and edges?
00907 
00908     if (inputType_ == MatrixAdapterType){
00909       if (objectOfInterest == defString ||
00910           objectOfInterest == std::string("matrix_rows") )
00911         graphFlags_.set(VERTICES_ARE_MATRIX_ROWS);
00912       else if (objectOfInterest == std::string("matrix_columns"))
00913         graphFlags_.set(VERTICES_ARE_MATRIX_COLUMNS);
00914       else if (objectOfInterest == std::string("matrix_nonzeros"))
00915         graphFlags_.set(VERTICES_ARE_MATRIX_NONZEROS);
00916     }
00917 
00918     else if (inputType_ == MeshAdapterType){
00919       if (objectOfInterest == defString ||
00920           objectOfInterest == std::string("mesh_nodes") )
00921         graphFlags_.set(VERTICES_ARE_MESH_NODES);
00922       else if (objectOfInterest == std::string("mesh_elements"))
00923         graphFlags_.set(VERTICES_ARE_MESH_ELEMENTS);
00924     } 
00925   }
00926   else if (modelType_ == IdentifierModelType){
00927 
00928     // Any special behaviors required by the algorithm?
00929     
00930     if (needConsecutiveGlobalIds)
00931       idFlags_.set(IDS_MUST_BE_GLOBALLY_CONSECUTIVE);
00932   }
00933   else if (modelType_ == CoordinateModelType){
00934 
00935     // Any special behaviors required by the algorithm?
00936     
00937     if (needConsecutiveGlobalIds)
00938       coordFlags_.set(IDS_MUST_BE_GLOBALLY_CONSECUTIVE);
00939   }
00940 
00941 
00942   if (  newData ||
00943        (modelType_ != previousModel) ||
00944        (graphFlags_ != previousGraphModelFlags) ||
00945        (coordFlags_ != previousCoordinateModelFlags) ||
00946        (idFlags_ != previousIdentifierModelFlags) ) {
00947 
00948     // Create the computational model.
00949     // Models are instantiated for base input adapter types (mesh,
00950     // matrix, graph, and so on).  We pass a pointer to the input
00951     // adapter, cast as the base input type.
00952 
00953     //KDD Not sure why this shadow declaration is needed
00954     //KDD Comment out for now; revisit later if problems.
00955     //KDD const Teuchos::ParameterList pl = this->envConst_->getParameters();
00956     //bool exceptionThrow = true;
00957 
00958     switch (modelType_) {
00959 
00960     case GraphModelType:
00961       this->env_->debug(DETAILED_STATUS, "    building graph model");
00962       this->graphModel_ = rcp(new GraphModel<base_adapter_t>(
00963         this->baseInputAdapter_, this->envConst_, problemComm_, graphFlags_));
00964 
00965       this->baseModel_ = rcp_implicit_cast<const Model<base_adapter_t> >(
00966         this->graphModel_);
00967 
00968       break;
00969 
00970     case HypergraphModelType:
00971       break;
00972   
00973     case CoordinateModelType:
00974       this->env_->debug(DETAILED_STATUS, "    building coordinate model");
00975       this->coordinateModel_ = rcp(new CoordinateModel<base_adapter_t>(
00976         this->baseInputAdapter_, this->envConst_, problemComm_, coordFlags_));
00977 
00979       // It's possible at this point that the Problem may want to
00980       // add problem parameters to the parameter list in the Environment.
00981       //
00982       // Since the parameters in the Environment have already been
00983       // validated in its constructor, a new Environment must be created:
00985       // Teuchos::RCP<const Teuchos::Comm<int> > oldComm = this->env_->comm_;
00986       // const ParameterList &oldParams = this->env_->getUnvalidatedParameters();
00987       //
00988       // ParameterList newParams = oldParams;
00989       // newParams.set("new_parameter", "new_value");
00990       //
00991       // ParameterList &newPartParams = newParams.sublist("partitioning");
00992       // newPartParams.set("new_partitioning_parameter", "its_value");
00993       //
00994       // this->env_ = rcp(new Environment(newParams, oldComm));
00996 /*
00997       if(algorithm == string("multijagged")){
00998 
00999           //int coordinateCnt = this->coordinateModel_->getCoordinateDim();
01000           //cout << coordinateCnt << " " << pl.getPtr<Array <int> >("pqParts")->size() << endl;
01001           //exceptionThrow = coordinateCnt == pl.getPtr<Array <int> >("pqParts")->size();
01002           int arraySize = pl.getPtr<Array <int> >("pqParts")->size() - 1;
01003           exceptionThrow = arraySize > 0;
01004           this->envConst_->localInputAssertion(__FILE__, __LINE__, "invalid length of cut lines. Size of cut lines should be at least 1.",
01005                         exceptionThrow, BASIC_ASSERTION);
01006 
01007 
01008           int totalPartCount = 1;
01009           for(int i = 0; i < arraySize; ++i){
01010             //cout <<  pl.getPtr<Array <int> >("pqParts")->getRawPtr()[i] << " ";
01011             totalPartCount *= pl.getPtr<Array <int> >("pqParts")->getRawPtr()[i];
01012 // TODO:  Using pointer in parameter list.   Ross says, "Bad."  Can't print it.
01013           }
01014           Teuchos::ParameterList newParams = pl;
01015 // TODO:  KDD I thought we got rid of sublists in the parameter list??
01016           Teuchos::ParameterList &parParams = newParams.sublist("partitioning");
01017 
01018 // TODO:  KDD Is there a more elegant solution here than changing the paramlist?
01019 // TODO:  How does this even work?  Where is newParams used?
01020 
01021           parParams.set("num_global_parts", totalPartCount);
01022 
01023 
01024           //cout << endl;
01025           Teuchos::RCP<const Teuchos::Comm<int> > oldComm = this->envConst_->comm_;
01026 
01027           //this->envConst_ = rcp(new Environment(newParams, oldComm));
01028 
01029 
01030       }
01031 */
01032 
01033       this->baseModel_ = rcp_implicit_cast<const Model<base_adapter_t> >(
01034         this->coordinateModel_);
01035       break;
01036 
01037     case IdentifierModelType:
01038       this->env_->debug(DETAILED_STATUS, "    building identifier model");
01039       this->identifierModel_ = rcp(new IdentifierModel<base_adapter_t>(
01040         this->baseInputAdapter_, this->envConst_, problemComm_, idFlags_));
01041 
01042       this->baseModel_ = rcp_implicit_cast<const Model<base_adapter_t> >(
01043         this->identifierModel_);
01044       break;
01045 
01046     default:
01047       cout << __func__ << " Invalid model" << modelType_ << endl;
01048       break;
01049     }
01050 
01051     this->env_->memory("After creating Model");
01052     this->env_->debug(DETAILED_STATUS, "createPartitioningProblem done");
01053   }
01054 
01055   /*
01056   Teuchos::RCP<const Teuchos::Comm<int> > oldComm = this->env_->comm_;
01057   const ParameterList &oldParams = this->env_->getUnvalidatedParameters();
01058 
01059   ParameterList newParams = oldParams;
01060   int totalPartCount = 1;
01061   const int *partNo = pl.getPtr<Array <int> >("pqParts")->getRawPtr();
01062 
01063   for (int i = 0; i < coordDim; ++i){
01064     totalPartCount *= partNo[i];
01065   }
01066   newParams.set("num_global_parts", totalPartCount);
01067 
01068 
01069   this->env_ = rcp(new Environment(newParams, oldComm));
01070 */
01071 }
01072 
01073 }  // namespace Zoltan2
01074 #endif