Zoltan2_AlgRCB_methods.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//   Zoltan2: A package of combinatorial algorithms for scientific computing
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Karen Devine      (kddevin@sandia.gov)
//                    Erik Boman        (egboman@sandia.gov)
//                    Siva Rajamanickam (srajama@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

#ifndef _ZOLTAN2_ALGRCB_METHODS_HPP_
#define _ZOLTAN2_ALGRCB_METHODS_HPP_

#include <Zoltan2_PartitioningSolution.hpp>
#include <Zoltan2_XpetraTraits.hpp>
#include <Zoltan2_Metric.hpp>

#include <Tpetra_Vector.hpp>

#include <sstream>
#include <string>
#include <cmath>
#include <bitset>

#include <Teuchos_DataAccess.hpp>

namespace Teuchos{

template <typename Ordinal, typename T>
  class Zoltan2_RCBOperation : public ValueTypeReductionOp<Ordinal,T>
{
private:
  Ordinal numSum_, numMin_, numMax_;

public:
  Zoltan2_RCBOperation():numSum_(0), numMin_(0), numMax_(0){}

  Zoltan2_RCBOperation(Ordinal nsum, Ordinal nmin, Ordinal nmax):
    numSum_(nsum), numMin_(nmin), numMax_(nmax){}
  
  void reduce( const Ordinal count, const T inBuffer[], T inoutBuffer[]) const
  {
    Ordinal next=0;
    for (Ordinal i=0; i < numSum_; i++, next++)
      inoutBuffer[next] += inBuffer[next];

    for (Ordinal i=0; i < numMin_; i++, next++)
      if (inoutBuffer[next] > inBuffer[next])
        inoutBuffer[next] = inBuffer[next];

    for (Ordinal i=0; i < numMax_; i++, next++)
      if (inoutBuffer[next] < inBuffer[next])
        inoutBuffer[next] = inBuffer[next];
  }

  ArrayRCP<T> getInitializedBuffer(T minExtreme, T maxExtreme) const
  {
    Ordinal len = numSum_+numMin_+numMax_;
    T *buf = new T [len];
    if (!buf)
      throw std::bad_alloc();

    Ordinal next=0;
    for (Ordinal i=0; i < numSum_; i++)
      buf[next++] = 0;
    for (Ordinal i=0; i < numMin_; i++)
      buf[next++] = minExtreme;
    for (Ordinal i=0; i < numMax_; i++)
      buf[next++] = maxExtreme;

    ArrayRCP<T> returnBuf = arcp(buf, 0, len, true);
    return returnBuf;
  }
};
} // namespace Teuchos

