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