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