namespace Zoltan2{

enum rcbParams {
  rcb_balanceCount,            
  rcb_balanceWeight,          
  rcb_minTotalWeight,      
  rcb_minMaximumWeight,  
  rcb_balanceTotalMaximum, 
  rcb_averageCuts,          
  rcb_rectilinear,    
  rcb_multiplePartSizeSpecs,  
  NUM_RCB_PARAMS
};

enum leftRightFlag{
  unsetFlag = 0xfd,       
  leftFlag = 0xfe,        
  rightFlag = 0xff        
};

// TODO: RCB has several helper methods.
// Do we want a convention for naming algorithm methods?  Do
// we want sub-namespaces for algorithms? 

template <typename Adapter>
  void getFractionLeft(
    const RCP<const Environment> &env,
    typename Adapter::part_t part0,
    typename Adapter::part_t part1,
    const RCP<PartitioningSolution<Adapter> > &solution,
    ArrayRCP<double> &fractionLeft,
    typename Adapter::part_t &numPartsLeftHalf)
{
  env->debug(DETAILED_STATUS, "Entering getFractionLeft");

  typedef typename Adapter::part_t part_t;

  part_t numParts = part1 - part0 + 1;
  // TODO In LDRD, substitute call to machine model for
  // TODO computation of numPartsLeftHalf below.
  numPartsLeftHalf = numParts / 2;
  part_t left0 = part0;   // First part in left half
  part_t left1 = left0 + numPartsLeftHalf - 1;  // Last part in left half
  part_t right0 = left1 + 1;  // First part in right half
  part_t right1 = part1;  // Last part in right half

  int nVecs = solution->getNumberOfCriteria();
  fractionLeft = arcp(new double [nVecs], 0, nVecs);

  for (int widx=0; widx<nVecs; widx++){
    if (solution->criteriaHasUniformPartSizes(widx)){
      fractionLeft[widx] = double(numPartsLeftHalf) / double(numParts);
    }
    else{
      fractionLeft[widx] = 0;
      for(int partId=left0; partId <= left1; partId++){
        fractionLeft[widx] += solution->getCriteriaPartSize(widx, partId);
      }
      double total = fractionLeft[widx];
      for(int partId=right0; partId <= right1; partId++){
        total += solution->getCriteriaPartSize(widx, partId);
      }
      fractionLeft[widx] /= total;
    }
  }
  env->debug(DETAILED_STATUS, "Exiting getFractionLeft");
}

template <typename mvector_t>
  void getCutDimension(
    const RCP<const Environment> &env,
    const RCP<Comm<int> > &comm,
    int coordDim,
    const RCP<mvector_t> &vectors,
    ArrayView<typename mvector_t::local_ordinal_type> index,
    int &dimension,                        // output
    typename mvector_t::scalar_type &minCoord,      // output
    typename mvector_t::scalar_type &maxCoord)      // output
{
  env->debug(DETAILED_STATUS, "Entering getCutDimension");
  typedef typename mvector_t::scalar_type scalar_t;
  typedef typename mvector_t::local_ordinal_type lno_t;

  int nprocs = comm->getSize();
  bool useIndices = index.size() > 0;
  lno_t numLocalCoords = 0;

  if (useIndices)
    numLocalCoords = index.size();
  else
    numLocalCoords = vectors->getLocalLength();

  Array<scalar_t> spans(coordDim*2);
  int next = 0;

  if (useIndices){
    if (numLocalCoords > 0){
      for (int dim=0; dim < coordDim; dim++){
        const scalar_t *val = vectors->getData(dim).getRawPtr();
        scalar_t v = val[index[0]];
        scalar_t min = v;
        scalar_t max = v;
        for (lno_t i=1; i < numLocalCoords; i++){
          v = val[index[i]];
          if (v > max)
            max = v;
          else if (v < min)
            min = v;
        }
  
        spans[next++] = min;
        spans[next++] = max * -1.0;
      }
    }
    else{
      for (int dim=0; dim < coordDim; dim++){
        spans[next++] = std::numeric_limits<scalar_t>::max();
        spans[next++] = std::numeric_limits<scalar_t>::min() * -1.0;
      }
    }
  }
  else{
    for (int dim=0; dim < coordDim; dim++){
      const scalar_t *val = vectors->getData(dim).getRawPtr();
      std::pair<scalar_t, scalar_t> minMax = 
        z2LocalMinMax<scalar_t>(val, numLocalCoords);
      spans[next++] = minMax.first;
      spans[next++] = minMax.second * -1.0;
    }
  }

  if (nprocs > 1){
    Array<scalar_t> globalValues(coordDim*2);
    try{
      reduceAll<int, scalar_t>(*comm, Teuchos::REDUCE_MIN, coordDim*2,
        spans.getRawPtr(), globalValues.getRawPtr());
    }
    Z2_THROW_OUTSIDE_ERROR(*env)

    for (int i=0; i < coordDim*2; i++)
      spans[i] = globalValues[i];
  }

  scalar_t maxSpan = 0;
  dimension = -1;
  next = 0;

  for (int dim=0; dim < coordDim; dim++){
    scalar_t min = spans[next++];
    scalar_t max = spans[next++] * -1.0;
    scalar_t newSpan = max - min;
    if (newSpan > maxSpan){
      maxSpan = newSpan;
      dimension = dim;
      minCoord = min;
      maxCoord = max;
    }
  }
  env->debug(DETAILED_STATUS, "Exiting getCutDimension");
}

template <typename mvector_t>
  void migrateData(
    const RCP<const Environment> &env, const RCP<Comm<int> > &comm,
    ArrayView<unsigned char> lrflags,
    RCP<mvector_t> &vectors,    // on return is the new data
    int &leftNumProcs)          // on return is num procs with left data
{
  env->debug(DETAILED_STATUS, "Entering migrateData");
  typedef typename mvector_t::scalar_type scalar_t;
  typedef typename mvector_t::local_ordinal_type lno_t;
  typedef typename mvector_t::global_ordinal_type gno_t;

  int nprocs = comm->getSize();
  size_t nobj = vectors->getLocalLength();
  size_t nGlobalObj = vectors->getGlobalLength();

  env->localBugAssertion(__FILE__, __LINE__, "migrateData input", 
    nprocs>1 && size_t(lrflags.size())==nobj, DEBUG_MODE_ASSERTION);

  gno_t myNumLeft= 0, numLeft;
  for (size_t i=0; i < nobj; i++)
    if (lrflags[i] == leftFlag)
      myNumLeft++;

  try{
    reduceAll<int, gno_t>(
      *comm, Teuchos::REDUCE_SUM, 1, &myNumLeft, &numLeft);
  }
  Z2_THROW_OUTSIDE_ERROR(*env)
    
  scalar_t leftFraction = scalar_t(numLeft)/scalar_t(nGlobalObj);
  leftNumProcs = static_cast<int>(floor((nprocs*leftFraction)+.5));

  if (leftNumProcs == 0)
    leftNumProcs++;
  else if (leftNumProcs == nprocs)
    leftNumProcs--;

  // Get a list of my new global numbers.

  Array<int> sendCount(nprocs, 0);
  Array<int> recvCount(nprocs, 0);
  Array<gno_t> sendBuf(nobj, 0);

  if (nobj > 0){
    int *procId = new int [nobj];
    env->localMemoryAssertion(__FILE__, __LINE__, nobj, procId) ;
    int leftProc0 = 0;
    int rightProc0 = leftProc0 + leftNumProcs;
  
    int nextLeftProc = leftProc0;
    int nextRightProc = rightProc0;
    int *p = procId;
  
    for (size_t i=0; i < nobj; i++){
      if (lrflags[i] == leftFlag){
        if (nextLeftProc == rightProc0)
          nextLeftProc = leftProc0;
  
        sendCount[nextLeftProc]++;
        *p++ = nextLeftProc++;
      }
      else{
        if (nextRightProc == nprocs)
          nextRightProc = rightProc0;
        sendCount[nextRightProc]++;
  
        *p++ = nextRightProc++;
      }
    }
  
    gno_t *sendOffset = new gno_t [nprocs];
    env->localMemoryAssertion(__FILE__, __LINE__, nprocs, sendOffset) ;
    sendOffset[0] = 0;
    for (int i=0; i < nprocs-1; i++)
      sendOffset[i+1] = sendOffset[i] + sendCount[i];

    ArrayView<const gno_t> gnoList = vectors->getMap()->getNodeElementList();

    for (size_t i=0; i < nobj; i++){
      int proc = procId[i];
      lno_t offset = sendOffset[proc]++;

      sendBuf[offset] = gnoList[i];
    }

    delete [] sendOffset;
    delete [] procId;
  }

  ArrayRCP<gno_t> recvBuf;

  try{
    AlltoAllv<gno_t>(*comm, *env,
      sendBuf(), sendCount(),
      recvBuf, recvCount());
  }
  Z2_FORWARD_EXCEPTIONS

  sendCount.clear();
  sendBuf.clear();

  // Migrate the multivector of data.

  gno_t numMyNewGnos = 0;
  for (int i=0; i < nprocs; i++)
    numMyNewGnos += recvCount[i];

  recvCount.clear();

  RCP<const mvector_t> newMultiVector;
  RCP<const mvector_t> constInput = rcp_const_cast<const mvector_t>(vectors);

  try{
    newMultiVector = XpetraTraits<mvector_t>::doMigration(
      constInput, numMyNewGnos, recvBuf.getRawPtr());
  }
  Z2_FORWARD_EXCEPTIONS

  vectors = rcp_const_cast<mvector_t>(newMultiVector);
  env->memory("Former problem data replaced with new data");
  env->debug(DETAILED_STATUS, "Exiting migrateData");
}

template <typename lno_t, typename scalar_t>
  scalar_t getCoordWeight(lno_t id, 
    multiCriteriaNorm mcnorm, 
    ArrayView<StridedData<lno_t, scalar_t> > weights)
{
  scalar_t wgtValue = 1.0;
  size_t nWeightsPerCoord = weights.size();
    
  if (nWeightsPerCoord > 1){
    Array<scalar_t> coordWeights(nWeightsPerCoord, 1.0);
    for (size_t widx=0; widx < nWeightsPerCoord; widx++){
      coordWeights[widx] = weights[widx][id];
    }

    wgtValue = normedWeight<scalar_t>(coordWeights.view(0,nWeightsPerCoord), mcnorm);
  }
  else if (nWeightsPerCoord > 0){
    wgtValue = weights[0][id];
  }

  return wgtValue;
}

template <typename scalar_t>
  bool emptyPartsCheck( const RCP<const Environment> &env,
   const ArrayView<double> fractionLeft, 
   scalar_t minCoord, scalar_t maxCoord,
   leftRightFlag &lrf,
   scalar_t &cutValue)
{
  env->debug(DETAILED_STATUS, "Entering emptyPartsCheck");
  // initialize return values
  lrf = leftFlag;
  cutValue = 0.0;

  size_t nVecs = fractionLeft.size();
  int numEmptyRight = 0, numEmptyLeft = 0;

  for (size_t i=0; i < nVecs; i++){
    if (fractionLeft[i] == 0.0)
      numEmptyLeft++;
    else if (fractionLeft[i] == 1.0)
      numEmptyRight++;
  }

  if (numEmptyRight + numEmptyLeft == 0)
    return false;

  env->localInputAssertion(__FILE__, __LINE__,
    "partitioning not solvable - conflicting part size criteria",
    (numEmptyRight==0) || (numEmptyLeft==0), BASIC_ASSERTION);

  if (numEmptyLeft > 0){
    lrf = leftFlag;
    cutValue = minCoord;
  }
  else{
    lrf = rightFlag;
    cutValue = maxCoord;
  }

  env->debug(DETAILED_STATUS, "Exiting emptyPartsCheck");
  return true;
}

template <typename lno_t, typename gno_t, typename scalar_t>
  void testCoordinatesOnRightBoundary(
    const RCP<const Environment> &env,
    const RCP<Comm<int> > &comm, 
    scalar_t totalWeightLeft,
    scalar_t targetLeftScalar,
    int cutLocation,
    ArrayView<scalar_t> localSums,
    ArrayView<scalar_t> globalSums,
    ArrayView<lno_t> index,
    ArrayView<StridedData<lno_t, scalar_t> > weights,
    multiCriteriaNorm mcnorm,
    ArrayView<unsigned char> lrFlags,
    scalar_t &globalWeightMovedRight)   // output
{
  env->debug(DETAILED_STATUS, "Entering testCoordinatesOnRightBoundary");
  int nprocs = comm->getSize();
  int rank = comm->getRank();

  globalWeightMovedRight = 0.0;

  scalar_t localBoundarySum = localSums[cutLocation];

  scalar_t total = totalWeightLeft;
  for (int i=0; i <= cutLocation; i++)
    total += globalSums[i];

  scalar_t totalMoveRight = total - targetLeftScalar;
  scalar_t localMoveRight = localBoundarySum;
  scalar_t actualWeightMovedRight = 0.0;
  
  Array<scalar_t> scansum(nprocs+1, 0.0);
  Teuchos::gatherAll<int, scalar_t>(*comm, 1, 
    &localBoundarySum, nprocs, scansum.getRawPtr()+1);
  for (int i=2; i<=nprocs; i++)
    scansum[i] += scansum[i-1];

  if (localBoundarySum > 0.0){

    scalar_t sumMine = scansum[rank]; // sum of ranks preceding me
    scalar_t diff = scansum[nprocs] - sumMine;
    localMoveRight = 0;

    if (diff <= totalMoveRight)
      localMoveRight = localBoundarySum;  // all
    else{
      scalar_t leftPart = diff - totalMoveRight;
      if (leftPart < localBoundarySum)
        localMoveRight = localBoundarySum - leftPart;
    }
  }

  if (localMoveRight > 0.0){

    bool moveAll =  (localMoveRight >= localBoundarySum);

    if (moveAll){
      actualWeightMovedRight = localBoundarySum;
      for (lno_t i=0; i < lrFlags.size(); i++){
        if (lrFlags[i] == cutLocation){
          lrFlags[i] = rightFlag;
        }
      }
    }
    else{
      int nWeightsPerCoord = weights.size();
      int useIndices = index.size() > 0;

      for (lno_t i=0; i < lrFlags.size(); i++){
        if (lrFlags[i] == cutLocation){
          lrFlags[i] = rightFlag;

          lno_t idx = (useIndices ? index[i] : i);

          actualWeightMovedRight += getCoordWeight<lno_t, scalar_t>(
            idx, mcnorm, weights.view(0, nWeightsPerCoord));

          if (actualWeightMovedRight >= localMoveRight)
            break;
        }
      } // next coordinate
    }
  }

  try{
    reduceAll<int, scalar_t>(
      *comm, Teuchos::REDUCE_SUM, 1, &actualWeightMovedRight,
      &globalWeightMovedRight);
  }
  Z2_THROW_OUTSIDE_ERROR(*env)

  env->debug(DETAILED_STATUS, "Exiting testCoordinatesOnRightBoundary");
  return;
}

template <typename mvector_t>
  void BSPfindCut(
    const RCP<const Environment> &env,
    const RCP<Comm<int> > &comm,
    const std::bitset<NUM_RCB_PARAMS> &params,
    int numTestCuts,
    typename mvector_t::scalar_type tolerance,
    int cutDim,
    int coordDim,
    int nWeightsPerCoord,
    const RCP<mvector_t> &vectors,
    ArrayView<typename mvector_t::local_ordinal_type> index,
    ArrayView<double> fractionLeft,
    typename mvector_t::scalar_type coordGlobalMin,
    typename mvector_t::scalar_type coordGlobalMax,
    typename mvector_t::scalar_type &cutValue,         // output
    ArrayView<unsigned char> lrFlags,                  // output
    typename mvector_t::scalar_type &totalWeightLeft,  // output
    typename mvector_t::scalar_type &totalWeightRight, // output
    typename mvector_t::local_ordinal_type &localCountLeft, // output
    typename mvector_t::scalar_type &imbalance)        // output
{
  env->debug(DETAILED_STATUS, "Entering BSPfindCut");

  // initialize output
  bool useIndices = index.size() > 0;
  int numAllCoords = vectors->getLocalLength();
  int numCoords = (useIndices ? index.size() : numAllCoords);
  cutValue = totalWeightLeft = totalWeightRight = imbalance = 0.0;
  localCountLeft = 0;

  for (int i=0; i < numCoords; i++)
    lrFlags[i] = unsetFlag;

  typedef typename mvector_t::scalar_type scalar_t;
  typedef typename mvector_t::local_ordinal_type lno_t;
  typedef typename mvector_t::global_ordinal_type gno_t;
  typedef StridedData<lno_t, scalar_t> input_t;

  bool multiplePartSizeSpecs = params.test(rcb_multiplePartSizeSpecs);
  bool rectilinear = params.test(rcb_rectilinear);
  bool averageCuts = params.test(rcb_averageCuts);

  // A coordinate is considered to be on a cut if it is within
  // this distance of the cut.

  double epsilon = (coordGlobalMax - coordGlobalMin) * 10e-9;

  // Find the coordinate values and weights.

  int nVecs = fractionLeft.size();

  if (env->getDebugLevel() >= DETAILED_STATUS){
    std::ostringstream info;
    info << "Num weights " << nWeightsPerCoord << ", Fraction left:";
    for (int i=0; i < nVecs; i++)
      info << " " << fractionLeft[i];
    info << endl << "Dimension " << cutDim << " [";
    info << coordGlobalMin << ", " << coordGlobalMax << "]";
    info << endl << "# test cuts " << numTestCuts;
    info << ", tolerance " << tolerance << std::endl;
    env->debug(DETAILED_STATUS, info.str());
  }

  const scalar_t *coordValue = vectors->getData(cutDim).getRawPtr();

  input_t *wgtinfo = new input_t [nWeightsPerCoord];
  env->localMemoryAssertion(__FILE__, __LINE__, nWeightsPerCoord, wgtinfo);
  ArrayRCP<input_t> weight(wgtinfo, 0, nWeightsPerCoord, true);

  for (int widx = 0, cidx=coordDim; widx < nWeightsPerCoord; widx++){
    weight[widx] = input_t(vectors->getData(cidx++), 1);
  }

  // Multicriteria norm

  multiCriteriaNorm mcnorm = normBalanceTotalMaximum;

  if (params.test(rcb_minMaximumWeight))
    mcnorm = normMinimizeMaximumWeight;
  else if (params.test(rcb_minTotalWeight))
    mcnorm = normMinimizeTotalWeight;
  
  // Goal is globally find one cut that comes close to leaving
  // partSizeLeft*totalWeight on the left side.

  Epetra_SerialDenseVector partSizeLeft(::Teuchos::View, fractionLeft.getRawPtr(), nVecs);
  
  // Where do we make the first test cuts?
  //
  //   min     1     2     3     max
  //
  // Space is [min, max].  
  // If there are three test cuts: 1, 2 and 3.
  //   4 regions: [min,1] [1,2] [2,3] [3,max].
  //   5 boundaries: min, 1, 2, 3, and max.

  int numRegions = numTestCuts + 1;
  int numBoundaries = numTestCuts + 2;
  int endBoundary = numBoundaries - 1;
  std::vector<scalar_t> boundaries(numBoundaries);
  std::vector<scalar_t> searchBoundaries(numBoundaries);

  int numSums = numBoundaries+numRegions;

  Teuchos::Zoltan2_RCBOperation<int, scalar_t> reductionOp(
    numSums,      // number of sums
    numRegions,   // number of mins
    numRegions);  // number of maxes

  typename std::vector<scalar_t>::iterator foundCut;

  bool done=false;
  bool fail=false;
  scalar_t min = coordGlobalMin;
  scalar_t max = coordGlobalMax;
  lno_t numRemaining = numCoords;
  lno_t prevNumRemaining = numCoords;
  size_t numGlobalPoints = vectors->getGlobalLength();
  size_t sanityCheck = numGlobalPoints;

  double totalWeight = 0;
  double targetLeftScalar = 0;
  double targetLeftNorm = 0;
  Epetra_SerialDenseVector targetLeftVector(nVecs);

  while (!done && !fail && sanityCheck--){

    // Create regions into which coordinates will be placed.
    scalar_t diff = (max - min) / numRegions;
    boundaries[0] = min;
    boundaries[endBoundary] = max;
    for (int i=0; i < endBoundary; i++){
      searchBoundaries[i+1] = boundaries[i+1] = boundaries[i] + diff;
    }

    // Move ends slightly so we catch points on boundary.
    searchBoundaries[0] = min - epsilon;
    searchBoundaries[endBoundary] = max + epsilon;

    // Save region and boundary sums, and region min and max.
 
    ArrayRCP<scalar_t> localSums = reductionOp.getInitializedBuffer(
      searchBoundaries[endBoundary], searchBoundaries[0]);
 
    ArrayRCP<scalar_t> globalSums = reductionOp.getInitializedBuffer(
      searchBoundaries[endBoundary], searchBoundaries[0]);

    scalar_t *sums = localSums.getRawPtr();
    scalar_t *regionMin = sums + numSums;
    scalar_t *regionMax = regionMin + numRegions;

    if (numRemaining > 0){

      // Assign each of my points to a region.
      // lower_bound() finds the first cut f, such that f >= coordValue[i].
      // So for now, objects that are on the cut boundary go into the
      // region on the "left" side.

      for (lno_t i=0; i < numCoords; i++){
  
        if (lrFlags[i] != unsetFlag)
          continue;

        int inRegion = 0;
        int idx = (useIndices ? index[i] : i);
        scalar_t value = coordValue[idx];

        if (numRegions > 2){
       
          foundCut = std::lower_bound(
            searchBoundaries.begin(), searchBoundaries.end(), value);
          
          env->localBugAssertion(__FILE__, __LINE__, "search cuts", 
            foundCut != searchBoundaries.end(), BASIC_ASSERTION);

          inRegion = foundCut - searchBoundaries.begin() - 1;        
        }
        else{
          if (value <= boundaries[1])
            inRegion = 0;
          else
            inRegion = 1;
        }

        int sumIdx = 1 + inRegion*2;

        if (value >= boundaries[inRegion+1]-epsilon){  
          // "on" right boundary of this region
          sumIdx++;
        }
        else if (inRegion==0 && (value < min + epsilon)){
          // in region 0 but far left boundary
          sumIdx--;
        }

        lrFlags[i] = (unsigned char)sumIdx;

        if (value < regionMin[inRegion])
          regionMin[inRegion] = value;
        if (value > regionMax[inRegion])
          regionMax[inRegion] = value;

        sums[sumIdx] += getCoordWeight<lno_t, scalar_t>(idx, 
              mcnorm, weight.view(0, nWeightsPerCoord));

      } // next coord
    }

    try{
      reduceAll<int, scalar_t>(*comm, reductionOp,
        localSums.size(), localSums.getRawPtr(), globalSums.getRawPtr());
    }
    Z2_THROW_OUTSIDE_ERROR(*env)

    sums = globalSums.getRawPtr();
    regionMin = sums + numSums;
    regionMax = regionMin + numRegions;

    if (env->getDebugLevel() >= DETAILED_STATUS){
      std::ostringstream info;
      info << "  Region " << min << " - " << max << std::endl;
      info << "  Remaining to classify: " << numRemaining << std::endl;
      info << "  Boundaries: ";
      for (int i=0; i < numBoundaries; i++)
        info << boundaries[i] << " ";
      info << endl << "  For searching: ";
      for (int i=0; i < numBoundaries; i++)
        info << searchBoundaries[i] << " ";
      info << endl << "  Global sums: ";
      double sumTotal=0;
      for (int i=0; i < numSums; i++){
        sumTotal += sums[i];
        info << sums[i] << " ";
      }
      info << " total: " << sumTotal << endl;
      env->debug(DETAILED_STATUS, info.str());
    }

    if (totalWeight == 0){   // first time through only

      for (int i=0; i < numSums; i++)
        totalWeight += sums[i];

      partSizeLeft.Scale(totalWeight);
      targetLeftVector = partSizeLeft;

      targetLeftScalar = targetLeftVector[0];
      targetLeftNorm = targetLeftVector.Norm2();
      totalWeightLeft = 0;
      totalWeightRight = 0;
    }

    int cutLocation=0;
    scalar_t testDiff=0, prevTestDiff=0, target=0;

    if (multiplePartSizeSpecs){
      // more complex: if we have multiple weights
      //   and the part sizes requested
      //   for each weight index differ, then we may not
      //   be able to reach the imbalance tolerance.
      //
      // TODO: discuss how to (or whether to) handle this case.
      // 
      // For now we look at this imbalance:
      // 
      //    |target - actual|^2 / |target|^2
  
      target = targetLeftNorm;
      Epetra_SerialDenseVector testVec(nVecs);
      for (int i=0; i < nVecs; i++)
        testVec[i] = totalWeightLeft;
      Epetra_SerialDenseVector diffVec = testVec;
      diffVec.Scale(-1.0);
      diffVec += targetLeftVector;

      testDiff = diffVec.Norm2(); // imbalance numerator
      prevTestDiff = testDiff;
      cutLocation= -1;

      while (++cutLocation< numSums){

        for (int i=0; i < nVecs; i++)
          testVec[i] += sums[cutLocation];
  
        diffVec = testVec;
        diffVec.Scale(-1.0);
        diffVec += targetLeftVector;
  
        testDiff = diffVec.Norm2();
        
        if (testDiff >= target)
          break;

        prevTestDiff = testDiff;
      }
    }
    else{    // the part sizes for each weight index are the same
  
      target = targetLeftScalar;
      testDiff = totalWeightLeft; 
      prevTestDiff = testDiff;
      cutLocation = -1;
  
      while (++cutLocation < numSums){
     
        testDiff += sums[cutLocation];
        if (testDiff >= target)
          break;
  
        prevTestDiff = testDiff;
      }
    }

    scalar_t diffLeftCut = target - prevTestDiff;
    scalar_t diffRightCut = testDiff - target;

    if (diffLeftCut < diffRightCut){
      imbalance = diffLeftCut / target;
      if (imbalance <= tolerance){
        env->debug(DETAILED_STATUS, "  Done, tolerance met");
        done = true;
        cutLocation--;
      }
    } 
    else{
      imbalance = diffRightCut / target;
      if (imbalance <= tolerance){
        env->debug(DETAILED_STATUS, "  Done, tolerance met");
        done = true;
     }
    }

    bool cutLocIsRegion = (cutLocation % 2 == 1);
    bool cutLocIsBoundary = !cutLocIsRegion;

    if (env->getDebugLevel() >= DETAILED_STATUS){
      std::ostringstream info;
      info << "  Best cut location: " << cutLocation;
      if (cutLocIsRegion) info << " just after a region." << std::endl;
      else info << " just after a boundary." << std::endl;
      env->debug(DETAILED_STATUS, info.str());
    }

    if (!done && cutLocIsBoundary){

      done = true;    // can not divide space any more

      env->debug(DETAILED_STATUS, "  Done, cutting at a region boundary");

      if (rectilinear){
        // Can not divide boundary points into two
        // different regions to achieve balance.
        fail = true;
      }
      else {
        // Maybe by moving some of the points right, we can
        // obtain a better balance.  If so, the lrFlag for
        // the points moved right will be updated here.

        scalar_t globalWeightMovedRight(0);

        testCoordinatesOnRightBoundary<lno_t, gno_t, scalar_t>(
          env, comm, totalWeightLeft, targetLeftScalar, cutLocation, 
          localSums.view(0, numSums), globalSums.view(0, numSums),
          index, weight.view(0, nWeightsPerCoord), mcnorm, lrFlags,
          globalWeightMovedRight);

        scalar_t newSum = testDiff - globalWeightMovedRight;
        globalSums[cutLocation] -= globalWeightMovedRight;

        if (newSum > target)
          imbalance = (target - newSum) / target;
        else
          imbalance = (newSum - target) / target;

        if (imbalance > tolerance)
          fail = true;
      }
    }

    int rightmostLeftNum=0, leftmostRightNum=0;

    if (!done){
      // Best cut is following a region.  Narrow down the boundaries.

      int regionNumber = (cutLocation - 1) / 2;

      min = regionMin[regionNumber];
      max = regionMax[regionNumber];

      rightmostLeftNum = cutLocation - 1;
      leftmostRightNum = cutLocation + 1;
    }
    else{
      rightmostLeftNum = cutLocation;
      leftmostRightNum = cutLocation + 1;

      if (cutLocIsRegion && averageCuts){
        int regionNumber = (cutLocation - 1) / 2;
        cutValue = (
          boundaries[regionNumber+1] + // boundary to right
          regionMax[regionNumber])     // rightmost point in region
          / 2.0;
      }
    }

    for (int i=0; i <= rightmostLeftNum; i++){
      totalWeightLeft += globalSums[i];
    }

    for (lno_t i=0; i < numCoords; i++){
      if (lrFlags[i] != leftFlag && lrFlags[i] != rightFlag){
        if (lrFlags[i] <= rightmostLeftNum){
          lrFlags[i] = leftFlag;
          numRemaining--;
          localCountLeft++;
        }
        else if (lrFlags[i] >= leftmostRightNum){
          lrFlags[i] = rightFlag;
          numRemaining--;
        }
        else{
          lrFlags[i] = unsetFlag;   // still to be determined
        }
      }
    }

    if (env->getDebugLevel() >= VERBOSE_DETAILED_STATUS && numCoords < 100){
      // For large numCoords, building this message
      // takes an extraordinarily long time.
      std::ostringstream ossLeft;
      std::ostringstream ossRight;
      ossLeft << "left: ";
      ossRight << "right: ";
      for (lno_t i=0; i < numCoords; i++){
        if (lrFlags[i] == unsetFlag)
          continue;
        lno_t idx = (useIndices ? index[i] : i);
        scalar_t val = coordValue[idx];
        if (lrFlags[i] == leftFlag)
          ossLeft << val << " ";
        else if (lrFlags[i] == rightFlag)
          ossRight << val << " ";
        else 
          env->localBugAssertion(__FILE__, __LINE__, 
            "left/right flags", false, BASIC_ASSERTION);
      }
      std::ostringstream msg;
      msg << ossLeft.str() << endl << ossRight.str() << std::endl;
      env->debug(VERBOSE_DETAILED_STATUS, msg.str());
    }

    prevNumRemaining = numRemaining;
  }  // while !done

  totalWeightRight = totalWeight - totalWeightLeft;

  if (fail)
    env->debug(BASIC_STATUS, "Warning: tolerance not achieved in sub-step");

  // TODO: If fail the warn that tolerance was not met.

  env->globalInputAssertion(__FILE__, __LINE__, "partitioning not solvable",
    done, DEBUG_MODE_ASSERTION, comm);

  env->memory("End of bisection");

  if (env->getDebugLevel() >= DETAILED_STATUS){
    std::ostringstream oss;
    oss << "Exiting BSPfindCut, ";
    oss << "# iterations: " << numGlobalPoints - sanityCheck;    
    env->debug(DETAILED_STATUS, oss.str());
  }

  return;
}

template <typename mvector_t, typename Adapter>
  void determineCut(
    const RCP<const Environment> &env,
    const RCP<Comm<int> > &comm,
    const std::bitset<NUM_RCB_PARAMS> &params,
    int numTestCuts,
    typename mvector_t::scalar_type tolerance,
    int coordDim, 
    int nWeightsPerCoord,
    const RCP<mvector_t> &vectors,
    multiCriteriaNorm mcnorm,
    const RCP<PartitioningSolution<Adapter> > &solution,
    typename Adapter::part_t part0, 
    typename Adapter::part_t part1,
    ArrayView<unsigned char> lrflags,  // output
    int &cutDimension,                  // output
    typename mvector_t::scalar_type &cutValue,   // output
    typename mvector_t::scalar_type &imbalance,  // output
    typename Adapter::part_t &numPartsLeftHalf,                  // output
    typename mvector_t::scalar_type &weightLeftHalf,  // output
    typename mvector_t::scalar_type &weightRightHalf  // output
    )
{
  env->debug(DETAILED_STATUS, "Entering determineCut");
  typedef typename Adapter::part_t part_t;
  typedef typename mvector_t::scalar_type scalar_t;
  typedef typename mvector_t::local_ordinal_type lno_t;
  typedef StridedData<lno_t, scalar_t> input_t;

  lno_t numLocalCoords = vectors->getLocalLength();

  // initialize return values
  cutDimension = 0;
  cutValue = imbalance = weightLeftHalf = weightRightHalf = 0.0;
  numPartsLeftHalf = 0;

  // Pick a cut direction.

  scalar_t globalMinCoord, globalMaxCoord;
  ArrayView<lno_t> emptyIndex;

  getCutDimension<mvector_t>(env, comm, coordDim, vectors, emptyIndex,
    cutDimension, globalMinCoord, globalMaxCoord);

  // Compute part sizes for the two parts.

  ArrayRCP<double> fractionLeft;

  getFractionLeft<Adapter>(env, part0, part1, solution,
    fractionLeft, numPartsLeftHalf);
  int nVecs = fractionLeft.size();

  // Special case of empty left or empty right.

  leftRightFlag lrf;

  bool emptyParts = emptyPartsCheck(env,
    fractionLeft.view(0, nVecs), // input
    globalMinCoord, globalMaxCoord,  // input
    lrf, cutValue);                  // output

  if (emptyParts){
    
    for (lno_t i=0; i < lrflags.size(); i++)
      lrflags[i] = lrf;

    imbalance = 0.0;                // perfect balance
    scalar_t totalWeight = 0.0;

    if (nWeightsPerCoord == 0)
      totalWeight = numLocalCoords;
    else if (nWeightsPerCoord == 1) {
      int wgtidx = coordDim;
      const scalar_t *val = vectors->getData(wgtidx).getRawPtr();
      for (lno_t i=0; i < numLocalCoords; i++)
        totalWeight += val[i];
    }
    else {  // need to add up total normed weight
      Array<input_t> wgts(nWeightsPerCoord);
      int wgtidx = coordDim;
      for (int i=0; i < nWeightsPerCoord; i++)
        wgts[i] = input_t(vectors->getData(wgtidx++), 1);

      part_t numParts, numNonemptyParts;
      ArrayRCP<MetricValues<scalar_t> > metrics;
      ArrayRCP<scalar_t> weightSums;
    
      globalSumsByPart<scalar_t, unsigned char, lno_t, part_t>(
        env, comm, lrflags, nWeightsPerCoord,
        wgts.view(0, nWeightsPerCoord), mcnorm,
        numParts, numNonemptyParts, metrics, weightSums);

      totalWeight = weightSums[1];  // sum normed weight
    }

    if (lrf == leftFlag){
      numPartsLeftHalf = part1 - part0 + 1;
      weightLeftHalf = totalWeight;
      weightRightHalf = 0;
    }
    else{
      numPartsLeftHalf = 0;
      weightRightHalf = totalWeight;
      weightLeftHalf = 0;
    }

    env->debug(DETAILED_STATUS, "Exiting determineCut");
    return;
  }

  // Divide the coordinates into balanced left and right
  // halves.

  ArrayView<lno_t> emptyIndices;
  lno_t localCountLeft;

  try{
    BSPfindCut<mvector_t>( env, comm,
      params, numTestCuts, tolerance,
      cutDimension, coordDim, nWeightsPerCoord, vectors, emptyIndices,
      fractionLeft.view(0, nVecs),
      globalMinCoord, globalMaxCoord,
      cutValue, lrflags.view(0, numLocalCoords),
      weightLeftHalf, weightRightHalf, localCountLeft, imbalance);
  }
  Z2_FORWARD_EXCEPTIONS
  env->debug(DETAILED_STATUS, "Exiting determineCut");
}


template <typename mvector_t, typename Adapter>
 void serialRCB(
    const RCP<const Environment> &env,
    int depth,
    const std::bitset<NUM_RCB_PARAMS> &params,
    int numTestCuts, 
    typename mvector_t::scalar_type tolerance, 
    int coordDim,
    int nWeightsPerCoord,
    const RCP<mvector_t> &vectors, 
    ArrayView<typename mvector_t::local_ordinal_type> index,
    const RCP<PartitioningSolution<Adapter> > &solution,
    typename Adapter::part_t part0, 
    typename Adapter::part_t part1,
    ArrayView<typename Adapter::part_t> partNum)   // output
{
  env->debug(DETAILED_STATUS, "Entering serialRCB");
  env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);

