Zoltan2
Zoltan2_AlgRCB.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_ALGRCB_HPP_
00051 #define _ZOLTAN2_ALGRCB_HPP_
00052 
00053 #include <Zoltan2_AlgRCB_methods.hpp>
00054 #include <Zoltan2_CoordinateModel.hpp>
00055 #include <Zoltan2_Exceptions.hpp>
00056 
00057 #include <Teuchos_ParameterList.hpp>
00058 
00059 namespace Zoltan2{
00060 
00082 template <typename Adapter>
00083 void AlgRCB(
00084   const RCP<const Environment> &env,
00085   const RCP<Comm<int> > &problemComm,
00086   const RCP<const CoordinateModel<
00087     typename Adapter::base_adapter_t> > &coords, 
00088   RCP<PartitioningSolution<Adapter> > &solution
00089 ) 
00090 {
00091 #ifndef INCLUDE_ZOLTAN2_EXPERIMENTAL
00092 
00093   Z2_THROW_EXPERIMENTAL("Zoltan2 RCB is strictly experimental software "
00094                         "due to performance problems in its use of Tpetra.")
00095 
00096 #else  // INCLUDE_ZOLTAN2_EXPERIMENTAL
00097 
00098   typedef typename Adapter::node_t node_t;
00099   typedef typename Adapter::lno_t lno_t;
00100   typedef typename Adapter::gno_t gno_t;
00101   typedef typename Adapter::part_t part_t;
00102   typedef typename Adapter::scalar_t scalar_t;
00103 
00104   // Make a copy of communicator because
00105   // we subdivide the communicator during the algorithm.
00106 
00107   RCP<Comm<int> > comm = problemComm->duplicate();
00108 
00109   std::bitset<NUM_RCB_PARAMS> params;
00110 
00111   env->debug(DETAILED_STATUS, "Entering AlgPartRCB");
00112 
00113   const Teuchos::ParameterList &pl = env->getParameters();
00114 
00115   env->timerStart(BOTH_TIMERS, "RCB set up");
00116 
00118   // Partitioning problem parameters of interest:
00119   //    objective
00120   //    imbalance_tolerance
00121 
00122   env->debug(DETAILED_STATUS, "Accessing parameters");
00123   multiCriteriaNorm mcnorm = normBalanceTotalMaximum;
00124   std::string obj;
00125 
00126   const Teuchos::ParameterEntry *pe = pl.getEntryPtr("partitioning_objective");
00127   if (pe)
00128     obj = pe->getValue(&obj);
00129 
00130   if (!pe){
00131     params.set(rcb_balanceWeight);
00132     mcnorm = normBalanceTotalMaximum;
00133   }
00134   else if (obj == std::string("balance_object_count")){
00135     params.set(rcb_balanceCount);
00136   }
00137   else if (obj == std::string("multicriteria_minimize_total_weight")){
00138     params.set(rcb_minTotalWeight);
00139     mcnorm = normMinimizeTotalWeight;
00140   }
00141   else if (obj == std::string("multicriteria_minimize_maximum_weight")){
00142     params.set(rcb_minMaximumWeight);
00143     mcnorm = normMinimizeMaximumWeight;
00144   }
00145   else if (obj == std::string("multicriteria_balance_total_maximum")){
00146     params.set(rcb_balanceTotalMaximum);
00147     mcnorm = normBalanceTotalMaximum;
00148   }
00149   else{
00150     params.set(rcb_balanceWeight);
00151     mcnorm = normBalanceTotalMaximum;
00152   }
00153 
00154   scalar_t imbalanceTolerance = .1;
00155   pe = pl.getEntryPtr("imbalance_tolerance");
00156   if (pe){
00157     double tol;
00158     tol = pe->getValue(&tol);
00159     imbalanceTolerance = tol - 1.0;
00160   }
00161 
00162   if (imbalanceTolerance <= 0)
00163     imbalanceTolerance = 10e-4;  // TODO - what's a good choice
00164 
00166   // Geometric partitioning problem parameters of interest:
00167   //    average_cuts
00168   //    rectilinear_blocks
00169   //    bisection_num_test_cuts (experimental)
00170 
00171   int val = 0;
00172   pe = pl.getEntryPtr("average_cuts");
00173   if (pe)
00174     val = pe->getValue(&val);
00175 
00176   if (val == 1)
00177     params.set(rcb_averageCuts);
00178 
00179   val = 0;
00180   pe = pl.getEntryPtr("rectilinear_blocks");
00181   if (pe)
00182     val = pe->getValue(&val);
00183 
00184   if (val == 1)
00185     params.set(rcb_rectilinearBlocks);
00186 
00187   int numTestCuts = 1;
00188   pe = pl.getEntryPtr("bisection_num_test_cuts");
00189   if (pe)
00190     numTestCuts = pe->getValue(&numTestCuts);
00191 
00193   // From the CoordinateModel we need:
00194   //    coordinate values
00195   //    coordinate weights, if any
00196   //    coordinate global Ids
00197 
00198   env->debug(DETAILED_STATUS, "Accessing coordinate model");
00199   typedef StridedData<lno_t, scalar_t> input_t;
00200 
00201   bool ignoreWeights = params.test(rcb_balanceCount);
00202 
00203   int coordDim = coords->getCoordinateDim();
00204 
00205   int nWeightsPerCoord = 0;
00206   if (!ignoreWeights) nWeightsPerCoord = coords->getNumWeightsPerCoordinate();
00207 
00208   size_t numLocalCoords = coords->getLocalNumCoordinates();
00209   global_size_t numGlobalCoords = coords->getGlobalNumCoordinates();
00210 
00211   ArrayView<const gno_t> gnos;
00212   ArrayView<input_t>     xyz;
00213   ArrayView<input_t>     wgts;
00214 
00215   coords->getCoordinates(gnos, xyz, wgts);
00216 
00217   Array<ArrayRCP<const scalar_t> > values(coordDim);
00218   for (int dim=0; dim < coordDim; dim++){
00219     ArrayRCP<const scalar_t> ar;
00220     xyz[dim].getInputArray(ar);
00221     values[dim] = ar;
00222   }
00223 
00224   env->debug(DETAILED_STATUS, "Storing weights");
00225 
00226   Array<ArrayRCP<const scalar_t> > weights(nWeightsPerCoord);
00227   for (int widx = 0; widx < nWeightsPerCoord; widx++){
00228     ArrayRCP<const scalar_t> ar;
00229     wgts[widx].getInputArray(ar);
00230     weights[widx] = ar;
00231   }
00232 
00233   if (env->doStatus() && (numGlobalCoords < 500)){
00234     std::ostringstream oss;
00235     oss << "Problem: ";
00236     for (size_t i=0; i < numLocalCoords; i++){
00237       oss << gnos[i] << " (";
00238       for (int dim=0; dim < coordDim; dim++)
00239         oss << (xyz[dim])[i] << " ";
00240       oss << ") ";
00241     }
00242 
00243     env->debug(VERBOSE_DETAILED_STATUS, oss.str());
00244   }
00245 
00247   // From the Solution we get part information.
00248   // If the part sizes for a given criteria are not uniform,
00249   // then they are values that sum to 1.0.
00250   env->debug(DETAILED_STATUS, "Getting part info");
00251 
00252   size_t numGlobalParts = solution->getTargetGlobalNumberOfParts();
00253 
00254   int nSizesPerPart = (nWeightsPerCoord ? nWeightsPerCoord : 1);
00255 
00256   Array<bool> uniformParts(nSizesPerPart);
00257   Array<ArrayRCP<scalar_t> > partSizes(nSizesPerPart);
00258 
00259   for (int widx = 0; widx < nSizesPerPart; widx++){
00260     if (solution->criteriaHasUniformPartSizes(widx)){
00261       uniformParts[widx] = true;
00262     }
00263     else{
00264       scalar_t *tmp = new scalar_t [numGlobalParts];
00265       env->localMemoryAssertion(__FILE__, __LINE__, numGlobalParts, tmp) ;
00266     
00267       for (size_t i=0; i < numGlobalParts; i++){
00268         tmp[i] = solution->getCriteriaPartSize(widx, i);
00269       }
00270 
00271       partSizes[widx] = arcp(tmp, 0, numGlobalParts);
00272     }
00273   }
00274 
00275   // It may not be possible to solve the partitioning problem
00276   // if we have multiple weights with part size
00277   // arrays that differ. So let's be aware of this possibility.
00278 
00279   bool multiplePartSizeSpecs = false;
00280 
00281   if (nSizesPerPart > 1){
00282     for (int widx1 = 0; widx1 < nSizesPerPart; widx1++)
00283       for (int widx2 = widx1+1; widx2 < nSizesPerPart; widx2++)
00284         if (!solution->criteriaHaveSamePartSizes(widx1, widx2)){
00285           multiplePartSizeSpecs = true;
00286           break;
00287         }
00288   }
00289   
00290   if (multiplePartSizeSpecs)
00291     params.set(rcb_multiplePartSizeSpecs);
00292 
00294   // Create the distributed data for the algorithm.
00295   //
00296   // It is a multivector containing one vector for each coordinate
00297   // dimension, plus a vector for each weight.
00298 
00299   env->debug(DETAILED_STATUS, "Creating multivec");
00300   typedef Tpetra::Map<lno_t, gno_t, node_t> map_t;
00301   typedef Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t> mvector_t;
00302 
00303   int multiVectorDim = coordDim + nWeightsPerCoord;
00304 
00305   gno_t gnoMin, gnoMax;
00306   coords->getIdentifierMap()->getGnoRange(gnoMin, gnoMax);
00307 
00308   RCP<map_t> map;
00309   try{
00310     map = rcp(new map_t(numGlobalCoords, gnos, gnoMin, comm));
00311   }
00312   Z2_THROW_OUTSIDE_ERROR(*env)
00313 
00314   typedef ArrayView<const scalar_t> coordList_t;
00315 
00316   coordList_t *avList = new coordList_t [multiVectorDim];
00317 
00318   for (int dim=0; dim < coordDim; dim++)
00319     avList[dim] = values[dim].view(0, numLocalCoords);
00320 
00321   for (int widx=0, idx=coordDim; widx < nWeightsPerCoord; widx++)
00322     avList[idx++] = weights[widx].view(0, numLocalCoords);
00323 
00324   ArrayRCP<const ArrayView<const scalar_t> > vectors =
00325     arcp(avList, 0, multiVectorDim);
00326 
00327   RCP<mvector_t> mvector;
00328 
00329   try{
00330     mvector = rcp(new mvector_t(
00331       map, vectors.view(0, multiVectorDim), multiVectorDim));
00332   }
00333   Z2_THROW_OUTSIDE_ERROR(*env)
00334 
00335   env->timerStop(BOTH_TIMERS, "RCB set up");
00336 
00339   // The algorithm
00342 
00343   env->debug(DETAILED_STATUS, "Beginning algorithm");
00344   part_t part0 = 0;
00345   part_t part1 = numGlobalParts-1;
00346   int sanityCheck = numGlobalParts;
00347   int groupSize = comm->getSize();
00348   int rank = comm->getRank();
00349 
00350   long imbalanceReductionFactor(1);
00351   long nparts = numGlobalParts;
00352   while ((nparts >>= 1) != 0) imbalanceReductionFactor++;
00353 
00354   imbalanceTolerance /= imbalanceReductionFactor;
00355 
00356   int iteration = 1;
00357 
00358   env->memory("RCB algorithm set up");
00359 
00360   env->timerStart(MACRO_TIMERS, "Parallel RCB");
00361 
00362   while (part1>part0 && groupSize>1 && numGlobalCoords>0 && sanityCheck--){
00363 
00365     // Which coordinates are left and which are right?
00366 
00367     Array<unsigned char> lrflags(numLocalCoords);
00368     scalar_t cutValue=0;  // TODO eventually save this for user
00369     int cutDimension=0;
00370     scalar_t imbalance=0, weightLeft=0, weightRight=0;
00371     part_t leftHalfNumParts=0;
00372 
00373     env->timerStart(MICRO_TIMERS, "Find cut", iteration, 2);
00374 
00375     try{
00376       determineCut<mvector_t, Adapter>(env, comm, 
00377         params, numTestCuts, imbalanceTolerance,
00378         coordDim, nWeightsPerCoord, mvector, mcnorm, solution, part0, part1,
00379         lrflags.view(0, numLocalCoords), 
00380         cutDimension, cutValue, imbalance, leftHalfNumParts,
00381         weightLeft, weightRight);
00382     }
00383     Z2_FORWARD_EXCEPTIONS
00384 
00385     env->timerStop(MICRO_TIMERS, "Find cut", iteration, 2);
00386 
00387     // Do we have empty left or right halves?
00388 
00389     bool skipLeft = (weightLeft == 0);
00390     bool skipRight = (weightRight == 0);
00391 
00393     // Migrate the multivector of data.
00394 
00395     int leftHalfNumProcs=0;
00396 
00397     env->timerStart(MICRO_TIMERS, "Migrate", iteration, 2);
00398 
00399     try{ // on return mvector has my new data
00400 
00401       migrateData<mvector_t>( env, comm, lrflags.view(0,numLocalCoords), 
00402         mvector, leftHalfNumProcs);
00403     }
00404     Z2_FORWARD_EXCEPTIONS
00405 
00406     env->timerStop(MICRO_TIMERS, "Migrate", iteration, 2);
00407 
00408     env->localBugAssertion(__FILE__, __LINE__, "num procs in half",
00409       leftHalfNumProcs > 0 && leftHalfNumProcs < groupSize,
00410       BASIC_ASSERTION);
00411 
00412     bool inLeftHalf = (rank < leftHalfNumProcs);
00413 
00414     if ((inLeftHalf && skipLeft) || (!inLeftHalf && skipRight)){
00415       groupSize = 1;
00416       numLocalCoords = 0;
00417       continue;
00418     }
00419 
00421     // Divide into two subgroups.
00422 
00423     env->timerStart(MICRO_TIMERS, "Create sub group, sub data", iteration, 2);
00424 
00425     int *ids = NULL;
00426 
00427     if (rank < leftHalfNumProcs){
00428       groupSize = leftHalfNumProcs;
00429       ids = new int [groupSize];
00430       env->localMemoryAssertion(__FILE__, __LINE__, groupSize, ids);
00431       for (int i=0; i < groupSize; i++)
00432         ids[i] = i;
00433       part1 = part0 + leftHalfNumParts - 1;
00434     }
00435     else {
00436       groupSize = comm->getSize() - leftHalfNumProcs;
00437       rank -= leftHalfNumProcs;
00438       ids = new int [groupSize];
00439       env->localMemoryAssertion(__FILE__, __LINE__, groupSize, ids);
00440       for (int i=0; i < groupSize; i++)
00441         ids[i] = i + leftHalfNumProcs;
00442       part0 += leftHalfNumParts;
00443     }
00444 
00445     ArrayView<const int> idView(ids, groupSize);
00446     comm = comm->createSubcommunicator(idView);
00447 
00448     delete [] ids;
00449 
00451     // Create a new multivector for my smaller group.
00452 
00453     ArrayView<const gno_t> gnoList = mvector->getMap()->getNodeElementList();
00454     size_t localSize = mvector->getLocalLength();
00455   
00456     // Tpetra will calculate the globalSize.
00457     size_t globalSize = Teuchos::OrdinalTraits<size_t>::invalid();
00458   
00459     RCP<map_t> subMap;
00460     try{
00461       subMap= rcp(new map_t(globalSize, gnoList, 0, comm));
00462     }
00463     Z2_THROW_OUTSIDE_ERROR(*env)
00464 
00465     coordList_t *avSubList = new coordList_t [multiVectorDim];
00466   
00467     for (int dim=0; dim < multiVectorDim; dim++)
00468       avSubList[dim] = mvector->getData(dim).view(0, localSize);
00469   
00470     ArrayRCP<const ArrayView<const scalar_t> > subVectors =
00471       arcp(avSubList, 0, multiVectorDim);
00472   
00473     RCP<mvector_t> subMvector;
00474   
00475     try{
00476       subMvector = rcp(new mvector_t(
00477         subMap, subVectors.view(0, multiVectorDim), multiVectorDim));
00478     }
00479     Z2_THROW_OUTSIDE_ERROR(*env)
00480 
00481     env->timerStop(MICRO_TIMERS, "Create sub group, sub data", iteration, 2);
00482   
00483     mvector = subMvector;
00484 
00485     numLocalCoords = mvector->getLocalLength();
00486     numGlobalCoords = mvector->getGlobalLength();
00487 
00488     iteration++;
00489 
00490     env->memory("New subgroup data created");
00491   } 
00492 
00493   env->timerStop(MACRO_TIMERS, "Parallel RCB");
00494 
00495   env->localBugAssertion(__FILE__, __LINE__, "partitioning failure", 
00496     sanityCheck, BASIC_ASSERTION);
00497 
00498   ArrayRCP<part_t> partId;
00499 
00500   if (numLocalCoords > 0){
00501     part_t *tmp = new part_t [numLocalCoords];
00502     env->localMemoryAssertion(__FILE__, __LINE__, numLocalCoords, tmp);
00503     partId = arcp(tmp, 0, numLocalCoords, true);
00504   }
00505 
00506   env->memory("Solution array created");
00507 
00508   if ((part1 > part0) && (numLocalCoords > 0)){ // Serial partitioning
00509 
00510     // scalar_t cutValue;   TODO
00511     // int cutDimension;
00512     // scalar_t imbalance;
00513 
00514     env->timerStart(MACRO_TIMERS, "Serial RCB");
00515 
00516     try{
00517       ArrayView<lno_t> emptyIndex;
00518 
00519       serialRCB<mvector_t, Adapter>(env, 1, params,
00520         numTestCuts, imbalanceTolerance,
00521         coordDim, nWeightsPerCoord, mvector, emptyIndex, solution,
00522         part0, part1, partId.view(0,numLocalCoords));
00523     }
00524     Z2_FORWARD_EXCEPTIONS
00525 
00526     env->timerStop(MACRO_TIMERS, "Serial RCB");
00527 
00528   }
00529   else{
00530     for (lno_t i=0; i < partId.size(); i++)
00531       partId[i] = part0;
00532   }
00533 
00535   // Done: update the solution
00536 
00537   ArrayRCP<const gno_t> gnoList = 
00538     arcpFromArrayView(mvector->getMap()->getNodeElementList());
00539 
00540   if (env->getDebugLevel() >= VERBOSE_DETAILED_STATUS && 
00541      (numGlobalCoords < 500)){
00542     std::ostringstream oss;
00543     oss << "Solution: ";
00544     for (gno_t i=0; i < gnoList.size(); i++)
00545       oss << gnoList[i] << " (" << partId[i] << ") ";
00546     
00547     env->debug(VERBOSE_DETAILED_STATUS, oss.str());
00548   }
00549 
00550   solution->setParts(gnoList, partId, false);
00551 #endif // INCLUDE_ZOLTAN2_EXPERIMENTAL
00552 }
00553 
00554 } // namespace Zoltan2
00555 
00556 #endif