Zoltan2
Zoltan2_AlgRCB_methods.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_METHODS_HPP_
00051 #define _ZOLTAN2_ALGRCB_METHODS_HPP_
00052 
00053 #include <Zoltan2_PartitioningSolution.hpp>
00054 #include <Zoltan2_XpetraTraits.hpp>
00055 #include <Zoltan2_Metric.hpp>
00056 
00057 #include <Tpetra_Vector.hpp>
00058 
00059 #include <sstream>
00060 #include <string>
00061 #include <cmath>
00062 #include <bitset>
00063 
00064 namespace Teuchos{
00065 
00072 template <typename Ordinal, typename T>
00073   class Zoltan2_RCBOperation : public ValueTypeReductionOp<Ordinal,T>
00074 {
00075 private:
00076   Ordinal numSum_, numMin_, numMax_;
00077 
00078 public:
00081   Zoltan2_RCBOperation():numSum_(0), numMin_(0), numMax_(0){}
00082 
00089   Zoltan2_RCBOperation(Ordinal nsum, Ordinal nmin, Ordinal nmax):
00090     numSum_(nsum), numMin_(nmin), numMax_(nmax){}
00091   
00094   void reduce( const Ordinal count, const T inBuffer[], T inoutBuffer[]) const
00095   {
00096     Ordinal next=0;
00097     for (Ordinal i=0; i < numSum_; i++, next++)
00098       inoutBuffer[next] += inBuffer[next];
00099 
00100     for (Ordinal i=0; i < numMin_; i++, next++)
00101       if (inoutBuffer[next] > inBuffer[next])
00102         inoutBuffer[next] = inBuffer[next];
00103 
00104     for (Ordinal i=0; i < numMax_; i++, next++)
00105       if (inoutBuffer[next] < inBuffer[next])
00106         inoutBuffer[next] = inBuffer[next];
00107   }
00108 
00111   ArrayRCP<T> getInitializedBuffer(T minExtreme, T maxExtreme) const
00112   {
00113     Ordinal len = numSum_+numMin_+numMax_;
00114     T *buf = new T [len];
00115     if (!buf)
00116       throw std::bad_alloc();
00117 
00118     Ordinal next=0;
00119     for (Ordinal i=0; i < numSum_; i++)
00120       buf[next++] = 0;
00121     for (Ordinal i=0; i < numMin_; i++)
00122       buf[next++] = minExtreme;
00123     for (Ordinal i=0; i < numMax_; i++)
00124       buf[next++] = maxExtreme;
00125 
00126     ArrayRCP<T> returnBuf = arcp(buf, 0, len, true);
00127     return returnBuf;
00128   }
00129 };
00130 } // namespace Teuchos
00131 
00132 namespace Zoltan2{
00133 
00136 enum rcbParams {
00137   rcb_balanceCount,            
00138   rcb_balanceWeight,          
00139   rcb_minTotalWeight,      
00140   rcb_minMaximumWeight,  
00141   rcb_balanceTotalMaximum, 
00142   rcb_averageCuts,          
00143   rcb_rectilinearBlocks,    
00144   rcb_multiplePartSizeSpecs,  
00145   NUM_RCB_PARAMS
00146 };
00147 
00155 enum leftRightFlag{
00156   unsetFlag = 0xfd,       
00157   leftFlag = 0xfe,        
00158   rightFlag = 0xff        
00159 };
00160 
00161 // TODO: RCB has several helper methods.
00162 // Do we want a convention for naming algorithm methods?  Do
00163 // we want sub-namespaces for algorithms? 
00164 
00178 template <typename Adapter>
00179   void getFractionLeft(
00180     const RCP<const Environment> &env,
00181     partId_t part0,
00182     partId_t part1,
00183     const RCP<PartitioningSolution<Adapter> > &solution,
00184     ArrayRCP<double> &fractionLeft,
00185     partId_t &numPartsLeftHalf)
00186 {
00187   env->debug(DETAILED_STATUS, "Entering getFractionLeft");
00188   partId_t numParts = part1 - part0 + 1;
00189   // TODO In LDRD, substitute call to machine model for
00190   // TODO computation of numPartsLeftHalf below.
00191   numPartsLeftHalf = numParts / 2;
00192   partId_t left0 = part0;   // First part in left half
00193   partId_t left1 = left0 + numPartsLeftHalf - 1;  // Last part in left half
00194   partId_t right0 = left1 + 1;  // First part in right half
00195   partId_t right1 = part1;  // Last part in right half
00196 
00197   int nVecs = solution->getNumberOfCriteria();
00198   fractionLeft = arcp(new double [nVecs], 0, nVecs);
00199 
00200   for (int widx=0; widx<nVecs; widx++){
00201     if (solution->criteriaHasUniformPartSizes(widx)){
00202       fractionLeft[widx] = double(numPartsLeftHalf) / double(numParts);
00203     }
00204     else{
00205       fractionLeft[widx] = 0;
00206       for(int partId=left0; partId <= left1; partId++){
00207         fractionLeft[widx] += solution->getCriteriaPartSize(widx, partId);
00208       }
00209       double total = fractionLeft[widx];
00210       for(int partId=right0; partId <= right1; partId++){
00211         total += solution->getCriteriaPartSize(widx, partId);
00212       }
00213       fractionLeft[widx] /= total;
00214     }
00215   }
00216   env->debug(DETAILED_STATUS, "Exiting getFractionLeft");
00217 }
00218 
00237 template <typename mvector_t>
00238   void getCutDimension(
00239     const RCP<const Environment> &env,
00240     const RCP<Comm<int> > &comm,
00241     int coordDim,
00242     const RCP<mvector_t> &vectors,
00243     ArrayView<typename mvector_t::local_ordinal_type> index,
00244     int &dimension,                        // output
00245     typename mvector_t::scalar_type &minCoord,      // output
00246     typename mvector_t::scalar_type &maxCoord)      // output
00247 {
00248   env->debug(DETAILED_STATUS, "Entering getCutDimension");
00249   typedef typename mvector_t::scalar_type scalar_t;
00250   typedef typename mvector_t::local_ordinal_type lno_t;
00251 
00252   int nprocs = comm->getSize();
00253   bool useIndices = index.size() > 0;
00254   lno_t numLocalCoords = 0;
00255 
00256   if (useIndices)
00257     numLocalCoords = index.size();
00258   else
00259     numLocalCoords = vectors->getLocalLength();
00260 
00261   Array<scalar_t> spans(coordDim*2);
00262   int next = 0;
00263 
00264   if (useIndices){
00265     if (numLocalCoords > 0){
00266       for (int dim=0; dim < coordDim; dim++){
00267         const scalar_t *val = vectors->getData(dim).getRawPtr();
00268         scalar_t v = val[index[0]];
00269         scalar_t min = v;
00270         scalar_t max = v;
00271         for (lno_t i=1; i < numLocalCoords; i++){
00272           v = val[index[i]];
00273           if (v > max)
00274             max = v;
00275           else if (v < min)
00276             min = v;
00277         }
00278   
00279         spans[next++] = min;
00280         spans[next++] = max * -1.0;
00281       }
00282     }
00283     else{
00284       for (int dim=0; dim < coordDim; dim++){
00285         spans[next++] = std::numeric_limits<scalar_t>::max();
00286         spans[next++] = std::numeric_limits<scalar_t>::min() * -1.0;
00287       }
00288     }
00289   }
00290   else{
00291     for (int dim=0; dim < coordDim; dim++){
00292       const scalar_t *val = vectors->getData(dim).getRawPtr();
00293       std::pair<scalar_t, scalar_t> minMax = 
00294         z2LocalMinMax<scalar_t>(val, numLocalCoords);
00295       spans[next++] = minMax.first;
00296       spans[next++] = minMax.second * -1.0;
00297     }
00298   }
00299 
00300   if (nprocs > 1){
00301     Array<scalar_t> globalValues(coordDim*2);
00302     try{
00303       reduceAll<int, scalar_t>(*comm, Teuchos::REDUCE_MIN, coordDim*2,
00304         spans.getRawPtr(), globalValues.getRawPtr());
00305     }
00306     Z2_THROW_OUTSIDE_ERROR(*env)
00307 
00308     for (int i=0; i < coordDim*2; i++)
00309       spans[i] = globalValues[i];
00310   }
00311 
00312   scalar_t maxSpan = 0;
00313   dimension = -1;
00314   next = 0;
00315 
00316   for (int dim=0; dim < coordDim; dim++){
00317     scalar_t min = spans[next++];
00318     scalar_t max = spans[next++] * -1.0;
00319     scalar_t newSpan = max - min;
00320     if (newSpan > maxSpan){
00321       maxSpan = newSpan;
00322       dimension = dim;
00323       minCoord = min;
00324       maxCoord = max;
00325     }
00326   }
00327   env->debug(DETAILED_STATUS, "Exiting getCutDimension");
00328 }
00329 
00346 template <typename mvector_t>
00347   void migrateData(
00348     const RCP<const Environment> &env, const RCP<Comm<int> > &comm,
00349     ArrayView<unsigned char> lrflags,
00350     RCP<mvector_t> &vectors,    // on return is the new data
00351     int &leftNumProcs)          // on return is num procs with left data
00352 {
00353   env->debug(DETAILED_STATUS, "Entering migrateData");
00354   typedef typename mvector_t::scalar_type scalar_t;
00355   typedef typename mvector_t::local_ordinal_type lno_t;
00356   typedef typename mvector_t::global_ordinal_type gno_t;
00357 
00358   int nprocs = comm->getSize();
00359   size_t nobj = vectors->getLocalLength();
00360   size_t nGlobalObj = vectors->getGlobalLength();
00361 
00362   env->localBugAssertion(__FILE__, __LINE__, "migrateData input", 
00363     nprocs>1 && size_t(lrflags.size())==nobj, DEBUG_MODE_ASSERTION);
00364 
00365   gno_t myNumLeft= 0, numLeft;
00366   for (size_t i=0; i < nobj; i++)
00367     if (lrflags[i] == leftFlag)
00368       myNumLeft++;
00369 
00370   try{
00371     reduceAll<int, gno_t>(
00372       *comm, Teuchos::REDUCE_SUM, 1, &myNumLeft, &numLeft);
00373   }
00374   Z2_THROW_OUTSIDE_ERROR(*env)
00375     
00376   scalar_t leftFraction = scalar_t(numLeft)/scalar_t(nGlobalObj);
00377   leftNumProcs = static_cast<int>(floor((nprocs*leftFraction)+.5));
00378 
00379   if (leftNumProcs == 0)
00380     leftNumProcs++;
00381   else if (leftNumProcs == nprocs)
00382     leftNumProcs--;
00383 
00385   // Get a list of my new global numbers.
00386 
00387   Array<int> sendCount(nprocs, 0);
00388   Array<int> recvCount(nprocs, 0);
00389   Array<gno_t> sendBuf(nobj, 0);
00390 
00391   if (nobj > 0){
00392     int *procId = new int [nobj];
00393     env->localMemoryAssertion(__FILE__, __LINE__, nobj, procId) ;
00394     int leftProc0 = 0;
00395     int rightProc0 = leftProc0 + leftNumProcs;
00396   
00397     int nextLeftProc = leftProc0;
00398     int nextRightProc = rightProc0;
00399     int *p = procId;
00400   
00401     for (size_t i=0; i < nobj; i++){
00402       if (lrflags[i] == leftFlag){
00403         if (nextLeftProc == rightProc0)
00404           nextLeftProc = leftProc0;
00405   
00406         sendCount[nextLeftProc]++;
00407         *p++ = nextLeftProc++;
00408       }
00409       else{
00410         if (nextRightProc == nprocs)
00411           nextRightProc = rightProc0;
00412         sendCount[nextRightProc]++;
00413   
00414         *p++ = nextRightProc++;
00415       }
00416     }
00417   
00418     gno_t *sendOffset = new gno_t [nprocs];
00419     env->localMemoryAssertion(__FILE__, __LINE__, nprocs, sendOffset) ;
00420     sendOffset[0] = 0;
00421     for (int i=0; i < nprocs-1; i++)
00422       sendOffset[i+1] = sendOffset[i] + sendCount[i];
00423 
00424     ArrayView<const gno_t> gnoList = vectors->getMap()->getNodeElementList();
00425 
00426     for (size_t i=0; i < nobj; i++){
00427       int proc = procId[i];
00428       lno_t offset = sendOffset[proc]++;
00429 
00430       sendBuf[offset] = gnoList[i];
00431     }
00432 
00433     delete [] sendOffset;
00434     delete [] procId;
00435   }
00436 
00437   ArrayRCP<gno_t> recvBuf;
00438 
00439   try{
00440     AlltoAllv<gno_t>(*comm, *env,
00441       sendBuf(), sendCount(),
00442       recvBuf, recvCount());
00443   }
00444   Z2_FORWARD_EXCEPTIONS
00445 
00446   sendCount.clear();
00447   sendBuf.clear();
00448 
00450   // Migrate the multivector of data.
00451 
00452   gno_t numMyNewGnos = 0;
00453   for (int i=0; i < nprocs; i++)
00454     numMyNewGnos += recvCount[i];
00455 
00456   recvCount.clear();
00457 
00458   RCP<const mvector_t> newMultiVector;
00459   RCP<const mvector_t> constInput = rcp_const_cast<const mvector_t>(vectors);
00460 
00461   try{
00462     newMultiVector = XpetraTraits<mvector_t>::doMigration(
00463       constInput, numMyNewGnos, recvBuf.getRawPtr());
00464   }
00465   Z2_FORWARD_EXCEPTIONS
00466 
00467   vectors = rcp_const_cast<mvector_t>(newMultiVector);
00468   env->memory("Former problem data replaced with new data");
00469   env->debug(DETAILED_STATUS, "Exiting migrateData");
00470 }
00471 
00472 template <typename lno_t, typename scalar_t>
00473   scalar_t getCoordWeight(lno_t id, 
00474     multiCriteriaNorm mcnorm, 
00475     ArrayView<StridedData<lno_t, scalar_t> > weights)
00476 {
00477   scalar_t wgtValue = 1.0;
00478   size_t nWeightsPerCoord = weights.size();
00479     
00480   if (nWeightsPerCoord > 1){
00481     Array<scalar_t> coordWeights(nWeightsPerCoord, 1.0);
00482     for (size_t widx=0; widx < nWeightsPerCoord; widx++){
00483       coordWeights[widx] = weights[widx][id];
00484     }
00485 
00486     wgtValue = normedWeight<scalar_t>(coordWeights.view(0,nWeightsPerCoord), mcnorm);
00487   }
00488   else if (nWeightsPerCoord > 0){
00489     wgtValue = weights[0][id];
00490   }
00491 
00492   return wgtValue;
00493 }
00494 
00512 template <typename scalar_t>
00513   bool emptyPartsCheck( const RCP<const Environment> &env,
00514    const ArrayView<double> fractionLeft, 
00515    scalar_t minCoord, scalar_t maxCoord,
00516    leftRightFlag &lrf,
00517    scalar_t &cutValue)
00518 {
00519   env->debug(DETAILED_STATUS, "Entering emptyPartsCheck");
00520   // initialize return values
00521   lrf = leftFlag;
00522   cutValue = 0.0;
00523 
00524   size_t nVecs = fractionLeft.size();
00525   int numEmptyRight = 0, numEmptyLeft = 0;
00526 
00527   for (size_t i=0; i < nVecs; i++){
00528     if (fractionLeft[i] == 0.0)
00529       numEmptyLeft++;
00530     else if (fractionLeft[i] == 1.0)
00531       numEmptyRight++;
00532   }
00533 
00534   if (numEmptyRight + numEmptyLeft == 0)
00535     return false;
00536 
00537   env->localInputAssertion(__FILE__, __LINE__,
00538     "partitioning not solvable - conflicting part size criteria",
00539     (numEmptyRight==0) || (numEmptyLeft==0), BASIC_ASSERTION);
00540 
00541   if (numEmptyLeft > 0){
00542     lrf = leftFlag;
00543     cutValue = minCoord;
00544   }
00545   else{
00546     lrf = rightFlag;
00547     cutValue = maxCoord;
00548   }
00549 
00550   env->debug(DETAILED_STATUS, "Exiting emptyPartsCheck");
00551   return true;
00552 }
00553 
00566 template <typename lno_t, typename gno_t, typename scalar_t>
00567   void testCoordinatesOnRightBoundary(
00568     const RCP<const Environment> &env,
00569     const RCP<Comm<int> > &comm, 
00570     scalar_t totalWeightLeft,
00571     scalar_t targetLeftScalar,
00572     int cutLocation,
00573     ArrayView<scalar_t> localSums,
00574     ArrayView<scalar_t> globalSums,
00575     ArrayView<lno_t> index,
00576     ArrayView<StridedData<lno_t, scalar_t> > weights,
00577     multiCriteriaNorm mcnorm,
00578     ArrayView<unsigned char> lrFlags,
00579     scalar_t &globalWeightMovedRight)   // output
00580 {
00581   env->debug(DETAILED_STATUS, "Entering testCoordinatesOnRightBoundary");
00582   int nprocs = comm->getSize();
00583   int rank = comm->getRank();
00584 
00585   globalWeightMovedRight = 0.0;
00586 
00587   scalar_t localBoundarySum = localSums[cutLocation];
00588 
00589   scalar_t total = totalWeightLeft;
00590   for (int i=0; i <= cutLocation; i++)
00591     total += globalSums[i];
00592 
00593   scalar_t totalMoveRight = total - targetLeftScalar;
00594   scalar_t localMoveRight = localBoundarySum;
00595   scalar_t actualWeightMovedRight = 0.0;
00596   
00597   Array<scalar_t> scansum(nprocs+1, 0.0);
00598   Teuchos::gatherAll<int, scalar_t>(*comm, 1, 
00599     &localBoundarySum, nprocs, scansum.getRawPtr()+1);
00600   for (int i=2; i<=nprocs; i++)
00601     scansum[i] += scansum[i-1];
00602 
00603   if (localBoundarySum > 0.0){
00604 
00605     scalar_t sumMine = scansum[rank]; // sum of ranks preceding me
00606     scalar_t diff = scansum[nprocs] - sumMine;
00607     localMoveRight = 0;
00608 
00609     if (diff <= totalMoveRight)
00610       localMoveRight = localBoundarySum;  // all
00611     else{
00612       scalar_t leftPart = diff - totalMoveRight;
00613       if (leftPart < localBoundarySum)
00614         localMoveRight = localBoundarySum - leftPart;
00615     }
00616   }
00617 
00618   if (localMoveRight > 0.0){
00619 
00620     bool moveAll =  (localMoveRight >= localBoundarySum);
00621 
00622     if (moveAll){
00623       actualWeightMovedRight = localBoundarySum;
00624       for (lno_t i=0; i < lrFlags.size(); i++){
00625         if (lrFlags[i] == cutLocation){
00626           lrFlags[i] = rightFlag;
00627         }
00628       }
00629     }
00630     else{
00631       int nWeightsPerCoord = weights.size();
00632       int useIndices = index.size() > 0;
00633 
00634       for (lno_t i=0; i < lrFlags.size(); i++){
00635         if (lrFlags[i] == cutLocation){
00636           lrFlags[i] = rightFlag;
00637 
00638           lno_t idx = (useIndices ? index[i] : i);
00639 
00640           actualWeightMovedRight += getCoordWeight<lno_t, scalar_t>(
00641             idx, mcnorm, weights.view(0, nWeightsPerCoord));
00642 
00643           if (actualWeightMovedRight >= localMoveRight)
00644             break;
00645         }
00646       } // next coordinate
00647     }
00648   }
00649 
00650   try{
00651     reduceAll<int, scalar_t>(
00652       *comm, Teuchos::REDUCE_SUM, 1, &actualWeightMovedRight,
00653       &globalWeightMovedRight);
00654   }
00655   Z2_THROW_OUTSIDE_ERROR(*env)
00656 
00657   env->debug(DETAILED_STATUS, "Exiting testCoordinatesOnRightBoundary");
00658   return;
00659 }
00660 
00707 template <typename mvector_t>
00708   void BSPfindCut(
00709     const RCP<const Environment> &env,
00710     const RCP<Comm<int> > &comm,
00711     const std::bitset<NUM_RCB_PARAMS> &params,
00712     int numTestCuts,
00713     typename mvector_t::scalar_type tolerance,
00714     int cutDim,
00715     int coordDim,
00716     int nWeightsPerCoord,
00717     const RCP<mvector_t> &vectors,
00718     ArrayView<typename mvector_t::local_ordinal_type> index,
00719     ArrayView<double> fractionLeft,
00720     typename mvector_t::scalar_type coordGlobalMin,
00721     typename mvector_t::scalar_type coordGlobalMax,
00722     typename mvector_t::scalar_type &cutValue,         // output
00723     ArrayView<unsigned char> lrFlags,                  // output
00724     typename mvector_t::scalar_type &totalWeightLeft,  // output
00725     typename mvector_t::scalar_type &totalWeightRight, // output
00726     typename mvector_t::local_ordinal_type &localCountLeft, // output
00727     typename mvector_t::scalar_type &imbalance)        // output
00728 {
00729   env->debug(DETAILED_STATUS, "Entering BSPfindCut");
00730 
00731   // initialize output
00732   bool useIndices = index.size() > 0;
00733   int numAllCoords = vectors->getLocalLength();
00734   int numCoords = (useIndices ? index.size() : numAllCoords);
00735   cutValue = totalWeightLeft = totalWeightRight = imbalance = 0.0;
00736   localCountLeft = 0;
00737 
00738   for (int i=0; i < numCoords; i++)
00739     lrFlags[i] = unsetFlag;
00740 
00741   typedef typename mvector_t::scalar_type scalar_t;
00742   typedef typename mvector_t::local_ordinal_type lno_t;
00743   typedef typename mvector_t::global_ordinal_type gno_t;
00744   typedef StridedData<lno_t, scalar_t> input_t;
00745 
00746   bool multiplePartSizeSpecs = params.test(rcb_multiplePartSizeSpecs);
00747   bool rectilinearBlocks = params.test(rcb_rectilinearBlocks);
00748   bool averageCuts = params.test(rcb_averageCuts);
00749 
00750   // A coordinate is considered to be on a cut if it is within
00751   // this distance of the cut.
00752 
00753   double epsilon = (coordGlobalMax - coordGlobalMin) * 10e-9;
00754 
00755   // Find the coordinate values and weights.
00756 
00757   int nVecs = fractionLeft.size();
00758 
00759   if (env->getDebugLevel() >= DETAILED_STATUS){
00760     std::ostringstream info;
00761     info << "Num weights " << nWeightsPerCoord << ", Fraction left:";
00762     for (int i=0; i < nVecs; i++)
00763       info << " " << fractionLeft[i];
00764     info << endl << "Dimension " << cutDim << " [";
00765     info << coordGlobalMin << ", " << coordGlobalMax << "]";
00766     info << endl << "# test cuts " << numTestCuts;
00767     info << ", tolerance " << tolerance << std::endl;
00768     env->debug(DETAILED_STATUS, info.str());
00769   }
00770 
00771   const scalar_t *coordValue = vectors->getData(cutDim).getRawPtr();
00772 
00773   input_t *wgtinfo = new input_t [nWeightsPerCoord];
00774   env->localMemoryAssertion(__FILE__, __LINE__, nWeightsPerCoord, wgtinfo);
00775   ArrayRCP<input_t> weight(wgtinfo, 0, nWeightsPerCoord, true);
00776 
00777   for (int widx = 0, cidx=coordDim; widx < nWeightsPerCoord; widx++){
00778     weight[widx] = input_t(vectors->getData(cidx++), 1);
00779   }
00780 
00781   // Multicriteria norm
00782 
00783   multiCriteriaNorm mcnorm = normBalanceTotalMaximum;
00784 
00785   if (params.test(rcb_minMaximumWeight))
00786     mcnorm = normMinimizeMaximumWeight;
00787   else if (params.test(rcb_minTotalWeight))
00788     mcnorm = normMinimizeTotalWeight;
00789   
00790   // Goal is globally find one cut that comes close to leaving
00791   // partSizeLeft*totalWeight on the left side.
00792 
00793   Epetra_SerialDenseVector partSizeLeft(View, fractionLeft.getRawPtr(), nVecs);
00794   
00795   // Where do we make the first test cuts?
00796   //
00797   //   min     1     2     3     max
00798   //
00799   // Space is [min, max].  
00800   // If there are three test cuts: 1, 2 and 3.
00801   //   4 regions: [min,1] [1,2] [2,3] [3,max].
00802   //   5 boundaries: min, 1, 2, 3, and max.
00803 
00804   int numRegions = numTestCuts + 1;
00805   int numBoundaries = numTestCuts + 2;
00806   int endBoundary = numBoundaries - 1;
00807   std::vector<scalar_t> boundaries(numBoundaries);
00808   std::vector<scalar_t> searchBoundaries(numBoundaries);
00809 
00810   int numSums = numBoundaries+numRegions;
00811 
00812   Teuchos::Zoltan2_RCBOperation<int, scalar_t> reductionOp(
00813     numSums,      // number of sums
00814     numRegions,   // number of mins
00815     numRegions);  // number of maxes
00816 
00817   typename std::vector<scalar_t>::iterator foundCut;
00818 
00819   bool done=false;
00820   bool fail=false;
00821   scalar_t min = coordGlobalMin;
00822   scalar_t max = coordGlobalMax;
00823   lno_t numRemaining = numCoords;
00824   lno_t prevNumRemaining = numCoords;
00825   size_t numGlobalPoints = vectors->getGlobalLength();
00826   size_t sanityCheck = numGlobalPoints;
00827 
00828   double totalWeight = 0;
00829   double targetLeftScalar = 0;
00830   double targetLeftNorm = 0;
00831   Epetra_SerialDenseVector targetLeftVector(nVecs);
00832 
00833   while (!done && !fail && sanityCheck--){
00834 
00835     // Create regions into which coordinates will be placed.
00836     scalar_t diff = (max - min) / numRegions;
00837     boundaries[0] = min;
00838     boundaries[endBoundary] = max;
00839     for (int i=0; i < endBoundary; i++){
00840       searchBoundaries[i+1] = boundaries[i+1] = boundaries[i] + diff;
00841     }
00842 
00843     // Move ends slightly so we catch points on boundary.
00844     searchBoundaries[0] = min - epsilon;
00845     searchBoundaries[endBoundary] = max + epsilon;
00846 
00847     // Save region and boundary sums, and region min and max.
00848  
00849     ArrayRCP<scalar_t> localSums = reductionOp.getInitializedBuffer(
00850       searchBoundaries[endBoundary], searchBoundaries[0]);
00851  
00852     ArrayRCP<scalar_t> globalSums = reductionOp.getInitializedBuffer(
00853       searchBoundaries[endBoundary], searchBoundaries[0]);
00854 
00855     scalar_t *sums = localSums.getRawPtr();
00856     scalar_t *regionMin = sums + numSums;
00857     scalar_t *regionMax = regionMin + numRegions;
00858 
00859     if (numRemaining > 0){
00860 
00861       // Assign each of my points to a region.
00862       // lower_bound() finds the first cut f, such that f >= coordValue[i].
00863       // So for now, objects that are on the cut boundary go into the
00864       // region on the "left" side.
00865 
00866       for (lno_t i=0; i < numCoords; i++){
00867   
00868         if (lrFlags[i] != unsetFlag)
00869           continue;
00870 
00871         int inRegion = 0;
00872         int idx = (useIndices ? index[i] : i);
00873         scalar_t value = coordValue[idx];
00874 
00875         if (numRegions > 2){
00876        
00877           foundCut = std::lower_bound(
00878             searchBoundaries.begin(), searchBoundaries.end(), value);
00879           
00880           env->localBugAssertion(__FILE__, __LINE__, "search cuts", 
00881             foundCut != searchBoundaries.end(), BASIC_ASSERTION);
00882 
00883           inRegion = foundCut - searchBoundaries.begin() - 1;        
00884         }
00885         else{
00886           if (value <= boundaries[1])
00887             inRegion = 0;
00888           else
00889             inRegion = 1;
00890         }
00891 
00892         int sumIdx = 1 + inRegion*2;
00893 
00894         if (value >= boundaries[inRegion+1]-epsilon){  
00895           // "on" right boundary of this region
00896           sumIdx++;
00897         }
00898         else if (inRegion==0 && (value < min + epsilon)){
00899           // in region 0 but far left boundary
00900           sumIdx--;
00901         }
00902 
00903         lrFlags[i] = (unsigned char)sumIdx;
00904 
00905         if (value < regionMin[inRegion])
00906           regionMin[inRegion] = value;
00907         if (value > regionMax[inRegion])
00908           regionMax[inRegion] = value;
00909 
00910         sums[sumIdx] += getCoordWeight<lno_t, scalar_t>(idx, 
00911               mcnorm, weight.view(0, nWeightsPerCoord));
00912 
00913       } // next coord
00914     }
00915 
00916     try{
00917       reduceAll<int, scalar_t>(*comm, reductionOp,
00918         localSums.size(), localSums.getRawPtr(), globalSums.getRawPtr());
00919     }
00920     Z2_THROW_OUTSIDE_ERROR(*env)
00921 
00922     sums = globalSums.getRawPtr();
00923     regionMin = sums + numSums;
00924     regionMax = regionMin + numRegions;
00925 
00926     if (env->getDebugLevel() >= DETAILED_STATUS){
00927       std::ostringstream info;
00928       info << "  Region " << min << " - " << max << std::endl;
00929       info << "  Remaining to classify: " << numRemaining << std::endl;
00930       info << "  Boundaries: ";
00931       for (int i=0; i < numBoundaries; i++)
00932         info << boundaries[i] << " ";
00933       info << endl << "  For searching: ";
00934       for (int i=0; i < numBoundaries; i++)
00935         info << searchBoundaries[i] << " ";
00936       info << endl << "  Global sums: ";
00937       double sumTotal=0;
00938       for (int i=0; i < numSums; i++){
00939         sumTotal += sums[i];
00940         info << sums[i] << " ";
00941       }
00942       info << " total: " << sumTotal << endl;
00943       env->debug(DETAILED_STATUS, info.str());
00944     }
00945 
00946     if (totalWeight == 0){   // first time through only
00947 
00948       for (int i=0; i < numSums; i++)
00949         totalWeight += sums[i];
00950 
00951       partSizeLeft.Scale(totalWeight);
00952       targetLeftVector = partSizeLeft;
00953 
00954       targetLeftScalar = targetLeftVector[0];
00955       targetLeftNorm = targetLeftVector.Norm2();
00956       totalWeightLeft = 0;
00957       totalWeightRight = 0;
00958     }
00959 
00960     int cutLocation=0;
00961     scalar_t testDiff=0, prevTestDiff=0, target=0;
00962 
00963     if (multiplePartSizeSpecs){
00964       // more complex: if we have multiple weights
00965       //   and the part sizes requested
00966       //   for each weight index differ, then we may not
00967       //   be able to reach the imbalance tolerance.
00968       //
00969       // TODO: discuss how to (or whether to) handle this case.
00970       // 
00971       // For now we look at this imbalance:
00972       // 
00973       //    |target - actual|^2 / |target|^2
00974   
00975       target = targetLeftNorm;
00976       Epetra_SerialDenseVector testVec(nVecs);
00977       for (int i=0; i < nVecs; i++)
00978         testVec[i] = totalWeightLeft;
00979       Epetra_SerialDenseVector diffVec = testVec;
00980       diffVec.Scale(-1.0);
00981       diffVec += targetLeftVector;
00982 
00983       testDiff = diffVec.Norm2(); // imbalance numerator
00984       prevTestDiff = testDiff;
00985       cutLocation= -1;
00986 
00987       while (++cutLocation< numSums){
00988 
00989         for (int i=0; i < nVecs; i++)
00990           testVec[i] += sums[cutLocation];
00991   
00992         diffVec = testVec;
00993         diffVec.Scale(-1.0);
00994         diffVec += targetLeftVector;
00995   
00996         testDiff = diffVec.Norm2();
00997         
00998         if (testDiff >= target)
00999           break;
01000 
01001         prevTestDiff = testDiff;
01002       }
01003     }
01004     else{    // the part sizes for each weight index are the same
01005   
01006       target = targetLeftScalar;
01007       testDiff = totalWeightLeft; 
01008       prevTestDiff = testDiff;
01009       cutLocation = -1;
01010   
01011       while (++cutLocation < numSums){
01012      
01013         testDiff += sums[cutLocation];
01014         if (testDiff >= target)
01015           break;
01016   
01017         prevTestDiff = testDiff;
01018       }
01019     }
01020 
01021     scalar_t diffLeftCut = target - prevTestDiff;
01022     scalar_t diffRightCut = testDiff - target;
01023 
01024     if (diffLeftCut < diffRightCut){
01025       imbalance = diffLeftCut / target;
01026       if (imbalance <= tolerance){
01027         env->debug(DETAILED_STATUS, "  Done, tolerance met");
01028         done = true;
01029         cutLocation--;
01030       }
01031     } 
01032     else{
01033       imbalance = diffRightCut / target;
01034       if (imbalance <= tolerance){
01035         env->debug(DETAILED_STATUS, "  Done, tolerance met");
01036         done = true;
01037      }
01038     }
01039 
01040     bool cutLocIsRegion = (cutLocation % 2 == 1);
01041     bool cutLocIsBoundary = !cutLocIsRegion;
01042 
01043     if (env->getDebugLevel() >= DETAILED_STATUS){
01044       std::ostringstream info;
01045       info << "  Best cut location: " << cutLocation;
01046       if (cutLocIsRegion) info << " just after a region." << std::endl;
01047       else info << " just after a boundary." << std::endl;
01048       env->debug(DETAILED_STATUS, info.str());
01049     }
01050 
01051     if (!done && cutLocIsBoundary){
01052 
01053       done = true;    // can not divide space any more
01054 
01055       env->debug(DETAILED_STATUS, "  Done, cutting at a region boundary");
01056 
01057       if (rectilinearBlocks){
01058         // Can not divide boundary points into two
01059         // different regions to achieve balance.
01060         fail = true;
01061       }
01062       else {
01063         // Maybe by moving some of the points right, we can
01064         // obtain a better balance.  If so, the lrFlag for
01065         // the points moved right will be updated here.
01066 
01067         scalar_t globalWeightMovedRight(0);
01068 
01069         testCoordinatesOnRightBoundary<lno_t, gno_t, scalar_t>(
01070           env, comm, totalWeightLeft, targetLeftScalar, cutLocation, 
01071           localSums.view(0, numSums), globalSums.view(0, numSums),
01072           index, weight.view(0, nWeightsPerCoord), mcnorm, lrFlags,
01073           globalWeightMovedRight);
01074 
01075         scalar_t newSum = testDiff - globalWeightMovedRight;
01076         globalSums[cutLocation] -= globalWeightMovedRight;
01077 
01078         if (newSum > target)
01079           imbalance = (target - newSum) / target;
01080         else
01081           imbalance = (newSum - target) / target;
01082 
01083         if (imbalance > tolerance)
01084           fail = true;
01085       }
01086     }
01087 
01088     int rightmostLeftNum=0, leftmostRightNum=0;
01089 
01090     if (!done){
01091       // Best cut is following a region.  Narrow down the boundaries.
01092 
01093       int regionNumber = (cutLocation - 1) / 2;
01094 
01095       min = regionMin[regionNumber];
01096       max = regionMax[regionNumber];
01097 
01098       rightmostLeftNum = cutLocation - 1;
01099       leftmostRightNum = cutLocation + 1;
01100     }
01101     else{
01102       rightmostLeftNum = cutLocation;
01103       leftmostRightNum = cutLocation + 1;
01104 
01105       if (cutLocIsRegion && averageCuts){
01106         int regionNumber = (cutLocation - 1) / 2;
01107         cutValue = (
01108           boundaries[regionNumber+1] + // boundary to right
01109           regionMax[regionNumber])     // rightmost point in region
01110           / 2.0;
01111       }
01112     }
01113 
01114     for (int i=0; i <= rightmostLeftNum; i++){
01115       totalWeightLeft += globalSums[i];
01116     }
01117 
01118     for (lno_t i=0; i < numCoords; i++){
01119       if (lrFlags[i] != leftFlag && lrFlags[i] != rightFlag){
01120         if (lrFlags[i] <= rightmostLeftNum){
01121           lrFlags[i] = leftFlag;
01122           numRemaining--;
01123           localCountLeft++;
01124         }
01125         else if (lrFlags[i] >= leftmostRightNum){
01126           lrFlags[i] = rightFlag;
01127           numRemaining--;
01128         }
01129         else{
01130           lrFlags[i] = unsetFlag;   // still to be determined
01131         }
01132       }
01133     }
01134 
01135     if (env->getDebugLevel() >= VERBOSE_DETAILED_STATUS && numCoords < 100){
01136       // For large numCoords, building this message
01137       // takes an extraordinarily long time.
01138       std::ostringstream ossLeft;
01139       std::ostringstream ossRight;
01140       ossLeft << "left: ";
01141       ossRight << "right: ";
01142       for (lno_t i=0; i < numCoords; i++){
01143         if (lrFlags[i] == unsetFlag)
01144           continue;
01145         lno_t idx = (useIndices ? index[i] : i);
01146         scalar_t val = coordValue[idx];
01147         if (lrFlags[i] == leftFlag)
01148           ossLeft << val << " ";
01149         else if (lrFlags[i] == rightFlag)
01150           ossRight << val << " ";
01151         else 
01152           env->localBugAssertion(__FILE__, __LINE__, 
01153             "left/right flags", false, BASIC_ASSERTION);
01154       }
01155       std::ostringstream msg;
01156       msg << ossLeft.str() << endl << ossRight.str() << std::endl;
01157       env->debug(VERBOSE_DETAILED_STATUS, msg.str());
01158     }
01159 
01160     prevNumRemaining = numRemaining;
01161   }  // while !done
01162 
01163   totalWeightRight = totalWeight - totalWeightLeft;
01164 
01165   if (fail)
01166     env->debug(BASIC_STATUS, "Warning: tolerance not achieved in sub-step");
01167 
01168   // TODO: If fail the warn that tolerance was not met.
01169 
01170   env->globalInputAssertion(__FILE__, __LINE__, "partitioning not solvable",
01171     done, DEBUG_MODE_ASSERTION, comm);
01172 
01173   env->memory("End of bisection");
01174 
01175   if (env->getDebugLevel() >= DETAILED_STATUS){
01176     std::ostringstream oss;
01177     oss << "Exiting BSPfindCut, ";
01178     oss << "# iterations: " << numGlobalPoints - sanityCheck;    
01179     env->debug(DETAILED_STATUS, oss.str());
01180   }
01181 
01182   return;
01183 }
01184 
01211 template <typename mvector_t, typename Adapter>
01212   void determineCut(
01213     const RCP<const Environment> &env,
01214     const RCP<Comm<int> > &comm,
01215     const std::bitset<NUM_RCB_PARAMS> &params,
01216     int numTestCuts,
01217     typename mvector_t::scalar_type tolerance,
01218     int coordDim, 
01219     int nWeightsPerCoord,
01220     const RCP<mvector_t> &vectors,
01221     multiCriteriaNorm mcnorm,
01222     const RCP<PartitioningSolution<Adapter> > &solution,
01223     partId_t part0, 
01224     partId_t part1,
01225     ArrayView<unsigned char> lrflags,  // output
01226     int &cutDimension,                  // output
01227     typename mvector_t::scalar_type &cutValue,   // output
01228     typename mvector_t::scalar_type &imbalance,  // output
01229     partId_t &numPartsLeftHalf,                  // output
01230     typename mvector_t::scalar_type &weightLeftHalf,  // output
01231     typename mvector_t::scalar_type &weightRightHalf  // output
01232     )
01233 {
01234   env->debug(DETAILED_STATUS, "Entering determineCut");
01235   typedef typename mvector_t::scalar_type scalar_t;
01236   typedef typename mvector_t::local_ordinal_type lno_t;
01237   typedef StridedData<lno_t, scalar_t> input_t;
01238 
01239   lno_t numLocalCoords = vectors->getLocalLength();
01240 
01241   // initialize return values
01242   cutDimension = 0;
01243   cutValue = imbalance = weightLeftHalf = weightRightHalf = 0.0;
01244   numPartsLeftHalf = 0;
01245 
01247   // Pick a cut direction.
01248 
01249   scalar_t globalMinCoord, globalMaxCoord;
01250   ArrayView<lno_t> emptyIndex;
01251 
01252   getCutDimension<mvector_t>(env, comm, coordDim, vectors, emptyIndex,
01253     cutDimension, globalMinCoord, globalMaxCoord);
01254 
01256   // Compute part sizes for the two parts.
01257 
01258   ArrayRCP<double> fractionLeft;
01259 
01260   getFractionLeft<Adapter>(env, part0, part1, solution,
01261     fractionLeft, numPartsLeftHalf);
01262   int nVecs = fractionLeft.size();
01263 
01264   // Special case of empty left or empty right.
01265 
01266   leftRightFlag lrf;
01267 
01268   bool emptyParts = emptyPartsCheck(env,
01269     fractionLeft.view(0, nVecs), // input
01270     globalMinCoord, globalMaxCoord,  // input
01271     lrf, cutValue);                  // output
01272 
01273   if (emptyParts){
01274     
01275     for (lno_t i=0; i < lrflags.size(); i++)
01276       lrflags[i] = lrf;
01277 
01278     imbalance = 0.0;                // perfect balance
01279     scalar_t totalWeight = 0.0;
01280 
01281     if (nWeightsPerCoord == 0)
01282       totalWeight = numLocalCoords;
01283     else if (nWeightsPerCoord == 1) {
01284       int wgtidx = coordDim;
01285       const scalar_t *val = vectors->getData(wgtidx).getRawPtr();
01286       for (lno_t i=0; i < numLocalCoords; i++)
01287         totalWeight += val[i];
01288     }
01289     else {  // need to add up total normed weight
01290       Array<input_t> wgts(nWeightsPerCoord);
01291       int wgtidx = coordDim;
01292       for (int i=0; i < nWeightsPerCoord; i++)
01293         wgts[i] = input_t(vectors->getData(wgtidx++), 1);
01294 
01295       partId_t numParts, numNonemptyParts;
01296       ArrayRCP<MetricValues<scalar_t> > metrics;
01297       ArrayRCP<scalar_t> weightSums;
01298     
01299       globalSumsByPart<scalar_t, unsigned char, lno_t>(
01300         env, comm, lrflags, nWeightsPerCoord,
01301         wgts.view(0, nWeightsPerCoord), mcnorm,
01302         numParts, numNonemptyParts, metrics, weightSums);
01303 
01304       totalWeight = weightSums[1];  // sum normed weight
01305     }
01306 
01307     if (lrf == leftFlag){
01308       numPartsLeftHalf = part1 - part0 + 1;
01309       weightLeftHalf = totalWeight;
01310       weightRightHalf = 0;
01311     }
01312     else{
01313       numPartsLeftHalf = 0;
01314       weightRightHalf = totalWeight;
01315       weightLeftHalf = 0;
01316     }
01317 
01318     env->debug(DETAILED_STATUS, "Exiting determineCut");
01319     return;
01320   }
01321 
01323   // Divide the coordinates into balanced left and right
01324   // halves.
01325 
01326   ArrayView<lno_t> emptyIndices;
01327   lno_t localCountLeft;
01328 
01329   try{
01330     BSPfindCut<mvector_t>( env, comm,
01331       params, numTestCuts, tolerance,
01332       cutDimension, coordDim, nWeightsPerCoord, vectors, emptyIndices,
01333       fractionLeft.view(0, nVecs),
01334       globalMinCoord, globalMaxCoord,
01335       cutValue, lrflags.view(0, numLocalCoords),
01336       weightLeftHalf, weightRightHalf, localCountLeft, imbalance);
01337   }
01338   Z2_FORWARD_EXCEPTIONS
01339   env->debug(DETAILED_STATUS, "Exiting determineCut");
01340 }
01341 
01342 
01365 template <typename mvector_t, typename Adapter>
01366  void serialRCB(
01367     const RCP<const Environment> &env,
01368     int depth,
01369     const std::bitset<NUM_RCB_PARAMS> &params,
01370     int numTestCuts, 
01371     typename mvector_t::scalar_type tolerance, 
01372     int coordDim,
01373     int nWeightsPerCoord,
01374     const RCP<mvector_t> &vectors, 
01375     ArrayView<typename mvector_t::local_ordinal_type> index,
01376     const RCP<PartitioningSolution<Adapter> > &solution,
01377     partId_t part0, 
01378     partId_t part1,
01379     ArrayView<partId_t> partNum)   // output
01380 {
01381   env->debug(DETAILED_STATUS, "Entering serialRCB");
01382   env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);
01383 
01384   typedef typename mvector_t::scalar_type scalar_t;
01385   typedef typename mvector_t::local_ordinal_type lno_t;
01386 
01387   RCP<Comm<int> > comm(new Teuchos::SerialComm<int>);  
01388 
01389   lno_t numLocalCoords=0;
01390   bool useIndices;
01391   bool firstCall=false;
01392 
01393   if (index.size() == 0){
01394     // First time through there are no indices.
01395     useIndices = false;
01396     numLocalCoords = vectors->getLocalLength();
01397     firstCall = true;
01398   }
01399   else{
01400     useIndices = true;
01401     numLocalCoords = index.size();
01402   }
01403 
01404   if (env->getDebugLevel() >= DETAILED_STATUS){
01405     std::ostringstream info;
01406     info << "  Number of coordinates: " << numLocalCoords << std::endl;
01407     info << "  Use index: " << useIndices << std::endl;
01408     env->debug(DETAILED_STATUS, info.str());
01409   }
01410 
01412   // Are we done?
01413 
01414   if (part1 == part0){
01415     if (useIndices)
01416       for (lno_t i=0; i < numLocalCoords; i++)
01417         partNum[index[i]] = part0;
01418     else
01419       for (lno_t i=0; i < numLocalCoords; i++)
01420         partNum[i] = part0;
01421 
01422     env->memory("serial RCB end");
01423     env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01424     env->debug(DETAILED_STATUS, "Exiting serialRCB");
01425 
01426     return;
01427   }
01428 
01430   // Pick a cut direction
01431 
01432   int cutDimension;
01433   scalar_t minCoord, maxCoord;
01434 
01435   try{
01436     getCutDimension(env, comm, coordDim, vectors, index,
01437       cutDimension, minCoord, maxCoord);
01438   }
01439   Z2_FORWARD_EXCEPTIONS
01440 
01442   // Compute relative sizes of the two halves.
01443 
01444   ArrayRCP<double> fractionLeft;
01445   partId_t numPartsLeftHalf;
01446 
01447   try{
01448     getFractionLeft<Adapter>(env, part0, part1, solution,
01449       fractionLeft, numPartsLeftHalf);
01450   }
01451   Z2_FORWARD_EXCEPTIONS
01452   int nVecs = fractionLeft.size();
01453 
01454   // Check for special case of empty left or empty right.
01455 
01456   scalar_t imbalance, cutValue;  //unused for now
01457   leftRightFlag lrf;
01458 
01459   bool emptyPart = emptyPartsCheck(env, fractionLeft.view(0, nVecs), 
01460     minCoord, maxCoord, lrf, cutValue);
01461 
01462   if (emptyPart){  // continue only on the side that is not empty
01463 
01464     if (lrf == rightFlag)  // right is empty
01465       part1 = part0 + numPartsLeftHalf - 1;
01466     else
01467       part0 = part0 + numPartsLeftHalf;
01468 
01469     if (part0 == part1){
01470       if (useIndices){
01471         for (lno_t i=0; i < numLocalCoords; i++){
01472           partNum[index[i]] = part0;
01473         }
01474       }
01475       else{
01476         for (lno_t i=0; i < numLocalCoords; i++){
01477           partNum[i] = part0;
01478         }
01479       }
01480       
01481       imbalance = 0.0;       // perfect
01482       env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01483       env->debug(DETAILED_STATUS, "Exiting serialRCB");
01484       return;
01485     }
01486   }
01487 
01489   // Divide into balanced left and right halves.
01490 
01491   scalar_t totalLeft=0, totalRight=0;
01492   lno_t localCountLeft=0, localCountRight=0;
01493   unsigned char *flags = new unsigned char [numLocalCoords];
01494   env->localMemoryAssertion(__FILE__, __LINE__, numLocalCoords, flags) ;
01495   ArrayRCP<unsigned char> lrflags(flags, 0, numLocalCoords, true);
01496 
01497   try{
01498     BSPfindCut<mvector_t>( env, comm,
01499       params, numTestCuts, tolerance,
01500       cutDimension, coordDim, nWeightsPerCoord, vectors, index,
01501       fractionLeft.view(0, nVecs),
01502       minCoord, maxCoord,
01503       cutValue, lrflags.view(0, numLocalCoords),
01504       totalLeft, totalRight, localCountLeft, imbalance);
01505   }
01506   Z2_FORWARD_EXCEPTIONS
01507 
01508   if (firstCall)
01509     env->memory("serial RCB start");
01510 
01512   // Adjust indices for left half and right half
01513 
01514  localCountRight = numLocalCoords - localCountLeft;
01515 
01516   if (localCountLeft){
01517 
01518     lno_t *newIndex = new lno_t [localCountLeft];
01519     env->localMemoryAssertion(__FILE__, __LINE__, localCountLeft, newIndex);
01520     ArrayView<lno_t> leftIndices(newIndex, localCountLeft);
01521 
01522     if (useIndices){
01523       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01524         if (lrflags[i] == leftFlag)
01525           newIndex[ii++] = index[i];
01526     }
01527     else{
01528       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01529         if (lrflags[i] == leftFlag)
01530           newIndex[ii++] = i;
01531     }
01532 
01533     partId_t newPart1 = part0 + numPartsLeftHalf - 1;
01534 
01535     env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01536 
01537     serialRCB<mvector_t, Adapter>(env, depth+1, params, numTestCuts, tolerance, 
01538       coordDim, nWeightsPerCoord, vectors, leftIndices,
01539       solution, part0, newPart1, partNum);
01540 
01541     env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);
01542 
01543     delete [] newIndex;
01544   }
01545 
01546 
01547   if (localCountRight){
01548 
01549     lno_t *newIndex = new lno_t [localCountRight];
01550     env->localMemoryAssertion(__FILE__, __LINE__, localCountRight, newIndex);
01551     ArrayView<lno_t> rightIndices(newIndex, localCountRight);
01552 
01553     if (useIndices){
01554       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01555         if (lrflags[i] == rightFlag){
01556           newIndex[ii++] = index[i];
01557         }
01558     }
01559     else{
01560       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01561         if (lrflags[i] == rightFlag)
01562           newIndex[ii++] = i;
01563     }
01564 
01565     partId_t newPart0 = part0 + numPartsLeftHalf;
01566 
01567     env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01568 
01569     serialRCB<mvector_t, Adapter>(env, depth+1, params, numTestCuts, tolerance, 
01570       coordDim, nWeightsPerCoord, vectors, rightIndices,
01571       solution, newPart0, part1, partNum);
01572 
01573     env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);
01574 
01575     delete [] newIndex;
01576   }
01577 
01578   env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01579   env->debug(DETAILED_STATUS, "Exiting serialRCB");
01580 }
01581 
01582 }// namespace Zoltan2
01583 
01584 #endif
01585