  typedef typename Adapter::part_t part_t;
  typedef typename mvector_t::scalar_type scalar_t;
  typedef typename mvector_t::local_ordinal_type lno_t;

  RCP<Comm<int> > comm(new Teuchos::SerialComm<int>);  

  lno_t numLocalCoords=0;
  bool useIndices;
  bool firstCall=false;

  if (index.size() == 0){
    // First time through there are no indices.
    useIndices = false;
    numLocalCoords = vectors->getLocalLength();
    firstCall = true;
  }
  else{
    useIndices = true;
    numLocalCoords = index.size();
  }

  if (env->getDebugLevel() >= DETAILED_STATUS){
    std::ostringstream info;
    info << "  Number of coordinates: " << numLocalCoords << std::endl;
    info << "  Use index: " << useIndices << std::endl;
    env->debug(DETAILED_STATUS, info.str());
  }

  // Are we done?

  if (part1 == part0){
    if (useIndices)
      for (lno_t i=0; i < numLocalCoords; i++)
        partNum[index[i]] = part0;
    else
      for (lno_t i=0; i < numLocalCoords; i++)
        partNum[i] = part0;

    env->memory("serial RCB end");
    env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
    env->debug(DETAILED_STATUS, "Exiting serialRCB");

    return;
  }

  // Pick a cut direction

