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