Zoltan2 Version of the Day
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, string("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 weightDim = solution->getNumberOfCriteria();
00198   fractionLeft = arcp(new double [weightDim], 0, weightDim);
00199 
00200   for (int wdim=0; wdim<weightDim; wdim++){
00201     if (solution->criteriaHasUniformPartSizes(wdim)){
00202       fractionLeft[wdim] = double(numPartsLeftHalf) / double(numParts);
00203     }
00204     else{
00205       fractionLeft[wdim] = 0;
00206       for(int partId=left0; partId <= left1; partId++){
00207         fractionLeft[wdim] += solution->getCriteriaPartSize(wdim, partId);
00208       }
00209       double total = fractionLeft[wdim];
00210       for(int partId=right0; partId <= right1; partId++){
00211         total += solution->getCriteriaPartSize(wdim, partId);
00212       }
00213       fractionLeft[wdim] /= total;
00214     }
00215   }
00216   env->debug(DETAILED_STATUS, string("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, string("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++] = numeric_limits<scalar_t>::max();
00286         spans[next++] = 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       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, string("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, string("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, string("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 weightDim = weights.size();
00479     
00480   if (weightDim > 1){
00481     Array<scalar_t> coordWeights(weightDim, 1.0);
00482     for (size_t wdim=0; wdim < weightDim; wdim++){
00483       if (weights[wdim].size() > 0)
00484         coordWeights[wdim] = weights[wdim][id];
00485     }
00486 
00487     wgtValue = normedWeight<scalar_t>(coordWeights.view(0,weightDim), mcnorm);
00488   }
00489   else if (weights[0].size() > 0){
00490     wgtValue = weights[0][id];
00491   }
00492 
00493   return wgtValue;
00494 }
00495 
00513 template <typename scalar_t>
00514   bool emptyPartsCheck( const RCP<const Environment> &env,
00515    const ArrayView<double> fractionLeft, 
00516    scalar_t minCoord, scalar_t maxCoord,
00517    leftRightFlag &lrf,
00518    scalar_t &cutValue)
00519 {
00520   env->debug(DETAILED_STATUS, string("Entering emptyPartsCheck"));
00521   // initialize return values
00522   lrf = leftFlag;
00523   cutValue = 0.0;
00524 
00525   size_t weightDim = fractionLeft.size();
00526   int numEmptyRight = 0, numEmptyLeft = 0;
00527 
00528   for (size_t i=0; i < weightDim; i++){
00529     if (fractionLeft[i] == 0.0)
00530       numEmptyLeft++;
00531     else if (fractionLeft[i] == 1.0)
00532       numEmptyRight++;
00533   }
00534 
00535   if (numEmptyRight + numEmptyLeft == 0)
00536     return false;
00537 
00538   env->localInputAssertion(__FILE__, __LINE__,
00539     "partitioning not solvable - conflicting part size criteria",
00540     (numEmptyRight==0) || (numEmptyLeft==0), BASIC_ASSERTION);
00541 
00542   if (numEmptyLeft > 0){
00543     lrf = leftFlag;
00544     cutValue = minCoord;
00545   }
00546   else{
00547     lrf = rightFlag;
00548     cutValue = maxCoord;
00549   }
00550 
00551   env->debug(DETAILED_STATUS, string("Exiting emptyPartsCheck"));
00552   return true;
00553 }
00554 
00567 template <typename lno_t, typename gno_t, typename scalar_t>
00568   void testCoordinatesOnRightBoundary(
00569     const RCP<const Environment> &env,
00570     const RCP<Comm<int> > &comm, 
00571     scalar_t totalWeightLeft,
00572     scalar_t targetLeftScalar,
00573     int cutLocation,
00574     ArrayView<scalar_t> localSums,
00575     ArrayView<scalar_t> globalSums,
00576     ArrayView<lno_t> index,
00577     ArrayView<StridedData<lno_t, scalar_t> > weights,
00578     multiCriteriaNorm mcnorm,
00579     ArrayView<unsigned char> lrFlags,
00580     scalar_t &globalWeightMovedRight)   // output
00581 {
00582   env->debug(DETAILED_STATUS, string("Entering testCoordinatesOnRightBoundary"));
00583   int nprocs = comm->getSize();
00584   int rank = comm->getRank();
00585 
00586   globalWeightMovedRight = 0.0;
00587 
00588   scalar_t localBoundarySum = localSums[cutLocation];
00589 
00590   scalar_t total = totalWeightLeft;
00591   for (int i=0; i <= cutLocation; i++)
00592     total += globalSums[i];
00593 
00594   scalar_t totalMoveRight = total - targetLeftScalar;
00595   scalar_t localMoveRight = localBoundarySum;
00596   scalar_t actualWeightMovedRight = 0.0;
00597   
00598   Array<scalar_t> scansum(nprocs+1, 0.0);
00599   Teuchos::gatherAll<int, scalar_t>(*comm, 1, 
00600     &localBoundarySum, nprocs, scansum.getRawPtr()+1);
00601   for (int i=2; i<=nprocs; i++)
00602     scansum[i] += scansum[i-1];
00603 
00604   if (localBoundarySum > 0.0){
00605 
00606     scalar_t sumMine = scansum[rank]; // sum of ranks preceding me
00607     scalar_t diff = scansum[nprocs] - sumMine;
00608     localMoveRight = 0;
00609 
00610     if (diff <= totalMoveRight)
00611       localMoveRight = localBoundarySum;  // all
00612     else{
00613       scalar_t leftPart = diff - totalMoveRight;
00614       if (leftPart < localBoundarySum)
00615         localMoveRight = localBoundarySum - leftPart;
00616     }
00617   }
00618 
00619   if (localMoveRight > 0.0){
00620 
00621     bool moveAll =  (localMoveRight >= localBoundarySum);
00622 
00623     if (moveAll){
00624       actualWeightMovedRight = localBoundarySum;
00625       for (lno_t i=0; i < lrFlags.size(); i++){
00626         if (lrFlags[i] == cutLocation){
00627           lrFlags[i] = rightFlag;
00628         }
00629       }
00630     }
00631     else{
00632       int weightDim = weights.size();
00633       int useIndices = index.size() > 0;
00634 
00635       for (lno_t i=0; i < lrFlags.size(); i++){
00636         if (lrFlags[i] == cutLocation){
00637           lrFlags[i] = rightFlag;
00638 
00639           lno_t idx = (useIndices ? index[i] : i);
00640 
00641           actualWeightMovedRight += getCoordWeight<lno_t, scalar_t>(
00642             idx, mcnorm, weights.view(0, weightDim));
00643 
00644           if (actualWeightMovedRight >= localMoveRight)
00645             break;
00646         }
00647       } // next coordinate
00648     }
00649   }
00650 
00651   try{
00652     reduceAll<int, scalar_t>(
00653       *comm, Teuchos::REDUCE_SUM, 1, &actualWeightMovedRight,
00654       &globalWeightMovedRight);
00655   }
00656   Z2_THROW_OUTSIDE_ERROR(*env)
00657 
00658   env->debug(DETAILED_STATUS, string("Exiting testCoordinatesOnRightBoundary"));
00659   return;
00660 }
00661 
00709 template <typename mvector_t>
00710   void BSPfindCut(
00711     const RCP<const Environment> &env,
00712     const RCP<Comm<int> > &comm,
00713     const std::bitset<NUM_RCB_PARAMS> &params,
00714     int numTestCuts,
00715     typename mvector_t::scalar_type tolerance,
00716     int cutDim,
00717     int coordDim,
00718     const RCP<mvector_t> &vectors,
00719     ArrayView<typename mvector_t::local_ordinal_type> index,
00720     ArrayView<double> fractionLeft,
00721     ArrayView<bool> uniformWeights,
00722     typename mvector_t::scalar_type coordGlobalMin,
00723     typename mvector_t::scalar_type coordGlobalMax,
00724     typename mvector_t::scalar_type &cutValue,         // output
00725     ArrayView<unsigned char> lrFlags,                  // output
00726     typename mvector_t::scalar_type &totalWeightLeft,  // output
00727     typename mvector_t::scalar_type &totalWeightRight, // output
00728     typename mvector_t::local_ordinal_type &localCountLeft, // output
00729     typename mvector_t::scalar_type &imbalance)        // output
00730 {
00731   env->debug(DETAILED_STATUS, string("Entering BSPfindCut"));
00732 
00733   // initialize output
00734   bool useIndices = index.size() > 0;
00735   int numAllCoords = vectors->getLocalLength();
00736   int numCoords = (useIndices ? index.size() : numAllCoords);
00737   cutValue = totalWeightLeft = totalWeightRight = imbalance = 0.0;
00738   localCountLeft = 0;
00739 
00740   for (int i=0; i < numCoords; i++)
00741     lrFlags[i] = unsetFlag;
00742 
00743   typedef typename mvector_t::scalar_type scalar_t;
00744   typedef typename mvector_t::local_ordinal_type lno_t;
00745   typedef typename mvector_t::global_ordinal_type gno_t;
00746   typedef StridedData<lno_t, scalar_t> input_t;
00747 
00748   bool multiplePartSizeSpecs = params.test(rcb_multiplePartSizeSpecs);
00749   bool rectilinearBlocks = params.test(rcb_rectilinearBlocks);
00750   bool averageCuts = params.test(rcb_averageCuts);
00751 
00752   // A coordinate is considered to be on a cut if it is within
00753   // this distance of the cut.
00754 
00755   double epsilon = (coordGlobalMax - coordGlobalMin) * 10e-9;
00756 
00757   // Find the coordinate values and weights.
00758 
00759   int weightDim = uniformWeights.size();   // at least one
00760   int numNonUniformWeights = 0;
00761   for (int i=0; i < weightDim; i++){
00762     if (!uniformWeights[i])
00763       numNonUniformWeights++;
00764   }
00765 
00766   if (env->getDebugLevel() >= DETAILED_STATUS){
00767     ostringstream info;
00768     info << "Weight dim " << weightDim << ", Fraction left:";
00769     for (int i=0; i < weightDim; i++)
00770       info << " " << fractionLeft[i];
00771     info << endl << "Dimension " << cutDim << " [";
00772     info << coordGlobalMin << ", " << coordGlobalMax << "]";
00773     info << endl << "# test cuts " << numTestCuts;
00774     info << ", tolerance " << tolerance << endl;
00775     env->debug(DETAILED_STATUS, info.str());
00776   }
00777 
00778   const scalar_t *coordValue = vectors->getData(cutDim).getRawPtr();
00779 
00780   // An empty input_t object implies uniform weights.
00781 
00782   input_t *wgtinfo = new input_t [weightDim];
00783   env->localMemoryAssertion(__FILE__, __LINE__, weightDim, wgtinfo);
00784   ArrayRCP<input_t> weight(wgtinfo, 0, weightDim, true);
00785 
00786   if (numNonUniformWeights > 0){
00787     for (int wdim = 0, widx=coordDim; wdim < weightDim; wdim++){
00788       if (!uniformWeights[wdim]){
00789         weight[wdim] = input_t(vectors->getData(widx++), 1);
00790       }
00791     }
00792   }
00793 
00794   // Multicriteria norm
00795 
00796   multiCriteriaNorm mcnorm = normBalanceTotalMaximum;
00797 
00798   if (params.test(rcb_minMaximumWeight))
00799     mcnorm = normMinimizeMaximumWeight;
00800   else if (params.test(rcb_minTotalWeight))
00801     mcnorm = normMinimizeTotalWeight;
00802   
00803   // Goal is globally find one cut that comes close to leaving
00804   // partSizeLeft*totalWeight on the left side.
00805 
00806   Epetra_SerialDenseVector partSizeLeft( 
00807     View, fractionLeft.getRawPtr(), weightDim);
00808   
00809   // Where do we make the first test cuts?
00810   //
00811   //   min     1     2     3     max
00812   //
00813   // Space is [min, max].  
00814   // If there are three test cuts: 1, 2 and 3.
00815   //   4 regions: [min,1] [1,2] [2,3] [3,max].
00816   //   5 boundaries: min, 1, 2, 3, and max.
00817 
00818   int numRegions = numTestCuts + 1;
00819   int numBoundaries = numTestCuts + 2;
00820   int endBoundary = numBoundaries - 1;
00821   vector<scalar_t> boundaries(numBoundaries);
00822   vector<scalar_t> searchBoundaries(numBoundaries);
00823 
00824   int numSums = numBoundaries+numRegions;
00825 
00826   Teuchos::Zoltan2_RCBOperation<int, scalar_t> reductionOp(
00827     numSums,      // number of sums
00828     numRegions,   // number of mins
00829     numRegions);  // number of maxes
00830 
00831   typename std::vector<scalar_t>::iterator foundCut;
00832 
00833   bool done=false;
00834   bool fail=false;
00835   scalar_t min = coordGlobalMin;
00836   scalar_t max = coordGlobalMax;
00837   lno_t numRemaining = numCoords;
00838   lno_t prevNumRemaining = numCoords;
00839   size_t numGlobalPoints = vectors->getGlobalLength();
00840   size_t sanityCheck = numGlobalPoints;
00841 
00842   double totalWeight = 0;
00843   double targetLeftScalar = 0;
00844   double targetLeftNorm = 0;
00845   Epetra_SerialDenseVector targetLeftVector(weightDim);
00846 
00847   while (!done && !fail && sanityCheck--){
00848 
00849     // Create regions into which coordinates will be placed.
00850     scalar_t diff = (max - min) / numRegions;
00851     boundaries[0] = min;
00852     boundaries[endBoundary] = max;
00853     for (int i=0; i < endBoundary; i++){
00854       searchBoundaries[i+1] = boundaries[i+1] = boundaries[i] + diff;
00855     }
00856 
00857     // Move ends slightly so we catch points on boundary.
00858     searchBoundaries[0] = min - epsilon;
00859     searchBoundaries[endBoundary] = max + epsilon;
00860 
00861     // Save region and boundary sums, and region min and max.
00862  
00863     ArrayRCP<scalar_t> localSums = reductionOp.getInitializedBuffer(
00864       searchBoundaries[endBoundary], searchBoundaries[0]);
00865  
00866     ArrayRCP<scalar_t> globalSums = reductionOp.getInitializedBuffer(
00867       searchBoundaries[endBoundary], searchBoundaries[0]);
00868 
00869     scalar_t *sums = localSums.getRawPtr();
00870     scalar_t *regionMin = sums + numSums;
00871     scalar_t *regionMax = regionMin + numRegions;
00872 
00873     if (numRemaining > 0){
00874 
00875       // Assign each of my points to a region.
00876       // lower_bound() finds the first cut f, such that f >= coordValue[i].
00877       // So for now, objects that are on the cut boundary go into the
00878       // region on the "left" side.
00879 
00880       for (lno_t i=0; i < numCoords; i++){
00881   
00882         if (lrFlags[i] != unsetFlag)
00883           continue;
00884 
00885         int inRegion = 0;
00886         int idx = (useIndices ? index[i] : i);
00887         scalar_t value = coordValue[idx];
00888 
00889         if (numRegions > 2){
00890        
00891           foundCut = std::lower_bound(
00892             searchBoundaries.begin(), searchBoundaries.end(), value);
00893           
00894           env->localBugAssertion(__FILE__, __LINE__, "search cuts", 
00895             foundCut != searchBoundaries.end(), BASIC_ASSERTION);
00896 
00897           inRegion = foundCut - searchBoundaries.begin() - 1;        
00898         }
00899         else{
00900           if (value <= boundaries[1])
00901             inRegion = 0;
00902           else
00903             inRegion = 1;
00904         }
00905 
00906         int sumIdx = 1 + inRegion*2;
00907 
00908         if (value >= boundaries[inRegion+1]-epsilon){  
00909           // "on" right boundary of this region
00910           sumIdx++;
00911         }
00912         else if (inRegion==0 && (value < min + epsilon)){
00913           // in region 0 but far left boundary
00914           sumIdx--;
00915         }
00916 
00917         lrFlags[i] = (unsigned char)sumIdx;
00918 
00919         if (value < regionMin[inRegion])
00920           regionMin[inRegion] = value;
00921         if (value > regionMax[inRegion])
00922           regionMax[inRegion] = value;
00923 
00924         sums[sumIdx] += getCoordWeight<lno_t, scalar_t>(idx, 
00925               mcnorm, weight.view(0, weightDim));
00926 
00927       } // next coord
00928     }
00929 
00930     try{
00931       reduceAll<int, scalar_t>(*comm, reductionOp,
00932         localSums.size(), localSums.getRawPtr(), globalSums.getRawPtr());
00933     }
00934     Z2_THROW_OUTSIDE_ERROR(*env)
00935 
00936     sums = globalSums.getRawPtr();
00937     regionMin = sums + numSums;
00938     regionMax = regionMin + numRegions;
00939 
00940     if (env->getDebugLevel() >= DETAILED_STATUS){
00941       ostringstream info;
00942       info << "  Region " << min << " - " << max << endl;
00943       info << "  Remaining to classify: " << numRemaining << endl;
00944       info << "  Boundaries: ";
00945       for (int i=0; i < numBoundaries; i++)
00946         info << boundaries[i] << " ";
00947       info << endl << "  For searching: ";
00948       for (int i=0; i < numBoundaries; i++)
00949         info << searchBoundaries[i] << " ";
00950       info << endl << "  Global sums: ";
00951       double sumTotal=0;
00952       for (int i=0; i < numSums; i++){
00953         sumTotal += sums[i];
00954         info << sums[i] << " ";
00955       }
00956       info << " total: " << sumTotal << endl;
00957       env->debug(DETAILED_STATUS, info.str());
00958     }
00959 
00960     if (totalWeight == 0){   // first time through only
00961 
00962       for (int i=0; i < numSums; i++)
00963         totalWeight += sums[i];
00964 
00965       partSizeLeft.Scale(totalWeight);
00966       targetLeftVector = partSizeLeft;
00967 
00968       targetLeftScalar = targetLeftVector[0];
00969       targetLeftNorm = targetLeftVector.Norm2();
00970       totalWeightLeft = 0;
00971       totalWeightRight = 0;
00972     }
00973 
00974     int cutLocation=0;
00975     scalar_t testDiff=0, prevTestDiff=0, target=0;
00976 
00977     if (multiplePartSizeSpecs){
00978       // more complex: if we have multiple weight dimensions, the
00979       //   weights are non-uniform, and the part sizes requested
00980       //   for each each weight dimension differ, then we may not
00981       //   be able to reach the imbalance tolerance.
00982       //
00983       // TODO: discuss how to (or whether to) handle this case.
00984       // 
00985       // For now we look at this imbalance:
00986       // 
00987       //    |target - actual|^2 / |target|^2
00988   
00989       target = targetLeftNorm;
00990       Epetra_SerialDenseVector testVec(weightDim);
00991       for (int i=0; i < weightDim; i++)
00992         testVec[i] = totalWeightLeft;
00993       Epetra_SerialDenseVector diffVec = testVec;
00994       diffVec.Scale(-1.0);
00995       diffVec += targetLeftVector;
00996 
00997       testDiff = diffVec.Norm2(); // imbalance numerator
00998       prevTestDiff = testDiff;
00999       cutLocation= -1;
01000 
01001       while (++cutLocation< numSums){
01002 
01003         for (int i=0; i < weightDim; i++)
01004           testVec[i] += sums[cutLocation];
01005   
01006         diffVec = testVec;
01007         diffVec.Scale(-1.0);
01008         diffVec += targetLeftVector;
01009   
01010         testDiff = diffVec.Norm2();
01011         
01012         if (testDiff >= target)
01013           break;
01014 
01015         prevTestDiff = testDiff;
01016       }
01017     }
01018     else{    // the part sizes for each weight dimension are the same
01019   
01020       target = targetLeftScalar;
01021       testDiff = totalWeightLeft; 
01022       prevTestDiff = testDiff;
01023       cutLocation = -1;
01024   
01025       while (++cutLocation < numSums){
01026      
01027         testDiff += sums[cutLocation];
01028         if (testDiff >= target)
01029           break;
01030   
01031         prevTestDiff = testDiff;
01032       }
01033     }
01034 
01035     scalar_t diffLeftCut = target - prevTestDiff;
01036     scalar_t diffRightCut = testDiff - target;
01037 
01038     if (diffLeftCut < diffRightCut){
01039       imbalance = diffLeftCut / target;
01040       if (imbalance <= tolerance){
01041         env->debug(DETAILED_STATUS, "  Done, tolerance met");
01042         done = true;
01043         cutLocation--;
01044       }
01045     } 
01046     else{
01047       imbalance = diffRightCut / target;
01048       if (imbalance <= tolerance){
01049         env->debug(DETAILED_STATUS, "  Done, tolerance met");
01050         done = true;
01051      }
01052     }
01053 
01054     bool cutLocIsRegion = (cutLocation % 2 == 1);
01055     bool cutLocIsBoundary = !cutLocIsRegion;
01056 
01057     if (env->getDebugLevel() >= DETAILED_STATUS){
01058       ostringstream info;
01059       info << "  Best cut location: " << cutLocation;
01060       if (cutLocIsRegion) info << " just after a region." << endl;
01061       else info << " just after a boundary." << endl;
01062       env->debug(DETAILED_STATUS, info.str());
01063     }
01064 
01065     if (!done && cutLocIsBoundary){
01066 
01067       done = true;    // can not divide space any more
01068 
01069       env->debug(DETAILED_STATUS, "  Done, cutting at a region boundary");
01070 
01071       if (rectilinearBlocks){
01072         // Can not divide boundary points into two
01073         // different regions to achieve balance.
01074         fail = true;
01075       }
01076       else {
01077         // Maybe by moving some of the points right, we can
01078         // obtain a better balance.  If so, the lrFlag for
01079         // the points moved right will be updated here.
01080 
01081         scalar_t globalWeightMovedRight(0);
01082 
01083         testCoordinatesOnRightBoundary<lno_t, gno_t, scalar_t>(
01084           env, comm, totalWeightLeft, targetLeftScalar, cutLocation, 
01085           localSums.view(0, numSums), globalSums.view(0, numSums),
01086           index, weight.view(0, weightDim), mcnorm, lrFlags,
01087           globalWeightMovedRight);
01088 
01089         scalar_t newSum = testDiff - globalWeightMovedRight;
01090         globalSums[cutLocation] -= globalWeightMovedRight;
01091 
01092         if (newSum > target)
01093           imbalance = (target - newSum) / target;
01094         else
01095           imbalance = (newSum - target) / target;
01096 
01097         if (imbalance > tolerance)
01098           fail = true;
01099       }
01100     }
01101 
01102     int rightmostLeftNum=0, leftmostRightNum=0;
01103 
01104     if (!done){
01105       // Best cut is following a region.  Narrow down the boundaries.
01106 
01107       int regionNumber = (cutLocation - 1) / 2;
01108 
01109       min = regionMin[regionNumber];
01110       max = regionMax[regionNumber];
01111 
01112       rightmostLeftNum = cutLocation - 1;
01113       leftmostRightNum = cutLocation + 1;
01114     }
01115     else{
01116       rightmostLeftNum = cutLocation;
01117       leftmostRightNum = cutLocation + 1;
01118 
01119       if (cutLocIsRegion && averageCuts){
01120         int regionNumber = (cutLocation - 1) / 2;
01121         cutValue = (
01122           boundaries[regionNumber+1] + // boundary to right
01123           regionMax[regionNumber])     // rightmost point in region
01124           / 2.0;
01125       }
01126     }
01127 
01128     for (int i=0; i <= rightmostLeftNum; i++){
01129       totalWeightLeft += globalSums[i];
01130     }
01131 
01132     for (lno_t i=0; i < numCoords; i++){
01133       if (lrFlags[i] != leftFlag && lrFlags[i] != rightFlag){
01134         if (lrFlags[i] <= rightmostLeftNum){
01135           lrFlags[i] = leftFlag;
01136           numRemaining--;
01137           localCountLeft++;
01138         }
01139         else if (lrFlags[i] >= leftmostRightNum){
01140           lrFlags[i] = rightFlag;
01141           numRemaining--;
01142         }
01143         else{
01144           lrFlags[i] = unsetFlag;   // still to be determined
01145         }
01146       }
01147     }
01148 
01149     if (env->getDebugLevel() >= VERBOSE_DETAILED_STATUS && numCoords < 100){
01150       // For large numCoords, building this message
01151       // takes an extraordinarily long time.
01152       ostringstream ossLeft;
01153       ostringstream ossRight;
01154       ossLeft << "left: ";
01155       ossRight << "right: ";
01156       for (lno_t i=0; i < numCoords; i++){
01157         if (lrFlags[i] == unsetFlag)
01158           continue;
01159         lno_t idx = (useIndices ? index[i] : i);
01160         scalar_t val = coordValue[idx];
01161         if (lrFlags[i] == leftFlag)
01162           ossLeft << val << " ";
01163         else if (lrFlags[i] == rightFlag)
01164           ossRight << val << " ";
01165         else 
01166           env->localBugAssertion(__FILE__, __LINE__, 
01167             "left/right flags", false, BASIC_ASSERTION);
01168       }
01169       ostringstream msg;
01170       msg << ossLeft.str() << endl << ossRight.str() << endl;
01171       env->debug(VERBOSE_DETAILED_STATUS, msg.str());
01172     }
01173 
01174     prevNumRemaining = numRemaining;
01175   }  // while !done
01176 
01177   totalWeightRight = totalWeight - totalWeightLeft;
01178 
01179   if (fail)
01180     env->debug(BASIC_STATUS, "Warning: tolerance not achieved in sub-step");
01181 
01182   // TODO: If fail the warn that tolerance was not met.
01183 
01184   env->globalInputAssertion(__FILE__, __LINE__, "partitioning not solvable",
01185     done, DEBUG_MODE_ASSERTION, comm);
01186 
01187   env->memory("End of bisection");
01188 
01189   if (env->getDebugLevel() >= DETAILED_STATUS){
01190     ostringstream oss;
01191     oss << "Exiting BSPfindCut, ";
01192     oss << "# iterations: " << numGlobalPoints - sanityCheck;    
01193     env->debug(DETAILED_STATUS, oss.str());
01194   }
01195 
01196   return;
01197 }
01198 
01226 template <typename mvector_t, typename Adapter>
01227   void determineCut(
01228     const RCP<const Environment> &env,
01229     const RCP<Comm<int> > &comm,
01230     const std::bitset<NUM_RCB_PARAMS> &params,
01231     int numTestCuts,
01232     typename mvector_t::scalar_type tolerance,
01233     int coordDim, 
01234     const RCP<mvector_t> &vectors,
01235     const ArrayView<bool> uniformWeights,
01236     multiCriteriaNorm mcnorm,
01237     const RCP<PartitioningSolution<Adapter> > &solution,
01238     partId_t part0, 
01239     partId_t part1,
01240     ArrayView<unsigned char> lrflags,  // output
01241     int &cutDimension,                  // output
01242     typename mvector_t::scalar_type &cutValue,   // output
01243     typename mvector_t::scalar_type &imbalance,  // output
01244     partId_t &numPartsLeftHalf,                  // output
01245     typename mvector_t::scalar_type &weightLeftHalf,  // output
01246     typename mvector_t::scalar_type &weightRightHalf  // output
01247     )
01248 {
01249   env->debug(DETAILED_STATUS, string("Entering determineCut"));
01250   typedef typename mvector_t::scalar_type scalar_t;
01251   typedef typename mvector_t::local_ordinal_type lno_t;
01252   typedef typename mvector_t::global_ordinal_type gno_t;
01253   typedef StridedData<lno_t, scalar_t> input_t;
01254 
01255   lno_t numLocalCoords = vectors->getLocalLength();
01256 
01257   // initialize return values
01258   cutDimension = 0;
01259   cutValue = imbalance = weightLeftHalf = weightRightHalf = 0.0;
01260   numPartsLeftHalf = 0;
01261 
01263   // Pick a cut direction.
01264 
01265   scalar_t globalMinCoord, globalMaxCoord;
01266   ArrayView<lno_t> emptyIndex;
01267 
01268   getCutDimension<mvector_t>(env, comm, coordDim, vectors, emptyIndex,
01269     cutDimension, globalMinCoord, globalMaxCoord);
01270 
01272   // Compute part sizes for the two parts.
01273 
01274   ArrayRCP<double> fractionLeft;
01275   size_t weightDim = uniformWeights.size();
01276 
01277   getFractionLeft<Adapter>(env, part0, part1, solution,
01278     fractionLeft, numPartsLeftHalf);
01279 
01280   // Special case of empty left or empty right.
01281 
01282   leftRightFlag lrf;
01283 
01284   bool emptyParts = emptyPartsCheck(env,
01285     fractionLeft.view(0, weightDim), // input
01286     globalMinCoord, globalMaxCoord,  // input
01287     lrf, cutValue);                  // output
01288 
01289   if (emptyParts){
01290     
01291     for (lno_t i=0; i < lrflags.size(); i++)
01292       lrflags[i] = lrf;
01293 
01294     imbalance = 0.0;                // perfect balance
01295     scalar_t totalWeight = 0.0;
01296     int numNonUniform = 0;
01297 
01298     for (size_t i=0; i < weightDim; i++)
01299       if (!uniformWeights[i])
01300         numNonUniform++;
01301 
01302     int wgt1 = vectors->getNumVectors() - numNonUniform;
01303 
01304     if (weightDim == 1){
01305       if (numNonUniform == 0)
01306         totalWeight = numLocalCoords;
01307       else{
01308         const scalar_t *val = vectors->getData(wgt1).getRawPtr();
01309         for (lno_t i=0; i < numLocalCoords; i++)
01310           totalWeight += val[i];
01311       }
01312     }
01313     else{  // need to add up total normed weight
01314       Array<input_t> wgts(weightDim);
01315       for (size_t i=0; i < weightDim; i++){
01316         if (!uniformWeights[i]){
01317           wgts[i] = input_t(vectors->getData(wgt1++), 1);
01318         }
01319       }
01320 
01321       partId_t numParts, numNonemptyParts;
01322       ArrayRCP<MetricValues<scalar_t> > metrics;
01323       ArrayRCP<scalar_t> weightSums;
01324     
01325       globalSumsByPart<scalar_t, unsigned char, lno_t>(
01326         env, comm, lrflags, 
01327         wgts.view(0, weightDim), mcnorm,
01328         numParts, numNonemptyParts, metrics, weightSums);
01329 
01330       totalWeight = weightSums[1];  // sum normed weight
01331     }
01332 
01333     if (lrf == leftFlag){
01334       numPartsLeftHalf = part1 - part0 + 1;
01335       weightLeftHalf = totalWeight;
01336       weightRightHalf = 0;
01337     }
01338     else{
01339       numPartsLeftHalf = 0;
01340       weightRightHalf = totalWeight;
01341       weightLeftHalf = 0;
01342     }
01343 
01344     env->debug(DETAILED_STATUS, string("Exiting determineCut"));
01345     return;
01346   }
01347 
01349   // Divide the coordinates into balanced left and right
01350   // halves.
01351 
01352   ArrayView<lno_t> emptyIndices;
01353   lno_t localCountLeft;
01354 
01355   try{
01356     BSPfindCut<mvector_t>( env, comm,
01357       params, numTestCuts, tolerance,
01358       cutDimension, coordDim, vectors, emptyIndices,
01359       fractionLeft.view(0, weightDim), uniformWeights.view(0, weightDim),
01360       globalMinCoord, globalMaxCoord,
01361       cutValue, lrflags.view(0, numLocalCoords),
01362       weightLeftHalf, weightRightHalf, localCountLeft, imbalance);
01363   }
01364   Z2_FORWARD_EXCEPTIONS
01365   env->debug(DETAILED_STATUS, string("Exiting determineCut"));
01366 }
01367 
01368 
01392 template <typename mvector_t, typename Adapter>
01393  void serialRCB(
01394     const RCP<const Environment> &env,
01395     int depth,
01396     const std::bitset<NUM_RCB_PARAMS> &params,
01397     int numTestCuts, 
01398     typename mvector_t::scalar_type tolerance, 
01399     int coordDim,
01400     const RCP<mvector_t> &vectors, 
01401     ArrayView<typename mvector_t::local_ordinal_type> index,
01402     const ArrayView<bool> uniformWeights,
01403     const RCP<PartitioningSolution<Adapter> > &solution,
01404     partId_t part0, 
01405     partId_t part1,
01406     ArrayView<partId_t> partNum)   // output
01407 {
01408   env->debug(DETAILED_STATUS, string("Entering serialRCB"));
01409   env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);
01410 
01411   typedef typename mvector_t::scalar_type scalar_t;
01412   typedef typename mvector_t::local_ordinal_type lno_t;
01413 
01414   RCP<Comm<int> > comm(new Teuchos::SerialComm<int>);  
01415 
01416   lno_t numLocalCoords=0;
01417   bool useIndices;
01418   bool firstCall=false;
01419 
01420   if (index.size() == 0){
01421     // First time through there are no indices.
01422     useIndices = false;
01423     numLocalCoords = vectors->getLocalLength();
01424     firstCall = true;
01425   }
01426   else{
01427     useIndices = true;
01428     numLocalCoords = index.size();
01429   }
01430 
01431   if (env->getDebugLevel() >= DETAILED_STATUS){
01432     ostringstream info;
01433     info << "  Number of coordinates: " << numLocalCoords << endl;
01434     info << "  Use index: " << useIndices << endl;
01435     env->debug(DETAILED_STATUS, info.str());
01436   }
01437 
01439   // Are we done?
01440 
01441   if (part1 == part0){
01442     if (useIndices)
01443       for (lno_t i=0; i < numLocalCoords; i++)
01444         partNum[index[i]] = part0;
01445     else
01446       for (lno_t i=0; i < numLocalCoords; i++)
01447         partNum[i] = part0;
01448 
01449     env->memory("serial RCB end");
01450     env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01451     env->debug(DETAILED_STATUS, string("Exiting serialRCB"));
01452 
01453     return;
01454   }
01455 
01457   // Pick a cut direction
01458 
01459   int cutDimension;
01460   scalar_t minCoord, maxCoord;
01461 
01462   try{
01463     getCutDimension(env, comm, coordDim, vectors, index,
01464       cutDimension, minCoord, maxCoord);
01465   }
01466   Z2_FORWARD_EXCEPTIONS
01467 
01469   // Compute relative sizes of the two halves.
01470 
01471   ArrayRCP<double> fractionLeft;
01472   partId_t numPartsLeftHalf;
01473 
01474   try{
01475     getFractionLeft<Adapter>(env, part0, part1, solution,
01476       fractionLeft, numPartsLeftHalf);
01477   }
01478   Z2_FORWARD_EXCEPTIONS
01479 
01480   // Check for special case of empty left or empty right.
01481 
01482   int weightDim = uniformWeights.size();
01483   scalar_t imbalance, cutValue;  //unused for now
01484   leftRightFlag lrf;
01485 
01486   bool emptyPart = emptyPartsCheck(env, fractionLeft.view(0, weightDim), 
01487     minCoord, maxCoord, lrf, cutValue);
01488 
01489   if (emptyPart){  // continue only on the side that is not empty
01490 
01491     if (lrf == rightFlag)  // right is empty
01492       part1 = part0 + numPartsLeftHalf - 1;
01493     else
01494       part0 = part0 + numPartsLeftHalf;
01495 
01496     if (part0 == part1){
01497       if (useIndices){
01498         for (lno_t i=0; i < numLocalCoords; i++){
01499           partNum[index[i]] = part0;
01500         }
01501       }
01502       else{
01503         for (lno_t i=0; i < numLocalCoords; i++){
01504           partNum[i] = part0;
01505         }
01506       }
01507       
01508       imbalance = 0.0;       // perfect
01509       env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01510       env->debug(DETAILED_STATUS, string("Exiting serialRCB"));
01511       return;
01512     }
01513   }
01514 
01516   // Divide into balanced left and right halves.
01517 
01518   scalar_t totalLeft=0, totalRight=0;
01519   lno_t localCountLeft=0, localCountRight=0;
01520   unsigned char *flags = new unsigned char [numLocalCoords];
01521   env->localMemoryAssertion(__FILE__, __LINE__, numLocalCoords, flags) ;
01522   ArrayRCP<unsigned char> lrflags(flags, 0, numLocalCoords, true);
01523 
01524   try{
01525     BSPfindCut<mvector_t>( env, comm,
01526       params, numTestCuts, tolerance,
01527       cutDimension, coordDim, vectors, index,
01528       fractionLeft.view(0, weightDim), uniformWeights.view(0, weightDim),
01529       minCoord, maxCoord,
01530       cutValue, lrflags.view(0, numLocalCoords),
01531       totalLeft, totalRight, localCountLeft, imbalance);
01532   }
01533   Z2_FORWARD_EXCEPTIONS
01534 
01535   if (firstCall)
01536     env->memory("serial RCB start");
01537 
01539   // Adjust indices for left half and right half
01540 
01541  localCountRight = numLocalCoords - localCountLeft;
01542 
01543   if (localCountLeft){
01544 
01545     lno_t *newIndex = new lno_t [localCountLeft];
01546     env->localMemoryAssertion(__FILE__, __LINE__, localCountLeft, newIndex);
01547     ArrayView<lno_t> leftIndices(newIndex, localCountLeft);
01548 
01549     if (useIndices){
01550       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01551         if (lrflags[i] == leftFlag)
01552           newIndex[ii++] = index[i];
01553     }
01554     else{
01555       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01556         if (lrflags[i] == leftFlag)
01557           newIndex[ii++] = i;
01558     }
01559 
01560     partId_t newPart1 = part0 + numPartsLeftHalf - 1;
01561 
01562     env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01563 
01564     serialRCB<mvector_t, Adapter>(env, depth+1, params, numTestCuts, tolerance, 
01565       coordDim, vectors, leftIndices,
01566       uniformWeights.view(0, weightDim), solution,
01567       part0, newPart1, partNum);
01568 
01569     env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);
01570 
01571     delete [] newIndex;
01572   }
01573 
01574 
01575   if (localCountRight){
01576 
01577     lno_t *newIndex = new lno_t [localCountRight];
01578     env->localMemoryAssertion(__FILE__, __LINE__, localCountRight, newIndex);
01579     ArrayView<lno_t> rightIndices(newIndex, localCountRight);
01580 
01581     if (useIndices){
01582       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01583         if (lrflags[i] == rightFlag){
01584           newIndex[ii++] = index[i];
01585         }
01586     }
01587     else{
01588       for (lno_t i=0, ii=0; i < numLocalCoords; i++)
01589         if (lrflags[i] == rightFlag)
01590           newIndex[ii++] = i;
01591     }
01592 
01593     partId_t newPart0 = part0 + numPartsLeftHalf;
01594 
01595     env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01596 
01597     serialRCB<mvector_t, Adapter>(env, depth+1, params, numTestCuts, tolerance, 
01598       coordDim, vectors, rightIndices,
01599       uniformWeights.view(0, weightDim), solution,
01600       newPart0, part1, partNum);
01601 
01602     env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);
01603 
01604     delete [] newIndex;
01605   }
01606 
01607   env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
01608   env->debug(DETAILED_STATUS, string("Exiting serialRCB"));
01609 }
01610 
01611 }// namespace Zoltan2
01612 
01613 #endif
01614