Actual source code: Topology.hh
petsc-3.4.2 2013-07-02

      1: #ifndef included_ALE_CoSieve_hh
      2: #define included_ALE_CoSieve_hh

      4: #ifndef  included_ALE_Sieve_hh
      5: #include <sieve/Sieve.hh>
      6: #endif

      8: namespace ALE {
      9:   // A Topology is a collection of Sieves
     10:   //   Each Sieve has a label, which we call a \emph{patch}
     11:   //   The collection itself we call a \emph{sheaf}
     12:   //   The main operation we provide in Topology is the creation of a \emph{label}
     13:   //     A label is a bidirectional mapping of Sieve points to integers, implemented with a Sifter
     14:   template<typename Patch_, typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::point_type> >
     15:   class Topology : public ALE::ParallelObject {
     16:   public:
     17:     typedef Patch_                                                patch_type;
     18:     typedef Sieve_                                                sieve_type;
     19:     typedef Alloc_                                                alloc_type;
     20:     typedef typename alloc_type::template rebind<int>::other      int_alloc_type;
     21:     typedef typename sieve_type::point_type                       point_type;
     22:     typedef typename alloc_type::template rebind<std::pair<const patch_type, Obj<sieve_type> > >::other       pair1_alloc_type;
     23:     typedef typename std::map<patch_type, Obj<sieve_type>, std::less<patch_type>, pair1_alloc_type>           sheaf_type;
     24:     typedef typename ALE::Sifter<int, point_type, int>                                                        patch_label_type;
     25:     typedef typename alloc_type::template rebind<patch_label_type>::other                                     patch_label_alloc_type;
     26:     typedef typename patch_label_alloc_type::pointer                                                          patch_label_ptr;
     27:     typedef typename alloc_type::template rebind<std::pair<const patch_type, Obj<patch_label_type> > >::other pair2_alloc_type;
     28:     typedef typename std::map<patch_type, Obj<patch_label_type>, std::less<patch_type>, pair2_alloc_type>     label_type;
     29:     typedef typename alloc_type::template rebind<std::pair<const patch_type, int> >::other                    pair3_alloc_type;
     30:     typedef typename std::map<patch_type, int, std::less<patch_type>, pair3_alloc_type>                       max_label_type;
     31:     typedef typename alloc_type::template rebind<std::pair<const std::string, label_type> >::other            pair4_alloc_type;
     32:     typedef typename std::map<const std::string, label_type, std::less<const std::string>, pair4_alloc_type>  labels_type;
     33:     typedef typename patch_label_type::supportSequence                       label_sequence;
     34:     typedef typename std::set<point_type, std::less<point_type>, alloc_type> point_set_type;
     35:     typedef typename alloc_type::template rebind<point_set_type>::other      point_set_alloc_type;
     36:     typedef typename point_set_alloc_type::pointer                           point_set_ptr;
     37:     typedef typename ALE::Sifter<int,point_type,point_type>                  send_overlap_type;
     38:     typedef typename alloc_type::template rebind<send_overlap_type>::other   send_overlap_alloc_type;
     39:     typedef typename send_overlap_alloc_type::pointer                        send_overlap_ptr;
     40:     typedef typename ALE::Sifter<point_type,int,point_type>                  recv_overlap_type;
     41:     typedef typename alloc_type::template rebind<recv_overlap_type>::other   recv_overlap_alloc_type;
     42:     typedef typename recv_overlap_alloc_type::pointer                        recv_overlap_ptr;
     43:   protected:
     44:     sheaf_type     _sheaf;
     45:     labels_type    _labels;
     46:     int            _maxHeight;
     47:     max_label_type _maxHeights;
     48:     int            _maxDepth;
     49:     max_label_type _maxDepths;
     50:     bool           _calculatedOverlap;
     51:     Obj<send_overlap_type> _sendOverlap;
     52:     Obj<recv_overlap_type> _recvOverlap;
     53:     Obj<send_overlap_type> _distSendOverlap;
     54:     Obj<recv_overlap_type> _distRecvOverlap;
     55:     // Work space
     56:     Obj<point_set_type>    _modifiedPoints;
     57:     alloc_type             _allocator;
     58:     int_alloc_type         _int_allocator;
     59:   public:
     60:     Topology(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1), _calculatedOverlap(false) {
     61:       send_overlap_ptr pSendOverlap = send_overlap_alloc_type(this->_allocator).allocate(1);
     62:       send_overlap_alloc_type(this->_allocator).construct(pSendOverlap, send_overlap_type(this->comm(), this->debug()));
     63:       this->_sendOverlap = Obj<send_overlap_type>(pSendOverlap, sizeof(send_overlap_type));
     64:       ///this->_sendOverlap    = new send_overlap_type(this->comm(), this->debug());
     65:       recv_overlap_ptr pRecvOverlap = recv_overlap_alloc_type(this->_allocator).allocate(1);
     66:       recv_overlap_alloc_type(this->_allocator).construct(pRecvOverlap, recv_overlap_type(this->comm(), this->debug()));
     67:       this->_recvOverlap = Obj<recv_overlap_type>(pRecvOverlap, sizeof(recv_overlap_type));
     68:       ///this->_recvOverlap    = new recv_overlap_type(this->comm(), this->debug());
     69:       point_set_ptr pModPoints = point_set_alloc_type(this->_allocator).allocate(1);
     70:       point_set_alloc_type(this->_allocator).construct(pModPoints, point_set_type());
     71:       this->_modifiedPoints = Obj<point_set_type>(pModPoints, sizeof(point_set_type));
     72:       ///this->_modifiedPoints = new point_set_type();
     73:     };
     74:     virtual ~Topology() {};
     75:   public: // Verifiers
     76:     void checkPatch(const patch_type& patch) {
     77:       if (this->_sheaf.find(patch) == this->_sheaf.end()) {
     78:         ostringstream msg;
     79:         msg << "Invalid topology patch: " << patch << std::endl;
     80:         throw ALE::Exception(msg.str().c_str());
     81:       }
     82:     };
     83:     void checkLabel(const std::string& name, const patch_type& patch) {
     84:       this->checkPatch(patch);
     85:       if ((this->_labels.find(name) == this->_labels.end()) || (this->_labels[name].find(patch) == this->_labels[name].end())) {
     86:         ostringstream msg;
     87:         msg << "Invalid label name: " << name << " for patch " << patch << std::endl;
     88:         throw ALE::Exception(msg.str().c_str());
     89:       }
     90:     };
     91:     bool hasPatch(const patch_type& patch) {
     92:       if (this->_sheaf.find(patch) != this->_sheaf.end()) {
     93:         return true;
     94:       }
     95:       return false;
     96:     };
     97:     bool hasLabel(const std::string& name, const patch_type& patch) {
     98:       if ((this->_labels.find(name) != this->_labels.end()) && (this->_labels[name].find(patch) != this->_labels[name].end())) {
     99:         return true;
    100:       }
    101:       return false;
    102:     };
    103:   public: // Accessors
    104:     const Obj<sieve_type>& getPatch(const patch_type& patch) {
    105:       this->checkPatch(patch);
    106:       return this->_sheaf[patch];
    107:     };
    108:     void setPatch(const patch_type& patch, const Obj<sieve_type>& sieve) {
    109:       this->_sheaf[patch] = sieve;
    110:     };
    111:     int getValue (const Obj<patch_label_type>& label, const point_type& point, const int defValue = 0) {
    112:       const Obj<typename patch_label_type::coneSequence>& cone = label->cone(point);

    114:       if (cone->size() == 0) return defValue;
    115:       return *cone->begin();
    116:     };
    117:     template<typename InputPoints>
    118:     int getMaxValue (const Obj<patch_label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
    119:       int maxValue = defValue;

    121:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
    122:         maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
    123:       }
    124:       return maxValue;
    125:     }
    126:     void setValue(const Obj<patch_label_type>& label, const point_type& point, const int value) {
    127:       label->setCone(value, point);
    128:     };
    129:     const Obj<patch_label_type>& createLabel(const patch_type& patch, const std::string& name) {
    130:       this->checkPatch(patch);
    131:       ///patch_label_ptr pLabel = patch_label_alloc_type(this->_allocator).allocate(1);
    132:       ///patch_label_alloc_type(this->_allocator).construct(pLabel, patch_label_type(this->comm(), this->debug()));
    133:       ///this->_labels[name][patch] = Obj<patch_label_type>(pLabel, sizeof(patch_label_type));
    134:       this->_labels[name][patch] = new patch_label_type(this->comm(), this->debug());
    135:       return this->_labels[name][patch];
    136:     };
    137:     const Obj<patch_label_type>& getLabel(const patch_type& patch, const std::string& name) {
    138:       this->checkLabel(name, patch);
    139:       return this->_labels[name][patch];
    140:     };
    141:     const Obj<label_sequence>& getLabelStratum(const patch_type& patch, const std::string& name, int value) {
    142:       this->checkLabel(name, patch);
    143:       return this->_labels[name][patch]->support(value);
    144:     };
    145:     const sheaf_type& getPatches() {
    146:       return this->_sheaf;
    147:     };
    148:     const labels_type& getLabels() {
    149:       return this->_labels;
    150:     };
    151:     void clear() {
    152:       this->_sheaf.clear();
    153:       this->_labels.clear();
    154:       this->_maxHeight = -1;
    155:       this->_maxHeights.clear();
    156:       this->_maxDepth = -1;
    157:       this->_maxDepths.clear();
    158:     };
    159:     const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
    160:     void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
    161:     const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
    162:     void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
    163:     const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
    164:     void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
    165:     const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
    166:     void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
    167:   public: // Stratification
    168:     template<class InputPoints>
    169:     void computeHeight(const Obj<patch_label_type>& height, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxHeight) {
    170:       this->_modifiedPoints->clear();

    172:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
    173:         // Compute the max height of the points in the support of p, and add 1
    174:         int h0 = this->getValue(height, *p_iter, -1);
    175:         int h1 = this->getMaxValue(height, sieve->support(*p_iter), -1) + 1;

    177:         if(h1 != h0) {
    178:           this->setValue(height, *p_iter, h1);
    179:           if (h1 > maxHeight) maxHeight = h1;
    180:           this->_modifiedPoints->insert(*p_iter);
    181:         }
    182:       }
    183:       // FIX: We would like to avoid the copy here with cone()
    184:       if(this->_modifiedPoints->size() > 0) {
    185:         this->computeHeight(height, sieve, sieve->cone(this->_modifiedPoints), maxHeight);
    186:       }
    187:     }
    188:     void computeHeights() {
    189:       const std::string name("height");

    191:       this->_maxHeight = -1;
    192:       for(typename sheaf_type::iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
    193:         const Obj<patch_label_type>& label = this->createLabel(s_iter->first, name);

    195:         this->_maxHeights[s_iter->first] = -1;
    196:         this->computeHeight(label, s_iter->second, s_iter->second->leaves(), this->_maxHeights[s_iter->first]);
    197:         if (this->_maxHeights[s_iter->first] > this->_maxHeight) this->_maxHeight = this->_maxHeights[s_iter->first];
    198:       }
    199:     }
    200:     int height() const {return this->_maxHeight;};
    201:     int height(const patch_type& patch) {
    202:       this->checkPatch(patch);
    203:       return this->_maxHeights[patch];
    204:     }
    205:     int height(const patch_type& patch, const point_type& point) {
    206:       return this->getValue(this->_labels["height"][patch], point, -1);
    207:     }
    208:     const Obj<label_sequence>& heightStratum(const patch_type& patch, int height) {
    209:       return this->getLabelStratum(patch, "height", height);
    210:     }
    211:     template<class InputPoints>
    212:     void computeDepth(const Obj<patch_label_type>& depth, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxDepth) {
    213:       this->_modifiedPoints->clear();

    215:       for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
    216:         // Compute the max depth of the points in the cone of p, and add 1
    217:         int d0 = this->getValue(depth, *p_iter, -1);
    218:         int d1 = this->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;

    220:         if(d1 != d0) {
    221:           this->setValue(depth, *p_iter, d1);
    222:           if (d1 > maxDepth) maxDepth = d1;
    223:           this->_modifiedPoints->insert(*p_iter);
    224:         }
    225:       }
    226:       // FIX: We would like to avoid the copy here with support()
    227:       if(this->_modifiedPoints->size() > 0) {
    228:         this->computeDepth(depth, sieve, sieve->support(this->_modifiedPoints), maxDepth);
    229:       }
    230:     }
    231:     void computeDepths() {
    232:       const std::string name("depth");

    234:       this->_maxDepth = -1;
    235:       for(typename sheaf_type::iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
    236:         const Obj<patch_label_type>& label = this->createLabel(s_iter->first, name);

    238:         this->_maxDepths[s_iter->first] = -1;
    239:         this->computeDepth(label, s_iter->second, s_iter->second->roots(), this->_maxDepths[s_iter->first]);
    240:         if (this->_maxDepths[s_iter->first] > this->_maxDepth) this->_maxDepth = this->_maxDepths[s_iter->first];
    241:       }
    242:     }
    243:     int depth() const {return this->_maxDepth;};
    244:     int depth(const patch_type& patch) {
    245:       this->checkPatch(patch);
    246:       return this->_maxDepths[patch];
    247:     }
    248:     int depth(const patch_type& patch, const point_type& point) {
    249:       return this->getValue(this->_labels["depth"][patch], point, -1);
    250:     }
    251:     const Obj<label_sequence>& depthStratum(const patch_type& patch, int depth) {
    252:       return this->getLabelStratum(patch, "depth", depth);
    253:     }
    256:     void stratify() {
    257:       ALE_LOG_EVENT_BEGIN;
    258:       this->computeHeights();
    259:       this->computeDepths();
    260:       ALE_LOG_EVENT_END;
    261:     }
    262:   public: // Viewers
    263:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
    264:       if (comm == MPI_COMM_NULL) {
    265:         comm = this->comm();
    266:       }
    267:       if (name == "") {
    268:         PetscPrintf(comm, "viewing a Topology\n");
    269:       } else {
    270:         PetscPrintf(comm, "viewing Topology '%s'\n", name.c_str());
    271:       }
    272:       PetscPrintf(comm, "  maximum height %d maximum depth %d\n", this->height(), this->depth());
    273:       for(typename sheaf_type::const_iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
    274:         ostringstream txt;

    276:         txt << "Patch " << s_iter->first;
    277:         s_iter->second->view(txt.str().c_str(), comm);
    278:         PetscPrintf(comm, "  maximum height %d maximum depth %d\n", this->height(s_iter->first), this->depth(s_iter->first));
    279:       }
    280:       for(typename labels_type::const_iterator l_iter = this->_labels.begin(); l_iter != this->_labels.end(); ++l_iter) {
    281:         PetscPrintf(comm, "  label %s constructed\n", l_iter->first.c_str());
    282:       }
    283:     }
    284:   public:
    285:     void constructOverlap(const patch_type& patch) {
    286:       if (this->_calculatedOverlap) return;
    287:       if (this->hasPatch(patch)) {
    288:         this->constructOverlap(this->getPatch(patch)->base(), this->_sendOverlap, this->_recvOverlap);
    289:         this->constructOverlap(this->getPatch(patch)->cap(), this->_sendOverlap, this->_recvOverlap);
    290:       }
    291:       if (this->debug()) {
    292:         this->_sendOverlap->view("Send overlap");
    293:         this->_recvOverlap->view("Receive overlap");
    294:       }
    295:       this->_calculatedOverlap = true;
    296:     }
    297:     template<typename Sequence>
    298:     void constructOverlap(const Obj<Sequence>& points, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
    299:       point_type *sendBuf = this->_allocator.allocate(points->size());
    300:       for(unsigned int i = 0; i < points->size(); ++i) {this->_allocator.construct(sendBuf+i, point_type());}
    301:       ///point_type *sendBuf = new point_type[points->size()];
    302:       int         size    = 0;
    303:       for(typename Sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
    304:         sendBuf[size++] = *l_iter;
    305:       }
    306:       int *sizes   = this->_int_allocator.allocate(this->commSize()+1);   // The number of points coming from each process
    307:       int *offsets = this->_int_allocator.allocate(this->commSize()+1); // Prefix sums for sizes
    308:       int *oldOffs = this->_int_allocator.allocate(this->commSize()+1); // Temporary storage
    309:       for(int i = 0; i < this->commSize()+1; ++i) {
    310:         this->_int_allocator.construct(sizes+i,   0);
    311:         this->_int_allocator.construct(offsets+i, 0);
    312:         this->_int_allocator.construct(oldOffs+i, 0);
    313:       }
    314:       ///int *sizes   = new int[this->commSize()];
    315:       ///int *offsets = new int[this->commSize()+1];
    316:       ///int *oldOffs = new int[this->commSize()+1];
    317:       point_type *remotePoints = NULL;            // The points from each process
    318:       int        *remoteRanks  = NULL;            // The rank and number of overlap points of each process that overlaps another

    320:       // Change to Allgather() for the correct binning algorithm
    321:       MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, this->comm());
    322:       if (this->commRank() == 0) {
    323:         offsets[0] = 0;
    324:         for(int p = 1; p <= this->commSize(); p++) {
    325:           offsets[p] = offsets[p-1] + sizes[p-1];
    326:         }
    327:         remotePoints = this->_allocator.allocate(offsets[this->commSize()]);
    328:         for(int i = 0; i < offsets[this->commSize()]; ++i) {this->_allocator.construct(remotePoints+i, point_type());}
    329:         ///remotePoints = new point_type[offsets[this->commSize()]];
    330:       }
    331:       MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, this->comm());
    332:       std::map<int, std::map<int, std::set<point_type> > > overlapInfo; // Maps (p,q) to their set of overlap points

    334:       if (this->commRank() == 0) {
    335:         for(int p = 0; p < this->commSize(); p++) {
    336:           std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]]);
    337:         }
    338:         for(int p = 0; p <= this->commSize(); p++) {
    339:           oldOffs[p] = offsets[p];
    340:         }
    341:         for(int p = 0; p < this->commSize(); p++) {
    342:           for(int q = p+1; q < this->commSize(); q++) {
    343:             std::set_intersection(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
    344:                                   &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
    345:                                   std::insert_iterator<std::set<point_type> >(overlapInfo[p][q], overlapInfo[p][q].begin()));
    346:             overlapInfo[q][p] = overlapInfo[p][q];
    347:           }
    348:           sizes[p]     = overlapInfo[p].size()*2;
    349:           offsets[p+1] = offsets[p] + sizes[p];
    350:         }
    351:         remoteRanks = this->_int_allocator.allocate(offsets[this->commSize()]);
    352:         for(int i = 0; i < offsets[this->commSize()]; ++i) {this->_int_allocator.construct(remoteRanks+i, 0);}
    353:         ///remoteRanks = new int[offsets[this->commSize()]];
    354:         int       k = 0;
    355:         for(int p = 0; p < this->commSize(); p++) {
    356:           for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
    357:             remoteRanks[k*2]   = r_iter->first;
    358:             remoteRanks[k*2+1] = r_iter->second.size();
    359:             k++;
    360:           }
    361:         }
    362:       }
    363:       int numOverlaps;                          // The number of processes overlapping this process
    364:       MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, this->comm());
    365:       int *overlapRanks = this->_int_allocator.allocate(numOverlaps);
    366:       for(int i = 0; i < numOverlaps; ++i) {this->_int_allocator.construct(overlapRanks+i, 0);}
    367:       ///int *overlapRanks = new int[numOverlaps]; // The rank and overlap size for each overlapping process
    368:       MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, this->comm());
    369:       point_type *sendPoints = NULL;            // The points to send to each process
    370:       if (this->commRank() == 0) {
    371:         for(int p = 0, k = 0; p < this->commSize(); p++) {
    372:           sizes[p] = 0;
    373:           for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
    374:             sizes[p] += remoteRanks[k*2+1];
    375:             k++;
    376:           }
    377:           offsets[p+1] = offsets[p] + sizes[p];
    378:         }
    379:         sendPoints = this->_allocator.allocate(offsets[this->commSize()]);
    380:         for(int i = 0; i < offsets[this->commSize()]; ++i) {this->_allocator.construct(sendPoints+i, point_type());}
    381:         ///sendPoints = new point_type[offsets[this->commSize()]];
    382:         for(int p = 0, k = 0; p < this->commSize(); p++) {
    383:           for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
    384:             int rank = r_iter->first;
    385:             for(typename std::set<point_type>::iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
    386:               sendPoints[k++] = *p_iter;
    387:             }
    388:           }
    389:         }
    390:       }
    391:       int numOverlapPoints = 0;
    392:       for(int r = 0; r < numOverlaps/2; r++) {
    393:         numOverlapPoints += overlapRanks[r*2+1];
    394:       }
    395:       point_type *overlapPoints = this->_allocator.allocate(numOverlapPoints);
    396:       for(int i = 0; i < numOverlapPoints; ++i) {this->_allocator.construct(overlapPoints+i, point_type());}
    397:       ///point_type *overlapPoints = new point_type[numOverlapPoints];
    398:       MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints, MPI_INT, 0, this->comm());

    400:       for(int r = 0, k = 0; r < numOverlaps/2; r++) {
    401:         int rank = overlapRanks[r*2];

    403:         for(int p = 0; p < overlapRanks[r*2+1]; p++) {
    404:           point_type point = overlapPoints[k++];

    406:           sendOverlap->addArrow(point, rank, point);
    407:           recvOverlap->addArrow(rank, point, point);
    408:         }
    409:       }

    411:       delete [] overlapPoints;
    412:       delete [] overlapRanks;
    413:       delete [] sizes;
    414:       delete [] offsets;
    415:       delete [] oldOffs;
    416:       if (this->commRank() == 0) {
    417:         delete [] remoteRanks;
    418:         delete [] remotePoints;
    419:         delete [] sendPoints;
    420:       }
    421:     }
    422:   };

    424:   // An Overlap is a Sifter describing the overlap of two Sieves
    425:   //   Each arrow is local point ---(remote point)---> remote rank right now
    426:   //     For XSifter, this should change to (local patch, local point) ---> (remote rank, remote patch, remote point)
    427: }

    429: #endif