  int cutDimension;
  scalar_t minCoord, maxCoord;

  try{
    getCutDimension(env, comm, coordDim, vectors, index,
      cutDimension, minCoord, maxCoord);
  }
  Z2_FORWARD_EXCEPTIONS

  // Compute relative sizes of the two halves.

  ArrayRCP<double> fractionLeft;
  part_t numPartsLeftHalf;

  try{
    getFractionLeft<Adapter>(env, part0, part1, solution,
      fractionLeft, numPartsLeftHalf);
  }
  Z2_FORWARD_EXCEPTIONS
  int nVecs = fractionLeft.size();

  // Check for special case of empty left or empty right.

  scalar_t imbalance, cutValue;  //unused for now
  leftRightFlag lrf;

  bool emptyPart = emptyPartsCheck(env, fractionLeft.view(0, nVecs), 
    minCoord, maxCoord, lrf, cutValue);

  if (emptyPart){  // continue only on the side that is not empty

    if (lrf == rightFlag)  // right is empty
      part1 = part0 + numPartsLeftHalf - 1;
    else
      part0 = part0 + numPartsLeftHalf;

    if (part0 == part1){
      if (useIndices){
        for (lno_t i=0; i < numLocalCoords; i++){
          partNum[index[i]] = part0;
        }
      }
      else{
        for (lno_t i=0; i < numLocalCoords; i++){
          partNum[i] = part0;
        }
      }
      
      imbalance = 0.0;       // perfect
      env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
      env->debug(DETAILED_STATUS, "Exiting serialRCB");
      return;
    }
  }

