Actual source code: ParallelMapping.hh
1: #ifndef included_ALE_ParallelMapping_hh
2: #define included_ALE_ParallelMapping_hh
4: #ifndef included_ALE_IField_hh
5: #include <IField.hh>
6: #endif
8: #ifndef included_ALE_Sections_hh
9: #include <Sections.hh>
10: #endif
12: #include <functional>
16: namespace ALE {
17: template<class _Tp>
18: struct Identity : public std::unary_function<_Tp,_Tp>
19: {
20: _Tp& operator()(_Tp& x) const {return x;}
21: const _Tp& operator()(const _Tp& x) const {return x;}
22: };
24: template<class _Tp>
25: struct IsEqual : public std::unary_function<_Tp, bool>, public std::binary_function<_Tp, _Tp, bool>
26: {
27: const _Tp& x;
28: IsEqual(const _Tp& x) : x(x) {};
29: bool operator()(_Tp& y) const {return x == y;}
30: const bool operator()(const _Tp& y) const {return x == y;}
31: bool operator()(_Tp& y, _Tp& dummy) const {return x == y;}
32: const bool operator()(const _Tp& y, const _Tp& dummy) const {return x == y;}
33: };
35: // Creates new global point names and renames local points globally
36: template<typename Point>
37: class PointFactory : ALE::ParallelObject {
38: public:
39: typedef Point point_type;
40: typedef std::map<point_type,point_type> renumbering_type;
41: typedef std::map<int,std::map<point_type,point_type> > remote_renumbering_type;
42: protected:
43: point_type originalMax;
44: point_type currentMax;
45: renumbering_type renumbering;
46: renumbering_type invRenumbering;
47: remote_renumbering_type remoteRenumbering;
48: protected:
49: PointFactory(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), originalMax(-1) {};
50: public:
51: ~PointFactory() {};
52: public:
53: static PointFactory& singleton(MPI_Comm comm, const point_type& maxPoint, const int debug = 0, bool cleanup = false) {
54: static PointFactory *_singleton = NULL;
56: if (cleanup) {
57: if (debug) {std::cout << "Destroying PointFactory" << std::endl;}
58: if (_singleton) {delete _singleton;}
59: _singleton = NULL;
60: } else if (_singleton == NULL) {
61: if (debug) {std::cout << "Creating new PointFactory" << std::endl;}
62: _singleton = new PointFactory(comm, debug);
63: _singleton->setMax(maxPoint);
64: }
65: return *_singleton;
66: };
67: void setMax(const point_type& maxPoint) {
68: PetscErrorCode MPI_Allreduce((void *) &maxPoint, &this->originalMax, 1, MPI_INT, MPI_MAX, this->comm());CHKERRXX(ierr);
69: ++this->originalMax;
70: this->currentMax = this->originalMax;
71: };
72: void clear() {
73: this->currentMax = this->originalMax;
74: this->renumbering.clear();
75: this->invRenumbering.clear();
76: };
77: public:
78: template<typename Iterator>
79: void renumberPoints(const Iterator& begin, const Iterator& end) {
80: renumberPoints(begin, end, Identity<typename Iterator::value_type>());
81: };
82: template<typename Iterator, typename KeyExtractor>
83: void renumberPoints(const Iterator& begin, const Iterator& end, const KeyExtractor& ex) {
84: int numPoints = 0, numGlobalPoints, firstPoint;
86: for(Iterator p_iter = begin; p_iter != end; ++p_iter) ++numPoints;
87: MPI_Allreduce(&numPoints, &numGlobalPoints, 1, MPI_INT, MPI_SUM, this->comm());
88: MPI_Scan(&numPoints, &firstPoint, 1, MPI_INT, MPI_SUM, this->comm());
89: firstPoint += this->currentMax - numPoints;
90: this->currentMax += numGlobalPoints;
91: for(Iterator p_iter = begin; p_iter != end; ++p_iter, ++firstPoint) {
92: if (this->debug()) {std::cout << "["<<this->commRank()<<"]: New point " << ex(*p_iter) << " --> " << firstPoint << std::endl;}
93: this->renumbering[firstPoint] = ex(*p_iter);
94: this->invRenumbering[ex(*p_iter)] = firstPoint;
95: }
96: };
97: public:
98: void constructRemoteRenumbering() {
99: const int localSize = this->renumbering.size();
100: int *remoteSizes = new int[this->commSize()];
101: int *localMap = new int[localSize*2];
102: int *recvCounts = new int[this->commSize()];
103: int *displs = new int[this->commSize()];
105: // Populate arrays
106: int r = 0;
107: for(typename renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter, ++r) {
108: localMap[r*2+0] = r_iter->first;
109: localMap[r*2+1] = r_iter->second;
110: }
111: // Communicate renumbering sizes
112: MPI_Allgather((void*) &localSize, 1, MPI_INT, remoteSizes, 1, MPI_INT, this->comm());
113: for(int p = 0; p < this->commSize(); ++p) {
114: recvCounts[p] = remoteSizes[p]*2;
115: if (p == 0) {
116: displs[p] = 0;
117: } else {
118: displs[p] = displs[p-1] + recvCounts[p-1];
119: }
120: }
121: int *remoteMaps = new int[displs[this->commSize()-1]+recvCounts[this->commSize()-1]];
122: // Communicate renumberings
123: MPI_Allgatherv(localMap, localSize*2, MPI_INT, remoteMaps, recvCounts, displs, MPI_INT, this->comm());
124: // Populate maps
125: for(int p = 0; p < this->commSize(); ++p) {
126: if (p == this->commRank()) continue;
127: int offset = displs[p];
129: for(int r = 0; r < remoteSizes[p]; ++r) {
130: this->remoteRenumbering[p][remoteMaps[r*2+0+offset]] = remoteMaps[r*2+1+offset];
131: if (this->debug()) {std::cout << "["<<this->commRank()<<"]: Remote renumbering["<<p<<"] " << remoteMaps[r*2+0+offset] << " --> " << remoteMaps[r*2+1+offset] << std::endl;}
132: }
133: }
134: // Cleanup
135: delete [] recvCounts;
136: delete [] displs;
137: delete [] localMap;
138: delete [] remoteMaps;
139: delete [] remoteSizes;
140: };
141: public:
142: // global point --> local point
143: renumbering_type& getRenumbering() {
144: return this->renumbering;
145: };
146: // local point --> global point
147: renumbering_type& getInvRenumbering() {
148: return this->invRenumbering;
149: };
150: // rank --> global point --> local point
151: remote_renumbering_type& getRemoteRenumbering() {
152: return this->remoteRenumbering;
153: };
154: };
156: // TODO: Check MPI return values and status of Waits
157: template<typename Value_>
158: class MPIMover : public ALE::ParallelObject {
159: public:
160: typedef Value_ value_type;
161: typedef size_t num_type;
162: typedef std::pair<num_type, const value_type*> move_type;
163: typedef std::map<int, move_type> moves_type;
164: typedef std::vector<MPI_Request> requests_type;
165: protected:
166: bool _createdType;
167: int _tag;
168: MPI_Datatype _datatype;
169: moves_type _sends;
170: moves_type _recvs;
171: requests_type _requests;
172: public:
173: MPIMover(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _createdType(0) {
174: this->_tag = this->getNewTag();
175: this->_datatype = this->getMPIDatatype();
176: };
177: MPIMover(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _createdType(0), _tag(tag) {
178: this->_datatype = this->getMPIDatatype();
179: };
180: MPIMover(MPI_Comm comm, const MPI_Datatype datatype, const int tag, const int debug) : ParallelObject(comm, debug), _createdType(0) {
181: if (tag == MPI_UNDEFINED) {
182: this->_tag = this->getNewTag();
183: } else {
184: this->_tag = tag;
185: }
186: if (datatype == MPI_DATATYPE_NULL) {
187: this->_datatype = this->getMPIDatatype();
188: } else {
189: this->_datatype = datatype;
190: }
191: };
192: ~MPIMover() {
193: if (_createdType) {
194: MPI_Type_free(&this->_datatype);
195: }
196: };
197: protected:
198: // TODO: Can rewrite this with template specialization?
199: MPI_Datatype getMPIDatatype() {
200: if (sizeof(value_type) == 1) {
201: return MPI_BYTE;
202: } else if (sizeof(value_type) == 2) {
203: return MPI_SHORT;
204: } else if (sizeof(value_type) == 4) {
205: return MPI_INT;
206: } else if (sizeof(value_type) == 8) {
207: return MPI_DOUBLE;
208: } else if (sizeof(value_type) == 28) {
209: int blen[2];
210: MPI_Aint indices[2];
211: MPI_Datatype oldtypes[2], newtype;
212: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
213: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
214: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
215: MPI_Type_commit(&newtype);
216: this->_createdType = true;
217: return newtype;
218: } else if (sizeof(value_type) == 32) {
219: int blen[2];
220: MPI_Aint indices[2];
221: MPI_Datatype oldtypes[2], newtype;
222: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_DOUBLE;
223: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
224: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
225: MPI_Type_commit(&newtype);
226: this->_createdType = true;
227: return newtype;
228: }
229: throw ALE::Exception("Cannot determine MPI type for value type");
230: };
231: int getNewTag() const {
232: static int tagKeyval = MPI_KEYVAL_INVALID;
233: int *tagvalp = NULL, *maxval, flg;
235: if (tagKeyval == MPI_KEYVAL_INVALID) {
236: tagvalp = (int *) malloc(sizeof(int));
237: MPI_Keyval_create(MPI_NULL_COPY_FN, Mesh_DelTag, &tagKeyval, (void *) NULL);
238: MPI_Attr_put(this->comm(), tagKeyval, tagvalp);
239: tagvalp[0] = 0;
240: }
241: MPI_Attr_get(this->comm(), tagKeyval, (void **) &tagvalp, &flg);
242: if (tagvalp[0] < 1) {
243: MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
244: tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
245: }
246: if (this->debug()) {
247: std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
248: }
249: return tagvalp[0]--;
250: };
251: void constructRequests() {
252: this->_requests.clear();
254: for(typename moves_type::const_iterator s_iter = this->_sends.begin(); s_iter != this->_sends.end(); ++s_iter) {
255: const int rank = s_iter->first;
256: const int num = s_iter->second.first;
257: void *data = (void *) s_iter->second.second;
258: MPI_Request request;
260: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << num << ") to " << rank << " tag " << this->_tag << std::endl;}
261: MPI_Send_init(data, num, this->_datatype, rank, this->_tag, this->comm(), &request);
262: this->_requests.push_back(request);
263: #if defined(PETSC_USE_LOG)
264: // PETSc logging
265: isend_ct++;
266: TypeSize(&isend_len, num, this->_datatype);
267: #endif
268: }
269: for(typename moves_type::const_iterator r_iter = this->_recvs.begin(); r_iter != this->_recvs.end(); ++r_iter) {
270: const int rank = r_iter->first;
271: const int num = r_iter->second.first;
272: void *data = (void *) r_iter->second.second;
273: MPI_Request request;
275: if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data (" << num << ") from " << rank << " tag " << this->_tag << std::endl;}
276: MPI_Recv_init(data, num, this->_datatype, rank, this->_tag, this->comm(), &request);
277: this->_requests.push_back(request);
278: #if defined(PETSC_USE_LOG)
279: // PETSc logging
280: irecv_ct++;
281: TypeSize(&irecv_len, num, this->_datatype);
282: #endif
283: }
284: };
285: public:
286: void send(const int rank, const int num, const value_type *data) {
287: this->_sends[rank] = move_type(num, data);
288: };
289: void recv(const int rank, const int num, const value_type *data) {
290: this->_recvs[rank] = move_type(num, data);
291: };
292: void start() {
293: this->constructRequests();
294: for(typename requests_type::const_iterator r_iter = this->_requests.begin(); r_iter != this->_requests.end(); ++r_iter) {
295: MPI_Request request = *r_iter;
297: MPI_Start(&request);
298: }
299: };
300: void end() {
301: MPI_Status status;
303: for(typename requests_type::const_iterator r_iter = this->_requests.begin(); r_iter != this->_requests.end(); ++r_iter) {
304: MPI_Request request = *r_iter;
306: MPI_Wait(&request, &status);
307: }
308: };
309: };
310: template<typename Alloc_ = malloc_allocator<int> >
311: class OverlapBuilder {
312: public:
313: typedef Alloc_ alloc_type;
314: protected:
315: template<typename T>
316: struct lessPair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
317: bool operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
318: return x.first < y.first;
319: }
320: };
321: template<typename T>
322: struct mergePair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
323: std::pair<T,std::pair<T,T> > operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
324: return std::pair<T,std::pair<T,T> >(x.first, std::pair<T,T>(x.second, y.second));
325: }
326: };
327: template<typename _InputIterator1, typename _InputIterator2, typename _OutputIterator, typename _Compare, typename _Merge>
328: static _OutputIterator set_intersection_merge(_InputIterator1 __first1, _InputIterator1 __last1,
329: _InputIterator2 __first2, _InputIterator2 __last2,
330: _OutputIterator __result, _Compare __comp, _Merge __merge)
331: {
332: while(__first1 != __last1 && __first2 != __last2) {
333: if (__comp(*__first1, *__first2))
334: ++__first1;
335: else if (__comp(*__first2, *__first1))
336: ++__first2;
337: else
338: {
339: *__result = __merge(*__first1, *__first2);
340: ++__first1;
341: ++__first2;
342: ++__result;
343: }
344: }
345: return __result;
346: };
347: public:
348: template<typename Sequence, typename Renumbering, typename SendOverlap, typename RecvOverlap>
349: static void constructOverlap(const Sequence& points, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
350: typedef typename SendOverlap::source_type point_type;
351: typedef std::pair<point_type,point_type> pointPair;
352: typedef std::pair<point_type,pointPair> pointTriple;
353: alloc_type allocator;
354: typename alloc_type::template rebind<point_type>::other point_allocator;
355: typename alloc_type::template rebind<pointPair>::other pointPair_allocator;
356: const MPI_Comm comm = sendOverlap->comm();
357: const int commSize = sendOverlap->commSize();
358: const int commRank = sendOverlap->commRank();
359: point_type *sendBuf = point_allocator.allocate(points.size()*2);
360: for(size_t i = 0; i < points.size()*2; ++i) {point_allocator.construct(sendBuf+i, point_type());}
361: int size = 0;
362: const int debug = sendOverlap->debug();
363: for(typename Sequence::const_iterator l_iter = points.begin(); l_iter != points.end(); ++l_iter) {
364: if (debug) {std::cout << "["<<commRank<<"]Send point["<<size<<"]: " << *l_iter << " " << renumbering[*l_iter] << std::endl;}
365: sendBuf[size++] = *l_iter;
366: sendBuf[size++] = renumbering[*l_iter];
367: }
368: int *sizes = allocator.allocate(commSize*3+2); // [size] The number of points coming from each process
369: for(int i = 0; i < commSize*3+2; ++i) {allocator.construct(sizes+i, 0);}
370: int *offsets = sizes+commSize; // [size+1] Prefix sums for sizes
371: int *oldOffs = offsets+commSize+1; // [size+1] Temporary storage
372: pointPair *remotePoints = NULL; // The points from each process
373: int *remoteRanks = NULL; // The rank and number of overlap points of each process that overlaps another
374: int numRemotePoints = 0;
375: int numRemoteRanks = 0;
377: // Change to Allgather() for the correct binning algorithm
378: MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, comm);
379: if (commRank == 0) {
380: offsets[0] = 0;
381: for(int p = 1; p <= commSize; p++) {
382: offsets[p] = offsets[p-1] + sizes[p-1];
383: }
384: numRemotePoints = offsets[commSize];
385: remotePoints = pointPair_allocator.allocate(numRemotePoints/2);
386: for(int i = 0; i < numRemotePoints/2; ++i) {pointPair_allocator.construct(remotePoints+i, pointPair());}
387: }
388: MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, comm);
389: for(size_t i = 0; i < points.size(); ++i) {point_allocator.destroy(sendBuf+i);}
390: point_allocator.deallocate(sendBuf, points.size());
391: std::map<int, std::map<int, std::set<pointTriple> > > overlapInfo; // Maps (p,q) to their set of overlap points
393: if (commRank == 0) {
394: for(int p = 0; p <= commSize; p++) {
395: offsets[p] /= 2;
396: }
397: for(int p = 0; p < commSize; p++) {
398: std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]], lessPair<point_type>());
399: }
400: for(int p = 0; p <= commSize; p++) {
401: oldOffs[p] = offsets[p];
402: }
403: for(int p = 0; p < commSize; p++) {
404: for(int q = 0; q < commSize; q++) {
405: if (p == q) continue;
406: set_intersection_merge(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
407: &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
408: std::insert_iterator<std::set<pointTriple> >(overlapInfo[p][q], overlapInfo[p][q].begin()),
409: lessPair<point_type>(), mergePair<point_type>());
410: }
411: sizes[p] = overlapInfo[p].size()*2;
412: offsets[p+1] = offsets[p] + sizes[p];
413: }
414: numRemoteRanks = offsets[commSize];
415: remoteRanks = allocator.allocate(numRemoteRanks);
416: for(int i = 0; i < numRemoteRanks; ++i) {allocator.construct(remoteRanks+i, 0);}
417: int k = 0;
418: for(int p = 0; p < commSize; p++) {
419: for(typename std::map<int, std::set<pointTriple> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
420: remoteRanks[k*2] = r_iter->first;
421: remoteRanks[k*2+1] = r_iter->second.size();
422: k++;
423: }
424: }
425: }
426: int numOverlaps; // The number of processes overlapping this process
427: MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, comm);
428: int *overlapRanks = allocator.allocate(numOverlaps); // The rank and overlap size for each overlapping process
429: for(int i = 0; i < numOverlaps; ++i) {allocator.construct(overlapRanks+i, 0);}
430: MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, comm);
431: point_type *sendPoints = NULL; // The points to send to each process
432: int numSendPoints = 0;
433: if (commRank == 0) {
434: for(int p = 0, k = 0; p < commSize; p++) {
435: sizes[p] = 0;
436: for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
437: sizes[p] += remoteRanks[k*2+1]*2;
438: k++;
439: }
440: offsets[p+1] = offsets[p] + sizes[p];
441: }
442: numSendPoints = offsets[commSize];
443: sendPoints = point_allocator.allocate(numSendPoints*2);
444: for(int i = 0; i < numSendPoints*2; ++i) {point_allocator.construct(sendPoints+i, point_type());}
445: for(int p = 0, k = 0; p < commSize; p++) {
446: for(typename std::map<int, std::set<pointTriple> >::const_iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
447: int rank = r_iter->first;
448: for(typename std::set<pointTriple>::const_iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
449: sendPoints[k++] = p_iter->first;
450: sendPoints[k++] = p_iter->second.second;
451: if (debug) {std::cout << "["<<commRank<<"]Sending points " << p_iter->first << " " << p_iter->second.second << " to rank " << rank << std::endl;}
452: }
453: }
454: }
455: }
456: int numOverlapPoints = 0;
457: for(int r = 0; r < numOverlaps/2; r++) {
458: numOverlapPoints += overlapRanks[r*2+1];
459: }
460: point_type *overlapPoints = point_allocator.allocate(numOverlapPoints*2);
461: for(int i = 0; i < numOverlapPoints*2; ++i) {point_allocator.construct(overlapPoints+i, point_type());}
462: MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints*2, MPI_INT, 0, comm);
464: for(int r = 0, k = 0; r < numOverlaps/2; r++) {
465: int rank = overlapRanks[r*2];
467: for(int p = 0; p < overlapRanks[r*2+1]; p++) {
468: point_type point = overlapPoints[k++];
469: point_type remotePoint = overlapPoints[k++];
471: if (debug) {std::cout << "["<<commRank<<"]Matched up remote point " << remotePoint << "("<<point<<") to local " << renumbering[point] << std::endl;}
472: sendOverlap->addArrow(renumbering[point], rank, remotePoint);
473: recvOverlap->addArrow(rank, renumbering[point], remotePoint);
474: }
475: }
477: for(int i = 0; i < numOverlapPoints; ++i) {point_allocator.destroy(overlapPoints+i);}
478: point_allocator.deallocate(overlapPoints, numOverlapPoints);
479: for(int i = 0; i < numOverlaps; ++i) {allocator.destroy(overlapRanks+i);}
480: allocator.deallocate(overlapRanks, numOverlaps);
481: for(int i = 0; i < commSize*3+2; ++i) {allocator.destroy(sizes+i);}
482: allocator.deallocate(sizes, commSize*3+2);
483: if (commRank == 0) {
484: for(int i = 0; i < numRemoteRanks; ++i) {allocator.destroy(remoteRanks+i);}
485: allocator.deallocate(remoteRanks, numRemoteRanks);
486: for(int i = 0; i < numRemotePoints; ++i) {pointPair_allocator.destroy(remotePoints+i);}
487: pointPair_allocator.deallocate(remotePoints, numRemotePoints);
488: for(int i = 0; i < numSendPoints; ++i) {point_allocator.destroy(sendPoints+i);}
489: point_allocator.deallocate(sendPoints, numSendPoints);
490: }
491: };
492: };
493: namespace Pullback {
494: class SimpleCopy {
495: public:
496: // Copy the overlap section to the related processes
497: // This version is for Constant sections, meaning the same, single value over all points
498: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
499: static void copyConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
500: MPIMover<char> pMover(sendSection->comm(), sendSection->debug());
501: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
502: std::map<int, char *> sendPoints;
503: std::map<int, char *> recvPoints;
504: typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
505: typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;
507: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
508: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
509: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
511: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
512: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
513: const int pSize = points->size();
514: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
515: char *v = sendAllocator.allocate(points->size());
516: int k = 0;
518: for(int i = 0; i < pSize; ++i) {sendAllocator.construct(v+i, 0);}
519: for(typename SendOverlap::coneSequence::iterator p_iter = points->begin(); p_iter != pEnd; ++p_iter, ++k) {
520: v[k] = (char) sendSection->hasPoint(*p_iter);
521: }
522: sendPoints[*r_iter] = v;
523: pMover.send(*r_iter, pSize, sendPoints[*r_iter]);
524: vMover.send(*r_iter, 2, sValues);
525: }
526: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
527: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
528: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
530: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
531: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
532: const int pSize = points->size();
533: char *v = recvAllocator.allocate(pSize);
535: for(int i = 0; i < pSize; ++i) {recvAllocator.construct(v+i, 0);}
536: recvPoints[*r_iter] = v;
537: pMover.recv(*r_iter, pSize, recvPoints[*r_iter]);
538: vMover.recv(*r_iter, 2, rValues);
539: }
540: pMover.start();
541: pMover.end();
542: vMover.start();
543: vMover.end();
544: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
545: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
546: const typename RecvOverlap::traits::supportSequence::iterator pEnd = points->end();
547: const char *v = recvPoints[*r_iter];
548: int k = 0;
550: for(typename RecvOverlap::traits::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter, ++k) {
551: if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(*r_iter, s_iter.color()));}
552: }
553: }
554: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
555: sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->cone(*r_iter)->size());
556: }
557: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
558: recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->support(*r_iter)->size());
559: }
560: };
561: // Specialize to ArrowSection
562: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
563: static void copyConstantArrow(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
564: MPIMover<char> pMover(sendSection->comm(), sendSection->debug());
565: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
566: std::map<int, char *> sendPoints;
567: std::map<int, char *> recvPoints;
568: typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
569: typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;
571: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
572: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
573: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
575: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
576: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
577: const int pSize = points->size();
578: const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
579: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
580: char *v = sendAllocator.allocate(pSize*pSize);
581: int k = 0;
583: for(size_t i = 0; i < pSize*pSize; ++i) {sendAllocator.construct(v+i, 0);}
584: for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
585: for(typename SendOverlap::coneSequence::iterator q_iter = pBegin; q_iter != pEnd; ++q_iter, ++k) {
586: v[k] = (char) sendSection->hasPoint(typename SendSection::point_type(*p_iter,*q_iter));
587: }
588: }
589: sendPoints[*r_iter] = v;
590: pMover.send(*r_iter, pSize*pSize, sendPoints[*r_iter]);
591: vMover.send(*r_iter, 2, sValues);
592: }
593: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
594: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
595: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
597: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
598: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
599: const int pSize = points->size();
600: char *v = recvAllocator.allocate(points->size()*points->size());
602: for(size_t i = 0; i < pSize*pSize; ++i) {recvAllocator.construct(v+i, 0);}
603: recvPoints[*r_iter] = v;
604: pMover.recv(*r_iter, pSize*pSize, recvPoints[*r_iter]);
605: vMover.recv(*r_iter, 2, rValues);
606: }
607: pMover.start();
608: pMover.end();
609: vMover.start();
610: vMover.end();
611: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
612: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
613: const typename RecvOverlap::traits::supportSequence::iterator pBegin = points->begin();
614: const typename RecvOverlap::traits::supportSequence::iterator pEnd = points->end();
615: const char *v = recvPoints[*r_iter];
616: int k = 0;
618: for(typename RecvOverlap::traits::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
619: for(typename RecvOverlap::traits::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter, ++k) {
620: if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(s_iter.color(),t_iter.color()));}
621: }
622: }
623: }
624: };
625: // Copy the overlap section to the related processes
626: // This version is for IConstant sections, meaning the same, single value over all points
627: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
628: static void copyIConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
629: MPIMover<typename SendSection::point_type> pMover(sendSection->comm(), sendSection->debug());
630: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
631: std::map<int, typename SendSection::point_type *> sendPoints;
632: std::map<int, typename SendSection::point_type *> recvPoints;
633: typename SendSection::alloc_type::template rebind<typename SendSection::point_type>::other sendAllocator;
634: typename RecvSection::alloc_type::template rebind<typename SendSection::point_type>::other recvAllocator;
636: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
637: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
638: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
640: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
641: typename SendSection::point_type *v = sendAllocator.allocate(2);
643: for(size_t i = 0; i < 2; ++i) {sendAllocator.construct(v+i, 0);}
644: v[0] = sendSection->getChart().min();
645: v[1] = sendSection->getChart().max();
646: sendPoints[*r_iter] = v;
647: pMover.send(*r_iter, 2, sendPoints[*r_iter]);
648: vMover.send(*r_iter, 2, sValues);
649: ///std::cout << "["<<sendOverlap->commRank()<<"]Sending chart (" << v[0] << ", " << v[1] << ") with values (" << sValues[0] << ", " << sValues[1] << ") to process " << *r_iter << std::endl;
650: }
651: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
652: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
653: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
655: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
656: typename SendSection::point_type *v = recvAllocator.allocate(2);
658: for(size_t i = 0; i < 2; ++i) {recvAllocator.construct(v+i, 0);}
659: recvPoints[*r_iter] = v;
660: pMover.recv(*r_iter, 2, recvPoints[*r_iter]);
661: vMover.recv(*r_iter, 2, rValues);
662: }
663: pMover.start();
664: pMover.end();
665: vMover.start();
666: vMover.end();
668: typename SendSection::point_type min = -1;
669: typename SendSection::point_type max = -1;
671: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
672: const typename RecvSection::point_type *v = recvPoints[*r_iter];
673: typename SendSection::point_type newMin = v[0];
674: typename SendSection::point_type newMax = v[1]-1;
675: //int pSize = 0;
677: ///std::cout << "["<<recvOverlap->commRank()<<"]Received chart (" << v[0] << ", " << v[1] << ") from process " << *r_iter << std::endl;
678: #if 0
679: // Translate to local numbering
680: if (recvOverlap->support(*r_iter)->size()) {
681: while(!pSize) {
682: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMin);
683: pSize = points->size();
684: if (pSize) {
685: newMin = *points->begin();
686: } else {
687: newMin++;
688: }
689: }
690: pSize = 0;
691: while(!pSize) {
692: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMax);
693: pSize = points->size();
694: if (pSize) {
695: newMax = *points->begin();
696: } else {
697: newMax--;
698: }
699: }
700: }
701: std::cout << "["<<recvOverlap->commRank()<<"]Translated to chart (" << newMin << ", " << newMax+1 << ") from process " << *r_iter << std::endl;
702: #endif
703: // Update chart
704: if (min < 0) {
705: min = newMin;
706: max = newMax+1;
707: } else {
708: min = std::min(min, newMin);
709: max = std::max(max, (typename SendSection::point_type) (newMax+1));
710: }
711: }
712: if (!rRanks->size()) {min = max = 0;}
713: recvSection->setChart(typename RecvSection::chart_type(min, max));
714: // TODO: deallocation
715: };
716: // Copy the overlap section to the related processes
717: // This version is for different sections, possibly with different data types
718: // TODO: Can cache MPIMover objects (like a VecScatter)
719: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
720: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection, const MPI_Datatype datatype = MPI_DATATYPE_NULL) {
721: typedef std::pair<int, typename SendSection::value_type *> allocPair;
722: typedef typename RecvSection::point_type recv_point_type;
723: const Obj<typename SendSection::atlas_type>& sendAtlas = sendSection->getAtlas();
724: const Obj<typename RecvSection::atlas_type>& recvAtlas = recvSection->getAtlas();
725: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
726: std::map<int, allocPair> sendValues;
727: std::map<int, allocPair> recvValues;
728: typename SendSection::alloc_type sendAllocator;
729: typename RecvSection::alloc_type recvAllocator;
731: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
732: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
733: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
735: // TODO: This should be const_iterator, but Sifter sucks
736: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
737: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
738: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
739: int numVals = 0;
741: // TODO: This should be const_iterator, but Sifter sucks
742: for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
743: numVals += sendSection->getFiberDimension(*c_iter);
744: }
745: typename SendSection::value_type *v = sendAllocator.allocate(numVals);
746: int k = 0;
748: for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
749: for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
750: const typename SendSection::value_type *vals = sendSection->restrictPoint(*c_iter);
752: for(int i = 0; i < sendSection->getFiberDimension(*c_iter); ++i, ++k) v[k] = vals[i];
753: }
754: sendValues[*r_iter] = allocPair(numVals, v);
755: vMover.send(*r_iter, numVals, v);
756: }
757: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
758: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
760: recvSection->allocatePoint();
761: // TODO: This should be const_iterator, but Sifter sucks
762: int maxVals = 0;
763: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
764: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
765: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
766: int numVals = 0;
768: // TODO: This should be const_iterator, but Sifter sucks
769: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
770: numVals += recvSection->getFiberDimension(recv_point_type(*r_iter, s_iter.color()));
771: }
772: typename SendSection::value_type *v = sendAllocator.allocate(numVals);
774: for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
775: recvValues[*r_iter] = allocPair(numVals, v);
776: vMover.recv(*r_iter, numVals, v);
777: maxVals = std::max(maxVals, numVals);
778: }
779: vMover.start();
780: vMover.end();
781: typename RecvSection::value_type *convertedValues = recvAllocator.allocate(maxVals);
782: for(int i = 0; i < maxVals; ++i) {recvAllocator.construct(convertedValues+i, 0);}
783: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
784: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
785: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
786: typename SendSection::value_type *v = recvValues[*r_iter].second;
787: const int numVals = recvValues[*r_iter].first;
788: int k = 0;
790: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
791: const int size = recvSection->getFiberDimension(recv_point_type(*r_iter, s_iter.color()));
793: for(int i = 0; i < size; ++i) {convertedValues[i] = (typename RecvSection::value_type) v[k+i];}
794: if (size) {recvSection->updatePoint(recv_point_type(*r_iter, s_iter.color()), convertedValues);}
795: k += size;
796: }
797: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
798: }
799: for(int i = 0; i < maxVals; ++i) {recvAllocator.destroy(convertedValues+i);}
800: recvAllocator.deallocate(convertedValues, maxVals);
801: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
802: typename SendSection::value_type *v = sendValues[*r_iter].second;
803: const int numVals = sendValues[*r_iter].first;
805: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
806: sendAllocator.deallocate(v, numVals);
807: }
808: //recvSection->view("After copy");
809: };
810: // Copy the overlap section to the related processes
811: // This version is for sections with the same type
812: template<typename SendOverlap, typename RecvOverlap, typename Section>
813: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& sendSection, const Obj<Section>& recvSection, const MPI_Datatype datatype = MPI_DATATYPE_NULL) {
814: typedef std::pair<int, typename Section::value_type *> allocPair;
815: const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
816: const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
817: MPIMover<typename Section::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
818: std::map<int, allocPair> sendValues;
819: std::map<int, allocPair> recvValues;
820: typename Section::alloc_type allocator;
822: ///sendAtlas->view("Send Atlas in same type copy()");
823: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
824: ///recvAtlas->view("Recv Atlas after same type copy()");
825: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
826: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
828: // TODO: This should be const_iterator, but Sifter sucks
829: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
830: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
831: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
832: int numVals = 0;
834: // TODO: This should be const_iterator, but Sifter sucks
835: for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
836: numVals += sendSection->getFiberDimension(*c_iter);
837: }
838: typename Section::value_type *v = allocator.allocate(numVals);
839: int k = 0;
841: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
842: for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
843: const typename Section::value_type *vals = sendSection->restrictPoint(*c_iter);
844: const int dim = sendSection->getFiberDimension(*c_iter);
846: for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
847: }
848: sendValues[*r_iter] = allocPair(numVals, v);
849: vMover.send(*r_iter, numVals, v);
850: }
851: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
852: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
854: recvSection->allocatePoint();
855: ///recvSection->view("Recv Section after same type copy() allocation");
856: // TODO: This should be const_iterator, but Sifter sucks
857: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
858: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
859: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
860: int numVals = 0;
862: // TODO: This should be const_iterator, but Sifter sucks
863: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
864: numVals += recvSection->getFiberDimension(s_iter.color());
865: }
866: typename Section::value_type *v = allocator.allocate(numVals);
868: recvValues[*r_iter] = allocPair(numVals, v);
869: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
870: vMover.recv(*r_iter, numVals, v);
871: }
872: vMover.start();
873: vMover.end();
874: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
875: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
876: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
877: typename Section::value_type *v = recvValues[*r_iter].second;
878: const int numVals = recvValues[*r_iter].first;
879: int k = 0;
881: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
882: const int size = recvSection->getFiberDimension(s_iter.color());
884: if (size) {recvSection->updatePoint(s_iter.color(), &v[k]);}
885: k += size;
886: }
887: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
888: allocator.deallocate(v, numVals);
889: }
890: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
891: typename Section::value_type *v = sendValues[*r_iter].second;
892: const int numVals = sendValues[*r_iter].first;
894: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
895: allocator.deallocate(v, numVals);
896: }
897: //recvSection->view("After copy");
898: };
899: // Specialize to ArrowSection
900: template<typename SendOverlap, typename RecvOverlap>
901: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& sendSection, const Obj<UniformSection<MinimalArrow<int,int>,int> >& recvSection, const MPI_Datatype datatype = MPI_DATATYPE_NULL) {
902: typedef UniformSection<MinimalArrow<int,int>,int> Section;
903: typedef std::pair<int, typename Section::value_type *> allocPair;
904: const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
905: const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
906: MPIMover<typename Section::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
907: std::map<int, allocPair> sendValues;
908: std::map<int, allocPair> recvValues;
909: typename Section::alloc_type allocator;
911: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
912: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
913: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
915: // TODO: This should be const_iterator, but Sifter sucks
916: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
917: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
918: const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
919: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
920: int numVals = 0;
922: // TODO: This should be const_iterator, but Sifter sucks
923: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
924: for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
925: numVals += sendSection->getFiberDimension(typename Section::point_type(*c_iter, *d_iter));
926: }
927: }
928: typename Section::value_type *v = allocator.allocate(numVals);
929: int k = 0;
931: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
932: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
933: for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
934: const typename Section::point_type arrow(*c_iter, *d_iter);
935: const typename Section::value_type *vals = sendSection->restrictPoint(arrow);
936: const int dim = sendSection->getFiberDimension(arrow);
938: for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
939: }
940: }
941: sendValues[*r_iter] = allocPair(numVals, v);
942: vMover.send(*r_iter, numVals, v);
943: }
944: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
945: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
947: recvSection->allocatePoint();
948: // TODO: This should be const_iterator, but Sifter sucks
949: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
950: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
951: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
952: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
953: int numVals = 0;
955: // TODO: This should be const_iterator, but Sifter sucks
956: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
957: for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
958: numVals += recvSection->getFiberDimension(typename Section::point_type(s_iter.color(), t_iter.color()));
959: }
960: }
961: typename Section::value_type *v = allocator.allocate(numVals);
963: recvValues[*r_iter] = allocPair(numVals, v);
964: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
965: vMover.recv(*r_iter, numVals, v);
966: }
967: vMover.start();
968: vMover.end();
969: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
970: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
971: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
972: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
973: typename Section::value_type *v = recvValues[*r_iter].second;
974: const int numVals = recvValues[*r_iter].first;
975: int k = 0;
977: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
978: for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
979: const typename Section::point_type arrow(s_iter.color(), t_iter.color());
980: const int size = recvSection->getFiberDimension(arrow);
982: if (size) {recvSection->updatePoint(arrow, &v[k]);}
983: k += size;
984: }
985: }
986: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
987: allocator.deallocate(v, numVals);
988: }
989: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
990: typename Section::value_type *v = sendValues[*r_iter].second;
991: const int numVals = sendValues[*r_iter].first;
993: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
994: allocator.deallocate(v, numVals);
995: }
996: //recvSection->view("After copy");
997: };
998: // Specialize to a ConstantSection
999: #if 0
1000: template<typename SendOverlap, typename RecvOverlap, typename Value>
1001: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& recvSection) {
1002: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1003: };
1004: template<typename SendOverlap, typename RecvOverlap, typename Value>
1005: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& recvSection) {
1006: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1007: };
1008: #else
1009: template<typename SendOverlap, typename RecvOverlap, typename Value>
1010: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, Value> >& recvSection) {
1011: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1012: };
1013: template<typename SendOverlap, typename RecvOverlap, typename Value>
1014: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, Value> >& recvSection) {
1015: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1016: };
1017: #endif
1018: // Specialize to an IConstantSection
1019: template<typename SendOverlap, typename RecvOverlap, typename Value>
1020: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& recvSection) {
1021: // Why doesn't this work?
1022: // This supposed to be a copy, BUT filtered through the sendOverlap
1023: // However, an IConstant section does not allow filtration of its
1024: // chart. Therefore, you end up with either
1025: //
1026: // a) Too many items in the chart, copied from the remote sendSection
1027: // b) A chart mapped to the local numbering, which we do not want
1028: copyIConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1029: };
1030: // Specialize to an BaseSection/ConstantSection pair
1031: #if 0
1032: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1033: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSection<Sieve_> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, int> >& recvSection) {
1034: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1035: };
1036: #else
1037: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1038: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSection<Sieve_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
1039: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1040: };
1041: #endif
1042: // Specialize to an BaseSectionV/ConstantSection pair
1043: #if 0
1044: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1045: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSectionV<Sieve_> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, int> >& recvSection) {
1046: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1047: };
1048: #else
1049: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1050: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSectionV<Sieve_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
1051: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1052: };
1053: #endif
1054: // Specialize to an LabelBaseSection/ConstantSection pair
1055: #if 0
1056: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
1057: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<LabelBaseSection<Sieve_, Label_> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, int> >& recvSection) {
1058: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1059: };
1060: #else
1061: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
1062: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<LabelBaseSection<Sieve_, Label_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
1063: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1064: };
1065: #endif
1066: // Specialize to a ConstantSection for ArrowSection
1067: template<typename SendOverlap, typename RecvOverlap, typename Value>
1068: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<ConstantSection<MinimalArrow<typename SendOverlap::source_type,typename SendOverlap::source_type>, Value> >& sendSection, const Obj<ConstantSection<MinimalArrow<typename SendOverlap::source_type,typename SendOverlap::source_type>, Value> >& recvSection) {
1069: copyConstantArrow(sendOverlap, recvOverlap, sendSection, recvSection);
1070: };
1071: };
1072: class BinaryFusion {
1073: public:
1074: template<typename OverlapSection, typename RecvOverlap, typename Section, typename Function>
1075: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, Function binaryOp) {
1076: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1077: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1079: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1080: // TODO: This must become a more general iterator over colors
1081: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1082: // Just taking the first value
1083: const typename Section::point_type& localPoint = *p_iter;
1084: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1085: const typename OverlapSection::value_type *overlapValues = overlapSection->restrictPoint(remotePoint);
1086: const typename Section::value_type *localValues = section->restrictPoint(localPoint);
1087: const int dim = section->getFiberDimension(localPoint);
1088: // TODO: optimize allocation
1089: typename Section::value_type *values = new typename Section::value_type[dim];
1091: for(int d = 0; d < dim; ++d) {
1092: values[d] = binaryOp(overlapValues[d], localValues[d]);
1093: }
1094: section->updatePoint(localPoint, values);
1095: delete [] values;
1096: }
1097: };
1098: };
1099: class ReplacementBinaryFusion {
1100: public:
1101: template<typename OverlapSection, typename RecvOverlap, typename Section>
1102: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1103: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1104: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1106: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1107: // TODO: This must become a more general iterator over colors
1108: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1109: // Just taking the first value
1110: const typename Section::point_type& localPoint = *p_iter;
1111: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1113: section->update(localPoint, overlapSection->restrictPoint(remotePoint));
1114: }
1115: };
1116: };
1117: class AdditiveBinaryFusion {
1118: public:
1119: template<typename OverlapSection, typename RecvOverlap, typename Section>
1120: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1121: typedef typename Section::point_type point_type;
1122: typedef typename OverlapSection::point_type overlap_point_type;
1123: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1124: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1126: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1127: // TODO: This must become a more general iterator over colors
1128: const typename Section::point_type& localPoint = *p_iter;
1129: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1130: const typename RecvOverlap::coneSequence::iterator cEnd = points->end();
1132: for(typename RecvOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != cEnd; ++c_iter) {
1133: const int rank = *c_iter;
1134: const point_type& remotePoint = c_iter.color();
1136: section->updateAddPoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1137: }
1138: }
1139: };
1140: };
1141: class InsertionBinaryFusion {
1142: public:
1143: // Insert the overlapSection values into section along recvOverlap
1144: template<typename OverlapSection, typename RecvOverlap, typename Section>
1145: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1146: typedef typename Section::point_type point_type;
1147: typedef typename OverlapSection::point_type overlap_point_type;
1148: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1149: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1151: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1152: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1153: const point_type& localPoint = *p_iter;
1154: const int rank = *points->begin();
1155: const point_type& remotePoint = points->begin().color();
1157: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1158: if (!section->hasPoint(localPoint)) {
1159: std::cout <<"["<<section->commRank()<<"]: Destination section does not have local point " << localPoint << " remote point " << remotePoint << " fiber dim " << overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)) << std::endl;
1160: }
1161: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1162: }
1163: }
1164: if (rPoints->size()) {section->allocatePoint();}
1165: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1166: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1167: const point_type& localPoint = *p_iter;
1168: const int rank = *points->begin();
1169: const point_type& remotePoint = points->begin().color();
1171: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1172: section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1173: }
1174: }
1175: };
1176: // Specialize to ArrowSection
1177: template<typename OverlapSection, typename RecvOverlap>
1178: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& section) {
1179: typedef UniformSection<MinimalArrow<int,int>,int> Section;
1180: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1181: const typename RecvOverlap::traits::baseSequence::iterator rBegin = rPoints->begin();
1182: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1184: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1185: const Obj<typename RecvOverlap::coneSequence>& sources = recvOverlap->cone(*p_iter);
1186: const typename RecvOverlap::target_type localSource = *p_iter;
1187: const typename RecvOverlap::target_type remoteSource = sources->begin().color();
1189: for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1190: const Obj<typename RecvOverlap::coneSequence>& targets = recvOverlap->cone(*q_iter);
1191: const typename RecvOverlap::target_type localTarget = *q_iter;
1192: const typename RecvOverlap::target_type remoteTarget = targets->begin().color();
1193: const typename Section::point_type localPoint(localSource, localTarget);
1194: const typename OverlapSection::point_type remotePoint(remoteSource, remoteTarget);
1196: if (overlapSection->hasPoint(remotePoint)) {section->setFiberDimension(localPoint, overlapSection->getFiberDimension(remotePoint));}
1197: }
1198: }
1199: if (rPoints->size()) {section->allocatePoint();}
1200: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1201: const Obj<typename RecvOverlap::coneSequence>& sources = recvOverlap->cone(*p_iter);
1202: const typename RecvOverlap::target_type localSource = *p_iter;
1203: const typename RecvOverlap::target_type remoteSource = sources->begin().color();
1205: for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1206: const Obj<typename RecvOverlap::coneSequence>& targets = recvOverlap->cone(*q_iter);
1207: const typename RecvOverlap::target_type localTarget = *q_iter;
1208: const typename RecvOverlap::target_type remoteTarget = targets->begin().color();
1209: const typename Section::point_type localPoint(localSource, localTarget);
1210: const typename OverlapSection::point_type remotePoint(remoteSource, remoteTarget);
1211:
1212: if (overlapSection->hasPoint(remotePoint)) {section->updatePoint(localPoint, overlapSection->restrictPoint(remotePoint));}
1213: }
1214: }
1215: };
1216: // Specialize to the Sieve
1217: template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1218: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<Sieve<Point,Point,int> >& sieve) {
1219: typedef typename OverlapSection::point_type overlap_point_type;
1220: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1221: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1223: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1224: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1225: const Point& localPoint = *p_iter;
1226: const int rank = *points->begin();
1227: const Point& remotePoint = points->begin().color();
1228: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1229: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1230: int c = 0;
1232: sieve->clearCone(localPoint);
1233: for(int i = 0; i < size; ++i, ++c) {sieve->addCone(renumbering[values[i]], localPoint, c);}
1234: }
1235: };
1236: // Specialize to the ISieve
1237: template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1238: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<IFSieve<Point> >& sieve) {
1239: typedef typename OverlapSection::point_type overlap_point_type;
1240: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1241: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1242: int maxSize = 0;
1244: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1245: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1246: const Point& localPoint = *p_iter;
1247: const int rank = *points->begin();
1248: const Point& remotePoint = points->begin().color();
1249: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1250: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1252: sieve->setConeSize(localPoint, size);
1253: ///for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i]], 1);}
1254: for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
1255: maxSize = std::max(maxSize, size);
1256: }
1257: sieve->allocate();
1258: ///typename OverlapSection::value_type *localValues = new typename OverlapSection::value_type[maxSize];
1259: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1260: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1262: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1263: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1264: const Point& localPoint = *p_iter;
1265: const int rank = *points->begin();
1266: const Point& remotePoint = points->begin().color();
1267: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1268: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1270: ///for(int i = 0; i < size; ++i) {localValues[i] = renumbering[values[i]];}
1271: for(int i = 0; i < size; ++i) {
1272: localValues[i] = renumbering[values[i].first];
1273: localOrientation[i] = values[i].second;
1274: }
1275: sieve->setCone(localValues, localPoint);
1276: sieve->setConeOrientation(localOrientation, localPoint);
1277: }
1278: delete [] localValues;
1279: delete [] localOrientation;
1280: };
1281: template<typename OverlapSection, typename RecvOverlap, typename Point>
1282: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<IFSieve<Point> >& sieve) {
1283: typedef typename OverlapSection::point_type overlap_point_type;
1284: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1285: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1286: int maxSize = 0;
1288: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1289: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1290: const Point& localPoint = *p_iter;
1291: const int rank = *points->begin();
1292: const Point& remotePoint = points->begin().color();
1293: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1294: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1296: sieve->setConeSize(localPoint, size);
1297: for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
1298: maxSize = std::max(maxSize, size);
1299: }
1300: sieve->allocate();
1301: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1302: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1304: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1305: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1306: const Point& localPoint = *p_iter;
1307: const int rank = *points->begin();
1308: const Point& remotePoint = points->begin().color();
1309: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1310: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1312: for(int i = 0; i < size; ++i) {
1313: localValues[i] = values[i].first;
1314: localOrientation[i] = values[i].second;
1315: }
1316: sieve->setCone(localValues, localPoint);
1317: sieve->setConeOrientation(localOrientation, localPoint);
1318: }
1319: delete [] localValues;
1320: delete [] localOrientation;
1321: };
1322: // Generic
1323: template<typename OverlapSection, typename RecvOverlap, typename Section, typename Bundle>
1324: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, const Obj<Bundle>& bundle) {
1325: typedef typename OverlapSection::point_type overlap_point_type;
1326: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1327: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1329: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1330: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1331: const typename Section::point_type& localPoint = *p_iter;
1332: const int rank = *points->begin();
1333: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1335: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1336: }
1337: bundle->allocate(section);
1338: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1339: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1340: const typename Section::point_type& localPoint = *p_iter;
1341: const int rank = *points->begin();
1342: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1344: section->update(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1345: }
1346: };
1347: };
1348: }
1349: }
1351: #endif