parallel/RCB.cpp

// -*- C++ -*-
//
//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//                                   Fehmi Cirak
//                        California Institute of Technology
//                           (C) 2004 All Rights Reserved
//
//  <LicenseText>
//
//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// EDITED BY MATT GALLOWAY FOR USE WITH DAPTA

#include "DAPTA.hpp"

#include <cmath>
#include <vector>
#include <algorithm>
#include <functional>
#include <cassert>
#include <cstdlib>

namespace DAPTA { // Define namespace DAPTA

template <typename V>
std::vector<int> RCBDecomposer<V>::decompose(unsigned numberOfPartitions, unsigned myNumber) {
   std::vector<int> myNodes;

    // head of the RCB tree
    RCBTree *head = new RCBTree(this->vertices.begin(), this->vertices.end(), 0);
        
    // maximum height of the RCB tree
    int maxLevel = (int)std::ceil(std::log((double)numberOfPartitions)/std::log(2.0));

    // RCB partitioning
    buildRCBTree(head, maxLevel);

    // assign the RCB nodes to the processors
    std::vector<PointIt> leafDataBegins, leafDataEnds;    
    leafDataBegins.reserve(numberOfPartitions*2);
    leafDataEnds.reserve(numberOfPartitions*2);
    collectLeafNodeData(head, leafDataBegins, leafDataEnds);
    
    std::vector<PointIt> collectProcBegins, collectProcEnds;
    collectProcBegins.reserve(leafDataBegins.size());
    collectProcEnds.reserve(leafDataEnds.size());
    
    // combine leafs so that there are as many subdomains as processors
    unsigned numLeafs = leafDataBegins.size();
    unsigned offSet=0;
    while (numLeafs!=numberOfPartitions) {
      collectProcBegins.push_back(leafDataBegins[offSet]);
      collectProcEnds.push_back(leafDataEnds[offSet+1]);
      offSet += 2;
      --numLeafs;
    }
    
    // just copy the remaining leafs
    for (unsigned j=offSet; j<leafDataBegins.size(); ++j) {
      collectProcBegins.push_back(leafDataBegins[j]);
      collectProcEnds.push_back(leafDataEnds[j]);
    }
    
    for (unsigned i=0; i<collectProcBegins.size(); ++i) {
      if (i!=myNumber) continue;
      for (PointIt it = collectProcBegins[i]; it!=collectProcEnds[i]; ++it) {
         myNodes.push_back((*it)->globalID);
      }
    }

    // delete the tree
    delete head;
   
   return myNodes;
}

// compute the longest axis
template <typename V>
int RCBDecomposer<V>::maxBoxInDimension(PointIt begin, PointIt end) {
   typename V::CoordType low[V::dimension];
   typename V::CoordType upper[V::dimension];
      
   assert(begin!=end);
   for (unsigned i=0; i<V::dimension; ++i) {
      low[i] = (*begin)->coords.at(i);
      upper[i] = low[i];
   }
   
   for (PointIt it=begin; it!=end; ++it) {
      for (unsigned i=0; i<V::dimension; ++i) {
         low[i] = std::min(low[i], (*it)->coords.at(i));
         upper[i] = std::max(upper[i], (*it)->coords.at(i));
      }
   }
   
   typename V::CoordType maxBox = upper[0]-low[0];
   unsigned maxDim = 0;
   for (unsigned i=1; i<V::dimension; ++i) {
      if ((upper[i]-low[i])>maxBox) {
         maxBox = upper[i]-low[i];
         maxDim = i;
      }
   }
   
   return maxDim;
}
   
// recursive partitioning function
template <typename V>
void RCBDecomposer<V>::buildRCBTree(RCBTree *tree, const int& maxLevel) {
   // compute the longest axis
   // todo: this computation can be done more efficiently
   int partAxis = maxBoxInDimension(tree->_dataBegin, tree->_dataEnd);
   
   // sort along the longest axis
   std::sort(tree->_dataBegin, tree->_dataEnd, VertexComparison<V>(partAxis));
   
   // create two new children and add to the tree
   int size = tree->_dataEnd-tree->_dataBegin;
   if (size<1) return;
   int middle = size/2;
   
   tree->_childL = new RCBTree(tree->_dataBegin, tree->_dataBegin+middle, 
                        tree->_level+1);
   tree->_childR = new RCBTree(tree->_dataBegin+middle, tree->_dataEnd, 
                        tree->_level+1);
   
   // recurse until maximum tree height is reached
   if (tree->_childL->_level<maxLevel) {
      buildRCBTree(tree->_childL, maxLevel);
      buildRCBTree(tree->_childR, maxLevel);
   }
      
   return;
}
    
// recursive leaf collecter
template <typename V>
void RCBDecomposer<V>::collectLeafNodeData(RCBTree *tree, std::vector<PointIt>& begins, 
                              std::vector<PointIt>& ends) {
   if (leaf(tree)) {
      assert(tree != NULL);
      begins.push_back(tree->_dataBegin);
      ends.push_back(tree->_dataEnd);
   } else {
      collectLeafNodeData(tree->_childL, begins, ends);     
      collectLeafNodeData(tree->_childR, begins, ends);     
   }
}

} // namespace DAPTA

---