  // Divide into balanced left and right halves.

  scalar_t totalLeft=0, totalRight=0;
  lno_t localCountLeft=0, localCountRight=0;
  unsigned char *flags = new unsigned char [numLocalCoords];
  env->localMemoryAssertion(__FILE__, __LINE__, numLocalCoords, flags) ;
  ArrayRCP<unsigned char> lrflags(flags, 0, numLocalCoords, true);

  try{
    BSPfindCut<mvector_t>( env, comm,
      params, numTestCuts, tolerance,
      cutDimension, coordDim, nWeightsPerCoord, vectors, index,
      fractionLeft.view(0, nVecs),
      minCoord, maxCoord,
      cutValue, lrflags.view(0, numLocalCoords),
      totalLeft, totalRight, localCountLeft, imbalance);
  }
  Z2_FORWARD_EXCEPTIONS

  if (firstCall)
    env->memory("serial RCB start");

  // Adjust indices for left half and right half

 localCountRight = numLocalCoords - localCountLeft;

  if (localCountLeft){

    lno_t *newIndex = new lno_t [localCountLeft];
    env->localMemoryAssertion(__FILE__, __LINE__, localCountLeft, newIndex);
    ArrayView<lno_t> leftIndices(newIndex, localCountLeft);

    if (useIndices){
      for (lno_t i=0, ii=0; i < numLocalCoords; i++)
        if (lrflags[i] == leftFlag)
          newIndex[ii++] = index[i];
    }
    else{
      for (lno_t i=0, ii=0; i < numLocalCoords; i++)
        if (lrflags[i] == leftFlag)
          newIndex[ii++] = i;
    }

    part_t newPart1 = part0 + numPartsLeftHalf - 1;

    env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);

