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