    serialRCB<mvector_t, Adapter>(env, depth+1, params, numTestCuts, tolerance, 
      coordDim, nWeightsPerCoord, vectors, leftIndices,
      solution, part0, newPart1, partNum);

    env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);

    delete [] newIndex;
  }


  if (localCountRight){

    lno_t *newIndex = new lno_t [localCountRight];
    env->localMemoryAssertion(__FILE__, __LINE__, localCountRight, newIndex);
    ArrayView<lno_t> rightIndices(newIndex, localCountRight);

    if (useIndices){
      for (lno_t i=0, ii=0; i < numLocalCoords; i++)
        if (lrflags[i] == rightFlag){
          newIndex[ii++] = index[i];
        }
    }
    else{
      for (lno_t i=0, ii=0; i < numLocalCoords; i++)
        if (lrflags[i] == rightFlag)
          newIndex[ii++] = i;
    }

    part_t newPart0 = part0 + numPartsLeftHalf;

    env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);

    serialRCB<mvector_t, Adapter>(env, depth+1, params, numTestCuts, tolerance, 
      coordDim, nWeightsPerCoord, vectors, rightIndices,
      solution, newPart0, part1, partNum);

    env->timerStart(MICRO_TIMERS, "serialRCB", depth, 2);

    delete [] newIndex;
  }

  env->timerStop(MICRO_TIMERS, "serialRCB", depth, 2);
  env->debug(DETAILED_STATUS, "Exiting serialRCB");
}

}// namespace Zoltan2

#endif