Actual source code: ParallelMapping.hh
1: #ifndef included_ALE_ParallelMapping_hh
2: #define included_ALE_ParallelMapping_hh
4: #ifndef included_ALE_BasicCommunication_hh
5: #include <sieve/BasicCommunication.hh>
6: #endif
8: #ifndef included_ALE_IField_hh
9: #include <sieve/IField.hh>
10: #endif
12: #ifndef included_ALE_Sections_hh
13: #include <sieve/Sections.hh>
14: #endif
16: #include <functional>
17: #include <valarray>
19: namespace ALE {
20: template<class _Tp>
21: struct Identity : public std::unary_function<_Tp,_Tp>
22: {
23: _Tp& operator()(_Tp& x) const {return x;}
24: const _Tp& operator()(const _Tp& x) const {return x;}
25: };
27: template<class _Tp>
28: struct IsEqual : public std::unary_function<_Tp, bool>, public std::binary_function<_Tp, _Tp, bool>
29: {
30: const _Tp& x;
31: IsEqual(const _Tp& x) : x(x) {};
32: bool operator()(_Tp& y) const {return x == y;}
33: bool operator()(const _Tp& y) const {return x == y;}
34: bool operator()(_Tp& y, _Tp& dummy) const {return x == y;}
35: bool operator()(const _Tp& y, const _Tp& dummy) const {return x == y;}
36: };
38: // Creates new global point names and renames local points globally
39: template<typename Point>
40: class PointFactory : ALE::ParallelObject {
41: public:
42: typedef Point point_type;
43: typedef std::map<point_type,point_type> renumbering_type;
44: typedef std::map<int,std::map<point_type,point_type> > remote_renumbering_type;
45: protected:
46: point_type originalMax;
47: point_type currentMax;
48: renumbering_type renumbering;
49: renumbering_type invRenumbering;
50: remote_renumbering_type remoteRenumbering;
51: protected:
52: PointFactory(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), originalMax(-1) {};
53: public:
54: ~PointFactory() {};
55: public:
56: static PointFactory& singleton(MPI_Comm comm, const point_type& maxPoint, const int debug = 0, bool cleanup = false) {
57: static PointFactory *_singleton = NULL;
59: if (cleanup) {
60: if (debug) {std::cout << "Destroying PointFactory" << std::endl;}
61: if (_singleton) {delete _singleton;}
62: _singleton = NULL;
63: } else if (_singleton == NULL) {
64: if (debug) {std::cout << "Creating new PointFactory" << std::endl;}
65: _singleton = new PointFactory(comm, debug);
66: _singleton->setMax(maxPoint);
67: }
68: return *_singleton;
69: };
70: void setMax(const point_type& maxPoint) {
71: PetscErrorCode MPI_Allreduce((void *) &maxPoint, &this->originalMax, 1, MPI_INT, MPI_MAX, this->comm());CHKERRXX(ierr);
72: ++this->originalMax;
73: this->currentMax = this->originalMax;
74: };
75: void clear() {
76: this->currentMax = this->originalMax;
77: this->renumbering.clear();
78: this->invRenumbering.clear();
79: };
80: public:
81: template<typename Iterator>
82: void renumberPoints(const Iterator& begin, const Iterator& end) {
83: renumberPoints(begin, end, Identity<typename Iterator::value_type>());
84: }
85: template<typename Iterator, typename KeyExtractor>
86: void renumberPoints(const Iterator& begin, const Iterator& end, const KeyExtractor& ex) {
87: int numPoints = 0, numGlobalPoints, firstPoint;
89: for(Iterator p_iter = begin; p_iter != end; ++p_iter) ++numPoints;
90: MPI_Allreduce(&numPoints, &numGlobalPoints, 1, MPI_INT, MPI_SUM, this->comm());
91: MPI_Scan(&numPoints, &firstPoint, 1, MPI_INT, MPI_SUM, this->comm());
92: firstPoint += this->currentMax - numPoints;
93: this->currentMax += numGlobalPoints;
94: for(Iterator p_iter = begin; p_iter != end; ++p_iter, ++firstPoint) {
95: if (this->debug()) {std::cout << "["<<this->commRank()<<"]: New point " << ex(*p_iter) << " --> " << firstPoint << std::endl;}
96: this->renumbering[firstPoint] = ex(*p_iter);
97: this->invRenumbering[ex(*p_iter)] = firstPoint;
98: }
99: }
100: public:
101: void constructRemoteRenumbering() {
102: const int localSize = this->renumbering.size();
103: int *remoteSizes = new int[this->commSize()];
104: int *localMap = new int[localSize*2];
105: int *recvCounts = new int[this->commSize()];
106: int *displs = new int[this->commSize()];
108: // Populate arrays
109: int r = 0;
110: for(typename renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter, ++r) {
111: localMap[r*2+0] = r_iter->first;
112: localMap[r*2+1] = r_iter->second;
113: }
114: // Communicate renumbering sizes
115: MPI_Allgather((void*) &localSize, 1, MPI_INT, remoteSizes, 1, MPI_INT, this->comm());
116: for(int p = 0; p < this->commSize(); ++p) {
117: recvCounts[p] = remoteSizes[p]*2;
118: if (p == 0) {
119: displs[p] = 0;
120: } else {
121: displs[p] = displs[p-1] + recvCounts[p-1];
122: }
123: }
124: int *remoteMaps = new int[displs[this->commSize()-1]+recvCounts[this->commSize()-1]];
125: // Communicate renumberings
126: MPI_Allgatherv(localMap, localSize*2, MPI_INT, remoteMaps, recvCounts, displs, MPI_INT, this->comm());
127: // Populate maps
128: for(int p = 0; p < this->commSize(); ++p) {
129: if (p == this->commRank()) continue;
130: int offset = displs[p];
132: for(int r = 0; r < remoteSizes[p]; ++r) {
133: this->remoteRenumbering[p][remoteMaps[r*2+0+offset]] = remoteMaps[r*2+1+offset];
134: if (this->debug()) {std::cout << "["<<this->commRank()<<"]: Remote renumbering["<<p<<"] " << remoteMaps[r*2+0+offset] << " --> " << remoteMaps[r*2+1+offset] << std::endl;}
135: }
136: }
137: // Cleanup
138: delete [] recvCounts;
139: delete [] displs;
140: delete [] localMap;
141: delete [] remoteMaps;
142: delete [] remoteSizes;
143: };
144: public:
145: // global point --> local point
146: renumbering_type& getRenumbering() {
147: return this->renumbering;
148: };
149: // local point --> global point
150: renumbering_type& getInvRenumbering() {
151: return this->invRenumbering;
152: };
153: // rank --> global point --> local point
154: remote_renumbering_type& getRemoteRenumbering() {
155: return this->remoteRenumbering;
156: };
157: };
159: template<typename Alloc_ = malloc_allocator<int> >
160: class OverlapBuilder {
161: public:
162: typedef Alloc_ alloc_type;
163: protected:
164: template<typename T>
165: struct lessPair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
166: bool operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
167: return x.first < y.first;
168: }
169: };
170: template<typename T>
171: struct mergePair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
172: std::pair<T,std::pair<T,T> > operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
173: return std::pair<T,std::pair<T,T> >(x.first, std::pair<T,T>(x.second, y.second));
174: }
175: };
176: template<typename _InputIterator1, typename _InputIterator2, typename _OutputIterator, typename _Compare, typename _Merge>
177: static _OutputIterator set_intersection_merge(_InputIterator1 __first1, _InputIterator1 __last1,
178: _InputIterator2 __first2, _InputIterator2 __last2,
179: _OutputIterator __result, _Compare __comp, _Merge __merge)
180: {
181: while(__first1 != __last1 && __first2 != __last2) {
182: if (__comp(*__first1, *__first2))
183: ++__first1;
184: else if (__comp(*__first2, *__first1))
185: ++__first2;
186: else
187: {
188: *__result = __merge(*__first1, *__first2);
189: ++__first1;
190: ++__first2;
191: ++__result;
192: }
193: }
194: return __result;
195: }
196: public:
197: template<typename Sequence, typename Renumbering, typename SendOverlap, typename RecvOverlap>
198: static void constructOverlap(const Sequence& points, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
199: typedef typename SendOverlap::source_type point_type;
200: typedef std::pair<point_type,point_type> pointPair;
201: typedef std::pair<point_type,pointPair> pointTriple;
202: alloc_type allocator;
203: typename alloc_type::template rebind<point_type>::other point_allocator;
204: typename alloc_type::template rebind<pointPair>::other pointPair_allocator;
205: const MPI_Comm comm = sendOverlap->comm();
206: const int commSize = sendOverlap->commSize();
207: const int commRank = sendOverlap->commRank();
208: point_type *sendBuf = point_allocator.allocate(points.size()*2);
209: for(size_t i = 0; i < points.size()*2; ++i) {point_allocator.construct(sendBuf+i, point_type());}
210: int size = 0;
211: const int debug = sendOverlap->debug();
212: for(typename Sequence::const_iterator l_iter = points.begin(); l_iter != points.end(); ++l_iter) {
213: if (debug) {std::cout << "["<<commRank<<"]Send point["<<size<<"]: " << *l_iter << " " << renumbering[*l_iter] << std::endl;}
214: sendBuf[size++] = *l_iter;
215: sendBuf[size++] = renumbering[*l_iter];
216: }
217: int *sizes = allocator.allocate(commSize*3+2); // [size] The number of points coming from each process
218: for(int i = 0; i < commSize*3+2; ++i) {allocator.construct(sizes+i, 0);}
219: int *offsets = sizes+commSize; // [size+1] Prefix sums for sizes
220: int *oldOffs = offsets+commSize+1; // [size+1] Temporary storage
221: pointPair *remotePoints = NULL; // The points from each process
222: int *remoteRanks = NULL; // The rank and number of overlap points of each process that overlaps another
223: int numRemotePoints = 0;
224: int numRemoteRanks = 0;
226: // Change to Allgather() for the correct binning algorithm
227: MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, comm);
228: if (commRank == 0) {
229: offsets[0] = 0;
230: for(int p = 1; p <= commSize; p++) {
231: offsets[p] = offsets[p-1] + sizes[p-1];
232: }
233: numRemotePoints = offsets[commSize];
234: remotePoints = pointPair_allocator.allocate(numRemotePoints/2);
235: for(int i = 0; i < numRemotePoints/2; ++i) {pointPair_allocator.construct(remotePoints+i, pointPair());}
236: }
237: MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, comm);
238: for(size_t i = 0; i < points.size(); ++i) {point_allocator.destroy(sendBuf+i);}
239: point_allocator.deallocate(sendBuf, points.size());
240: std::map<int, std::map<int, std::set<pointTriple> > > overlapInfo; // Maps (p,q) to their set of overlap points
242: if (commRank == 0) {
243: for(int p = 0; p <= commSize; p++) {
244: offsets[p] /= 2;
245: }
246: for(int p = 0; p < commSize; p++) {
247: std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]], lessPair<point_type>());
248: }
249: for(int p = 0; p <= commSize; p++) {
250: oldOffs[p] = offsets[p];
251: }
252: for(int p = 0; p < commSize; p++) {
253: for(int q = 0; q < commSize; q++) {
254: if (p == q) continue;
255: set_intersection_merge(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
256: &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
257: std::insert_iterator<std::set<pointTriple> >(overlapInfo[p][q], overlapInfo[p][q].begin()),
258: lessPair<point_type>(), mergePair<point_type>());
259: }
260: sizes[p] = overlapInfo[p].size()*2;
261: offsets[p+1] = offsets[p] + sizes[p];
262: }
263: numRemoteRanks = offsets[commSize];
264: remoteRanks = allocator.allocate(numRemoteRanks);
265: for(int i = 0; i < numRemoteRanks; ++i) {allocator.construct(remoteRanks+i, 0);}
266: int k = 0;
267: for(int p = 0; p < commSize; p++) {
268: for(typename std::map<int, std::set<pointTriple> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
269: remoteRanks[k*2] = r_iter->first;
270: remoteRanks[k*2+1] = r_iter->second.size();
271: k++;
272: }
273: }
274: }
275: int numOverlaps; // The number of processes overlapping this process
276: MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, comm);
277: int *overlapRanks = allocator.allocate(numOverlaps); // The rank and overlap size for each overlapping process
278: for(int i = 0; i < numOverlaps; ++i) {allocator.construct(overlapRanks+i, 0);}
279: MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, comm);
280: point_type *sendPoints = NULL; // The points to send to each process
281: int numSendPoints = 0;
282: if (commRank == 0) {
283: for(int p = 0, k = 0; p < commSize; p++) {
284: sizes[p] = 0;
285: for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
286: sizes[p] += remoteRanks[k*2+1]*2;
287: k++;
288: }
289: offsets[p+1] = offsets[p] + sizes[p];
290: }
291: numSendPoints = offsets[commSize];
292: sendPoints = point_allocator.allocate(numSendPoints*2);
293: for(int i = 0; i < numSendPoints*2; ++i) {point_allocator.construct(sendPoints+i, point_type());}
294: for(int p = 0, k = 0; p < commSize; p++) {
295: for(typename std::map<int, std::set<pointTriple> >::const_iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
296: int rank = r_iter->first;
297: for(typename std::set<pointTriple>::const_iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
298: sendPoints[k++] = p_iter->first;
299: sendPoints[k++] = p_iter->second.second;
300: if (debug) {std::cout << "["<<commRank<<"]Sending points " << p_iter->first << " " << p_iter->second.second << " to rank " << rank << std::endl;}
301: }
302: }
303: }
304: }
305: int numOverlapPoints = 0;
306: for(int r = 0; r < numOverlaps/2; r++) {
307: numOverlapPoints += overlapRanks[r*2+1];
308: }
309: point_type *overlapPoints = point_allocator.allocate(numOverlapPoints*2);
310: for(int i = 0; i < numOverlapPoints*2; ++i) {point_allocator.construct(overlapPoints+i, point_type());}
311: MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints*2, MPI_INT, 0, comm);
313: for(int r = 0, k = 0; r < numOverlaps/2; r++) {
314: int rank = overlapRanks[r*2];
316: for(int p = 0; p < overlapRanks[r*2+1]; p++) {
317: point_type point = overlapPoints[k++];
318: point_type remotePoint = overlapPoints[k++];
320: if (debug) {std::cout << "["<<commRank<<"]Matched up remote point " << remotePoint << "("<<point<<") to local " << renumbering[point] << std::endl;}
321: sendOverlap->addArrow(renumbering[point], rank, remotePoint);
322: recvOverlap->addArrow(rank, renumbering[point], remotePoint);
323: }
324: }
326: for(int i = 0; i < numOverlapPoints; ++i) {point_allocator.destroy(overlapPoints+i);}
327: point_allocator.deallocate(overlapPoints, numOverlapPoints);
328: for(int i = 0; i < numOverlaps; ++i) {allocator.destroy(overlapRanks+i);}
329: allocator.deallocate(overlapRanks, numOverlaps);
330: for(int i = 0; i < commSize*3+2; ++i) {allocator.destroy(sizes+i);}
331: allocator.deallocate(sizes, commSize*3+2);
332: if (commRank == 0) {
333: for(int i = 0; i < numRemoteRanks; ++i) {allocator.destroy(remoteRanks+i);}
334: allocator.deallocate(remoteRanks, numRemoteRanks);
335: for(int i = 0; i < numRemotePoints; ++i) {pointPair_allocator.destroy(remotePoints+i);}
336: pointPair_allocator.deallocate(remotePoints, numRemotePoints);
337: for(int i = 0; i < numSendPoints; ++i) {point_allocator.destroy(sendPoints+i);}
338: point_allocator.deallocate(sendPoints, numSendPoints);
339: }
340: // TODO: Rewrite above to use optimized construction
341: sendOverlap->assemble();
342: recvOverlap->assemble();
343: }
344: };
345: namespace Pullback {
346: class SimpleCopy {
347: public:
348: // Copy the overlap section to the related processes
349: // This version is for Constant sections, meaning the same, single value over all points
350: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
351: static void copyConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
352: MPIMover<char> pMover(sendSection->comm(), sendSection->debug());
353: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
354: std::map<int, char *> sendPoints;
355: std::map<int, char *> recvPoints;
356: typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
357: typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;
359: const typename SendOverlap::baseSequence::iterator sBegin = sendOverlap->baseBegin();
360: const typename SendOverlap::baseSequence::iterator sEnd = sendOverlap->baseEnd();
361: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
363: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
364: const int pSize = sendOverlap->getConeSize(*r_iter);
365: const typename SendOverlap::coneSequence::iterator pBegin = sendOverlap->coneBegin(*r_iter);
366: const typename SendOverlap::coneSequence::iterator pEnd = sendOverlap->coneEnd(*r_iter);
367: char *v = sendAllocator.allocate(pSize);
368: int k = 0;
370: for(int i = 0; i < pSize; ++i) {sendAllocator.construct(v+i, 0);}
371: for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter, ++k) {
372: v[k] = (char) sendSection->hasPoint(*p_iter);
373: }
374: sendPoints[*r_iter] = v;
375: pMover.send(*r_iter, pSize, sendPoints[*r_iter]);
376: vMover.send(*r_iter, 2, sValues);
377: }
378: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
379: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
380: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
382: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
383: const int pSize = recvOverlap->getSupportSize(*r_iter);
384: char *v = recvAllocator.allocate(pSize);
386: for(int i = 0; i < pSize; ++i) {recvAllocator.construct(v+i, 0);}
387: recvPoints[*r_iter] = v;
388: pMover.recv(*r_iter, pSize, recvPoints[*r_iter]);
389: vMover.recv(*r_iter, 2, rValues);
390: }
391: pMover.start();
392: pMover.end();
393: vMover.start();
394: vMover.end();
395: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
396: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
397: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
398: const char *v = recvPoints[*r_iter];
399: int k = 0;
401: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter, ++k) {
402: if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(*r_iter, s_iter.color()));}
403: }
404: }
405: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
406: sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->getConeSize(*r_iter));
407: }
408: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
409: recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->getSupportSize(*r_iter));
410: }
411: }
412: // Specialize to ArrowSection
413: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
414: static void copyConstantArrow(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
415: MPIMover<char> pMover(sendSection->comm(), sendSection->debug());
416: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
417: std::map<int, char *> sendPoints;
418: std::map<int, char *> recvPoints;
419: typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
420: typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;
422: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
423: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
424: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
426: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
427: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
428: const int pSize = sendOverlap->getConeSize(*r_iter);
429: const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
430: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
431: char *v = sendAllocator.allocate(pSize*pSize);
432: int k = 0;
434: for(size_t i = 0; i < pSize*pSize; ++i) {sendAllocator.construct(v+i, 0);}
435: for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
436: for(typename SendOverlap::coneSequence::iterator q_iter = pBegin; q_iter != pEnd; ++q_iter, ++k) {
437: v[k] = (char) sendSection->hasPoint(typename SendSection::point_type(*p_iter,*q_iter));
438: }
439: }
440: sendPoints[*r_iter] = v;
441: pMover.send(*r_iter, pSize*pSize, sendPoints[*r_iter]);
442: vMover.send(*r_iter, 2, sValues);
443: }
444: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
445: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
446: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
448: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
449: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
450: const int pSize = recvOverlap->getSupportSize(*r_iter);
451: char *v = recvAllocator.allocate(pSize*pSize);
453: for(size_t i = 0; i < pSize*pSize; ++i) {recvAllocator.construct(v+i, 0);}
454: recvPoints[*r_iter] = v;
455: pMover.recv(*r_iter, pSize*pSize, recvPoints[*r_iter]);
456: vMover.recv(*r_iter, 2, rValues);
457: }
458: pMover.start();
459: pMover.end();
460: vMover.start();
461: vMover.end();
462: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
463: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
464: const typename RecvOverlap::traits::supportSequence::iterator pBegin = points->begin();
465: const typename RecvOverlap::traits::supportSequence::iterator pEnd = points->end();
466: const char *v = recvPoints[*r_iter];
467: int k = 0;
469: for(typename RecvOverlap::traits::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
470: for(typename RecvOverlap::traits::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter, ++k) {
471: if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(s_iter.color(),t_iter.color()));}
472: }
473: }
474: }
475: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
476: sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->getConeSize(*r_iter)*sendOverlap->getConeSize(*r_iter));
477: }
478: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
479: recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->getSupportSize(*r_iter)*recvOverlap->getSupportSize(*r_iter));
480: }
481: }
482: // Copy the overlap section to the related processes
483: // This version is for IConstant sections, meaning the same, single value over all points
484: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
485: static void copyIConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
486: MPIMover<typename SendSection::point_type> pMover(sendSection->comm(), sendSection->debug());
487: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
488: std::map<int, typename SendSection::point_type *> sendPoints;
489: std::map<int, typename SendSection::point_type *> recvPoints;
490: typename SendSection::alloc_type::template rebind<typename SendSection::point_type>::other sendAllocator;
491: typename RecvSection::alloc_type::template rebind<typename SendSection::point_type>::other recvAllocator;
493: const Obj<typename SendOverlap::baseSequence> sRanks = sendOverlap->base();
494: const typename SendOverlap::baseSequence::iterator sEnd = sRanks->end();
495: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
497: for(typename SendOverlap::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
498: typename SendSection::point_type *v = sendAllocator.allocate(2);
500: for(size_t i = 0; i < 2; ++i) {sendAllocator.construct(v+i, 0);}
501: v[0] = sendSection->getChart().min();
502: v[1] = sendSection->getChart().max();
503: sendPoints[*r_iter] = v;
504: pMover.send(*r_iter, 2, sendPoints[*r_iter]);
505: vMover.send(*r_iter, 2, sValues);
506: ///std::cout << "["<<sendOverlap->commRank()<<"]Sending chart (" << v[0] << ", " << v[1] << ") with values (" << sValues[0] << ", " << sValues[1] << ") to process " << *r_iter << std::endl;
507: }
508: const Obj<typename RecvOverlap::capSequence> rRanks = recvOverlap->cap();
509: const typename RecvOverlap::capSequence::iterator rEnd = rRanks->end();
510: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
512: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
513: typename SendSection::point_type *v = recvAllocator.allocate(2);
515: for(size_t i = 0; i < 2; ++i) {recvAllocator.construct(v+i, 0);}
516: recvPoints[*r_iter] = v;
517: pMover.recv(*r_iter, 2, recvPoints[*r_iter]);
518: vMover.recv(*r_iter, 2, rValues);
519: }
520: pMover.start();
521: pMover.end();
522: vMover.start();
523: vMover.end();
525: typename SendSection::point_type min = -1;
526: typename SendSection::point_type max = -1;
528: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
529: const typename RecvSection::point_type *v = recvPoints[*r_iter];
530: typename SendSection::point_type newMin = v[0];
531: typename SendSection::point_type newMax = v[1]-1;
532: //int pSize = 0;
534: ///std::cout << "["<<recvOverlap->commRank()<<"]Received chart (" << v[0] << ", " << v[1] << ") from process " << *r_iter << std::endl;
535: #if 0
536: // Translate to local numbering
537: if (recvOverlap->support(*r_iter)->size()) {
538: while(!pSize) {
539: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMin);
540: pSize = points->size();
541: if (pSize) {
542: newMin = *points->begin();
543: } else {
544: newMin++;
545: }
546: }
547: pSize = 0;
548: while(!pSize) {
549: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMax);
550: pSize = points->size();
551: if (pSize) {
552: newMax = *points->begin();
553: } else {
554: newMax--;
555: }
556: }
557: }
558: std::cout << "["<<recvOverlap->commRank()<<"]Translated to chart (" << newMin << ", " << newMax+1 << ") from process " << *r_iter << std::endl;
559: #endif
560: // Update chart
561: if (min < 0) {
562: min = newMin;
563: max = newMax+1;
564: } else {
565: min = std::min(min, newMin);
566: max = std::max(max, (typename SendSection::point_type) (newMax+1));
567: }
568: }
569: if (!rRanks->size()) {min = max = 0;}
570: recvSection->setChart(typename RecvSection::chart_type(min, max));
571: for(typename SendOverlap::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
572: sendAllocator.deallocate(sendPoints[*r_iter], 2);
573: }
574: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
575: recvAllocator.deallocate(recvPoints[*r_iter], 2);
576: }
577: }
578: // Copy the overlap section to the related processes
579: // This version is for different sections, possibly with different data types
580: // TODO: Can cache MPIMover objects (like a VecScatter)
581: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
582: 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) {
583: typedef std::pair<int, typename SendSection::value_type *> allocPair;
584: typedef typename RecvSection::point_type recv_point_type;
585: const Obj<typename SendSection::atlas_type>& sendAtlas = sendSection->getAtlas();
586: const Obj<typename RecvSection::atlas_type>& recvAtlas = recvSection->getAtlas();
587: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
588: std::map<int, allocPair> sendValues;
589: std::map<int, allocPair> recvValues;
590: typename SendSection::alloc_type sendAllocator;
591: typename RecvSection::alloc_type recvAllocator;
593: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
594: const typename SendOverlap::baseSequence::iterator sBegin = sendOverlap->baseBegin();
595: const typename SendOverlap::baseSequence::iterator sEnd = sendOverlap->baseEnd();
597: // TODO: This should be const_iterator, but Sifter sucks
598: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
599: const typename SendOverlap::coneSequence::iterator pBegin = sendOverlap->coneBegin(*r_iter);
600: const typename SendOverlap::coneSequence::iterator pEnd = sendOverlap->coneEnd(*r_iter);
601: const int numPoints = sendOverlap->getConeSize(*r_iter);
602: std::valarray<typename SendOverlap::source_type> sortedPoints(numPoints);
603: int numVals = 0, p = 0;
605: // TODO: This should be const_iterator, but Sifter sucks
606: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter, ++p) {
607: numVals += sendSection->getFiberDimension(*c_iter);
608: sortedPoints[p] = *c_iter;
609: }
610: typename SendSection::value_type *v = sendAllocator.allocate(numVals);
611: int k = 0;
613: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
614: for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
615: for(p = 0; p < numPoints; ++p) {
616: const typename SendSection::value_type *vals = sendSection->restrictPoint(sortedPoints[p]);
618: for(int i = 0; i < sendSection->getFiberDimension(sortedPoints[p]); ++i, ++k) v[k] = vals[i];
619: }
620: sendValues[*r_iter] = allocPair(numVals, v);
621: vMover.send(*r_iter, numVals, v);
622: }
623: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
624: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
626: recvSection->allocatePoint();
627: // TODO: This should be const_iterator, but Sifter sucks
628: int maxVals = 0;
629: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
630: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
631: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
632: int numVals = 0;
634: // TODO: This should be const_iterator, but Sifter sucks
635: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
636: numVals += recvSection->getFiberDimension(recv_point_type(*r_iter, s_iter.color()));
637: }
638: typename SendSection::value_type *v = sendAllocator.allocate(numVals);
640: for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
641: recvValues[*r_iter] = allocPair(numVals, v);
642: vMover.recv(*r_iter, numVals, v);
643: maxVals = std::max(maxVals, numVals);
644: }
645: vMover.start();
646: vMover.end();
647: typename RecvSection::value_type *convertedValues = recvAllocator.allocate(maxVals);
648: for(int i = 0; i < maxVals; ++i) {recvAllocator.construct(convertedValues+i, 0);}
649: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
650: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
651: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
652: const int numPoints = recvOverlap->getSupportSize(*r_iter);
653: std::valarray<typename RecvOverlap::color_type> sortedPoints(numPoints);
654: typename SendSection::value_type *v = recvValues[*r_iter].second;
655: const int numVals = recvValues[*r_iter].first;
656: int k = 0, p = 0;
658: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter, ++p) {
659: sortedPoints[p] = s_iter.color();
660: }
661: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
663: //for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
664: for(p = 0; p < numPoints; ++p) {
665: const int size = recvSection->getFiberDimension(recv_point_type(*r_iter, sortedPoints[p]));
667: for(int i = 0; i < size; ++i) {convertedValues[i] = (typename RecvSection::value_type) v[k+i];}
668: if (size) {recvSection->updatePoint(recv_point_type(*r_iter, sortedPoints[p]), convertedValues);}
669: k += size;
670: }
671: assert(k == numVals);
672: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
673: }
674: for(int i = 0; i < maxVals; ++i) {recvAllocator.destroy(convertedValues+i);}
675: recvAllocator.deallocate(convertedValues, maxVals);
676: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
677: typename SendSection::value_type *v = sendValues[*r_iter].second;
678: const int numVals = sendValues[*r_iter].first;
680: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
681: sendAllocator.deallocate(v, numVals);
682: }
683: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
684: typename SendSection::value_type *v = recvValues[*r_iter].second;
685: const int numVals = recvValues[*r_iter].first;
687: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
688: sendAllocator.deallocate(v, numVals);
689: }
690: //recvSection->view("After copy");
691: }
692: // Copy the overlap section to the related processes
693: // This version is for sections with the same type
694: template<typename SendOverlap, typename RecvOverlap, typename Section>
695: 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) {
696: typedef std::pair<int, typename Section::value_type *> allocPair;
697: const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
698: const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
699: MPIMover<typename Section::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
700: std::map<int, allocPair> sendValues;
701: std::map<int, allocPair> recvValues;
702: typename Section::alloc_type allocator;
704: ///sendAtlas->view("Send Atlas in same type copy()");
705: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
706: ///recvAtlas->view("Recv Atlas after same type copy()");
707: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
708: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
710: // TODO: This should be const_iterator, but Sifter sucks
711: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
712: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
713: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
714: const int numPoints = sendOverlap->getConeSize(*r_iter);
715: std::valarray<typename SendOverlap::source_type> sortedPoints(numPoints);
716: int numVals = 0, p = 0;
718: // TODO: This should be const_iterator, but Sifter sucks
719: for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter, ++p) {
720: numVals += sendSection->getFiberDimension(*c_iter);
721: sortedPoints[p] = *c_iter;
722: }
723: typename Section::value_type *v = allocator.allocate(numVals);
724: int k = 0;
726: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
727: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
728: //for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
729: for(p = 0; p < numPoints; ++p) {
730: const typename Section::value_type *vals = sendSection->restrictPoint(sortedPoints[p]);
731: const int dim = sendSection->getFiberDimension(sortedPoints[p]);
733: for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
734: }
735: sendValues[*r_iter] = allocPair(numVals, v);
736: vMover.send(*r_iter, numVals, v);
737: }
738: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
739: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
741: recvSection->allocatePoint();
742: ///recvSection->view("Recv Section after same type copy() allocation");
743: // TODO: This should be const_iterator, but Sifter sucks
744: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
745: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
746: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
747: int numVals = 0;
749: // TODO: This should be const_iterator, but Sifter sucks
750: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
751: numVals += recvSection->getFiberDimension(s_iter.color());
752: }
753: typename Section::value_type *v = allocator.allocate(numVals);
755: recvValues[*r_iter] = allocPair(numVals, v);
756: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
757: vMover.recv(*r_iter, numVals, v);
758: }
759: vMover.start();
760: vMover.end();
761: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
762: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
763: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
764: const int numPoints = recvOverlap->getSupportSize(*r_iter);
765: std::valarray<typename RecvOverlap::color_type> sortedPoints(numPoints);
766: typename Section::value_type *v = recvValues[*r_iter].second;
767: const int numVals = recvValues[*r_iter].first;
768: int k = 0, p = 0;
770: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter, ++p) {
771: sortedPoints[p] = s_iter.color();
772: }
773: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
775: //for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
776: for(p = 0; p < numPoints; ++p) {
777: const int size = recvSection->getFiberDimension(sortedPoints[p]);
779: if (size) {recvSection->updatePoint(sortedPoints[p], &v[k]);}
780: k += size;
781: }
782: assert(k == numVals);
783: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
784: allocator.deallocate(v, numVals);
785: }
786: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
787: typename Section::value_type *v = sendValues[*r_iter].second;
788: const int numVals = sendValues[*r_iter].first;
790: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
791: allocator.deallocate(v, numVals);
792: }
793: //recvSection->view("After copy");
794: }
795: // Specialize to ArrowSection
796: template<typename SendOverlap, typename RecvOverlap>
797: 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) {
798: typedef UniformSection<MinimalArrow<int,int>,int> Section;
799: typedef std::pair<int, typename Section::value_type *> allocPair;
800: const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
801: const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
802: MPIMover<typename Section::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
803: std::map<int, allocPair> sendValues;
804: std::map<int, allocPair> recvValues;
805: typename Section::alloc_type allocator;
807: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
808: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
809: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
811: // TODO: This should be const_iterator, but Sifter sucks
812: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
813: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
814: const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
815: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
816: int numVals = 0;
818: // TODO: This should be const_iterator, but Sifter sucks
819: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
820: for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
821: numVals += sendSection->getFiberDimension(typename Section::point_type(*c_iter, *d_iter));
822: }
823: }
824: typename Section::value_type *v = allocator.allocate(numVals);
825: int k = 0;
827: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
828: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
829: for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
830: const typename Section::point_type arrow(*c_iter, *d_iter);
831: const typename Section::value_type *vals = sendSection->restrictPoint(arrow);
832: const int dim = sendSection->getFiberDimension(arrow);
834: for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
835: }
836: }
837: sendValues[*r_iter] = allocPair(numVals, v);
838: vMover.send(*r_iter, numVals, v);
839: }
840: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
841: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
843: recvSection->allocatePoint();
844: // TODO: This should be const_iterator, but Sifter sucks
845: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
846: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
847: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
848: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
849: int numVals = 0;
851: // TODO: This should be const_iterator, but Sifter sucks
852: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
853: for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
854: numVals += recvSection->getFiberDimension(typename Section::point_type(s_iter.color(), t_iter.color()));
855: }
856: }
857: typename Section::value_type *v = allocator.allocate(numVals);
859: recvValues[*r_iter] = allocPair(numVals, v);
860: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
861: vMover.recv(*r_iter, numVals, v);
862: }
863: vMover.start();
864: vMover.end();
865: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
866: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
867: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
868: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
869: typename Section::value_type *v = recvValues[*r_iter].second;
870: const int numVals = recvValues[*r_iter].first;
871: int k = 0;
873: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
874: for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
875: const typename Section::point_type arrow(s_iter.color(), t_iter.color());
876: const int size = recvSection->getFiberDimension(arrow);
878: if (size) {recvSection->updatePoint(arrow, &v[k]);}
879: k += size;
880: }
881: }
882: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
883: allocator.deallocate(v, numVals);
884: }
885: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
886: typename Section::value_type *v = sendValues[*r_iter].second;
887: const int numVals = sendValues[*r_iter].first;
889: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
890: allocator.deallocate(v, numVals);
891: }
892: //recvSection->view("After copy");
893: }
894: // Specialize to a ConstantSection
895: #if 0
896: template<typename SendOverlap, typename RecvOverlap, typename Value>
897: 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) {
898: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
899: };
900: template<typename SendOverlap, typename RecvOverlap, typename Value>
901: 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) {
902: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
903: };
904: #else
905: template<typename SendOverlap, typename RecvOverlap, typename Value>
906: 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) {
907: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
908: }
909: template<typename SendOverlap, typename RecvOverlap, typename Value>
910: 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) {
911: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
912: }
913: #endif
914: // Specialize to an IConstantSection
915: template<typename SendOverlap, typename RecvOverlap, typename Value>
916: 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) {
917: // Why doesn't this work?
918: // This supposed to be a copy, BUT filtered through the sendOverlap
919: // However, an IConstant section does not allow filtration of its
920: // chart. Therefore, you end up with either
921: //
922: // a) Too many items in the chart, copied from the remote sendSection
923: // b) A chart mapped to the local numbering, which we do not want
924: copyIConstant(sendOverlap, recvOverlap, sendSection, recvSection);
925: }
926: // Specialize to an BaseSection/ConstantSection pair
927: #if 0
928: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
929: 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) {
930: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
931: };
932: #else
933: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
934: 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) {
935: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
936: }
937: #endif
938: // Specialize to an BaseSectionV/ConstantSection pair
939: #if 0
940: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
941: 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) {
942: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
943: };
944: #else
945: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
946: 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) {
947: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
948: }
949: #endif
950: // Specialize to an LabelBaseSection/ConstantSection pair
951: #if 0
952: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
953: 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) {
954: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
955: };
956: #else
957: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
958: 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) {
959: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
960: }
961: #endif
962: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
963: static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<LabelBaseSectionV<Sieve_, Label_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
964: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
965: }
966: // Specialize to a ConstantSection for ArrowSection
967: template<typename SendOverlap, typename RecvOverlap, typename Value>
968: 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) {
969: copyConstantArrow(sendOverlap, recvOverlap, sendSection, recvSection);
970: }
971: };
972: class BinaryFusion {
973: public:
974: template<typename OverlapSection, typename RecvOverlap, typename Section, typename Function>
975: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, Function binaryOp) {
976: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
977: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
979: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
980: // TODO: This must become a more general iterator over colors
981: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
982: // Just taking the first value
983: const typename Section::point_type& localPoint = *p_iter;
984: const typename OverlapSection::point_type& remotePoint = points->begin().color();
985: const typename OverlapSection::value_type *overlapValues = overlapSection->restrictPoint(remotePoint);
986: const typename Section::value_type *localValues = section->restrictPoint(localPoint);
987: const int dim = section->getFiberDimension(localPoint);
988: // TODO: optimize allocation
989: typename Section::value_type *values = new typename Section::value_type[dim];
991: for(int d = 0; d < dim; ++d) {
992: values[d] = binaryOp(overlapValues[d], localValues[d]);
993: }
994: section->updatePoint(localPoint, values);
995: delete [] values;
996: }
997: }
998: };
999: class ReplacementBinaryFusion {
1000: public:
1001: template<typename OverlapSection, typename RecvOverlap, typename Section>
1002: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1003: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1004: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1006: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1007: // TODO: This must become a more general iterator over colors
1008: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1009: // Just taking the first value
1010: const typename Section::point_type& localPoint = *p_iter;
1011: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1013: section->update(localPoint, overlapSection->restrictPoint(remotePoint));
1014: }
1015: }
1016: };
1017: class AdditiveBinaryFusion {
1018: public:
1019: template<typename OverlapSection, typename RecvOverlap, typename Section>
1020: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1021: typedef typename Section::point_type point_type;
1022: typedef typename OverlapSection::point_type overlap_point_type;
1023: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
1024: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
1026: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1027: const int rank = *r_iter;
1028: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1029: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1031: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1032: const point_type& localPoint = *p_iter;
1033: const point_type& remotePoint = p_iter.color();
1035: section->updateAddPoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1036: }
1037: }
1038: }
1039: };
1040: class InsertionBinaryFusion {
1041: public:
1042: // Insert the overlapSection values into section along recvOverlap
1043: template<typename OverlapSection, typename RecvOverlap, typename Section>
1044: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1045: typedef typename Section::point_type point_type;
1046: typedef typename OverlapSection::point_type overlap_point_type;
1047: #if 0
1048: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1049: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1051: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1052: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1053: const point_type& localPoint = *p_iter;
1054: const int rank = *points->begin();
1055: const point_type& remotePoint = points->begin().color();
1057: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1058: if (!section->hasPoint(localPoint)) {
1059: 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;
1060: }
1061: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1062: }
1063: }
1064: if (rPoints->size()) {section->allocatePoint();}
1065: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1066: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1067: const point_type& localPoint = *p_iter;
1068: const int rank = *points->begin();
1069: const point_type& remotePoint = points->begin().color();
1071: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1072: section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1073: }
1074: }
1075: #else
1076: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
1077: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
1078: bool hasPoints = false;
1080: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1081: const int rank = *r_iter;
1082: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1083: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1085: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1086: const point_type& localPoint = *p_iter;
1087: const point_type& remotePoint = p_iter.color();
1089: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1090: if (!section->hasPoint(localPoint)) {
1091: 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;
1092: }
1093: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1094: }
1095: hasPoints = true;
1096: }
1097: }
1098: if (hasPoints) {section->allocatePoint();}
1099: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1100: const int rank = *r_iter;
1101: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1102: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1104: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1105: const point_type& localPoint = *p_iter;
1106: const point_type& remotePoint = p_iter.color();
1108: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1109: section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1110: }
1111: }
1112: }
1113: #endif
1114: }
1115: // Specialize to ArrowSection
1116: template<typename OverlapSection, typename RecvOverlap>
1117: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& section) {
1118: typedef UniformSection<MinimalArrow<int,int>,int> Section;
1119: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1120: const typename RecvOverlap::traits::baseSequence::iterator rBegin = rPoints->begin();
1121: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1123: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1124: const Obj<typename RecvOverlap::coneSequence>& sources = recvOverlap->cone(*p_iter);
1125: const typename RecvOverlap::target_type localSource = *p_iter;
1126: const typename RecvOverlap::target_type remoteSource = sources->begin().color();
1128: for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1129: const Obj<typename RecvOverlap::coneSequence>& targets = recvOverlap->cone(*q_iter);
1130: const typename RecvOverlap::target_type localTarget = *q_iter;
1131: const typename RecvOverlap::target_type remoteTarget = targets->begin().color();
1132: const typename Section::point_type localPoint(localSource, localTarget);
1133: const typename OverlapSection::point_type remotePoint(remoteSource, remoteTarget);
1135: if (overlapSection->hasPoint(remotePoint)) {section->setFiberDimension(localPoint, overlapSection->getFiberDimension(remotePoint));}
1136: }
1137: }
1138: if (rPoints->size()) {section->allocatePoint();}
1139: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1140: const Obj<typename RecvOverlap::coneSequence>& sources = recvOverlap->cone(*p_iter);
1141: const typename RecvOverlap::target_type localSource = *p_iter;
1142: const typename RecvOverlap::target_type remoteSource = sources->begin().color();
1144: for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1145: const Obj<typename RecvOverlap::coneSequence>& targets = recvOverlap->cone(*q_iter);
1146: const typename RecvOverlap::target_type localTarget = *q_iter;
1147: const typename RecvOverlap::target_type remoteTarget = targets->begin().color();
1148: const typename Section::point_type localPoint(localSource, localTarget);
1149: const typename OverlapSection::point_type remotePoint(remoteSource, remoteTarget);
1151: if (overlapSection->hasPoint(remotePoint)) {section->updatePoint(localPoint, overlapSection->restrictPoint(remotePoint));}
1152: }
1153: }
1154: }
1155: // Specialize to the Sieve
1156: template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1157: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<Sieve<Point,Point,int> >& sieve) {
1158: typedef typename OverlapSection::point_type overlap_point_type;
1159: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1160: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1162: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1163: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1164: const Point& localPoint = *p_iter;
1165: const int rank = *points->begin();
1166: const Point& remotePoint = points->begin().color();
1167: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1168: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1169: int c = 0;
1171: sieve->clearCone(localPoint);
1172: for(int i = 0; i < size; ++i, ++c) {sieve->addCone(renumbering[values[i]], localPoint, c);}
1173: }
1174: }
1175: // Specialize to the ISieve
1176: template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1177: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<IFSieve<Point> >& sieve) {
1178: typedef typename OverlapSection::point_type overlap_point_type;
1179: #if 0
1180: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1181: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1182: int maxSize = 0;
1184: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1185: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1186: const Point& localPoint = *p_iter;
1187: const int rank = *points->begin();
1188: const Point& remotePoint = points->begin().color();
1189: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1190: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1192: sieve->setConeSize(localPoint, size);
1193: ///for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i]], 1);}
1194: for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
1195: maxSize = std::max(maxSize, size);
1196: }
1197: sieve->allocate();
1198: ///typename OverlapSection::value_type *localValues = new typename OverlapSection::value_type[maxSize];
1199: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1200: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1202: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1203: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1204: const Point& localPoint = *p_iter;
1205: const int rank = *points->begin();
1206: const Point& remotePoint = points->begin().color();
1207: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1208: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1210: ///for(int i = 0; i < size; ++i) {localValues[i] = renumbering[values[i]];}
1211: for(int i = 0; i < size; ++i) {
1212: localValues[i] = renumbering[values[i].first];
1213: localOrientation[i] = values[i].second;
1214: }
1215: sieve->setCone(localValues, localPoint);
1216: sieve->setConeOrientation(localOrientation, localPoint);
1217: }
1218: delete [] localValues;
1219: delete [] localOrientation;
1220: #else
1221: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
1222: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
1223: int maxSize = 0;
1225: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1226: const int rank = *r_iter;
1227: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1228: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1230: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1231: const Point& localPoint = *p_iter;
1232: const Point& remotePoint = p_iter.color();
1233: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1234: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1236: sieve->setConeSize(localPoint, size);
1237: for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
1238: maxSize = std::max(maxSize, size);
1239: }
1240: }
1241: sieve->allocate();
1242: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1243: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1245: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1246: const int rank = *r_iter;
1247: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1248: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1250: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1251: const Point& localPoint = *p_iter;
1252: const Point& remotePoint = p_iter.color();
1253: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1254: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1256: for(int i = 0; i < size; ++i) {
1257: localValues[i] = renumbering[values[i].first];
1258: localOrientation[i] = values[i].second;
1259: }
1260: sieve->setCone(localValues, localPoint);
1261: sieve->setConeOrientation(localOrientation, localPoint);
1262: }
1263: }
1264: delete [] localValues;
1265: delete [] localOrientation;
1266: #endif
1267: }
1268: template<typename OverlapSection, typename RecvOverlap, typename Point>
1269: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<IFSieve<Point> >& sieve) {
1270: typedef typename OverlapSection::point_type overlap_point_type;
1271: #if 0
1272: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1273: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1274: int maxSize = 0;
1276: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1277: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1278: const Point& localPoint = *p_iter;
1279: const int rank = *points->begin();
1280: const Point& remotePoint = points->begin().color();
1281: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1282: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1284: sieve->setConeSize(localPoint, size);
1285: for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
1286: maxSize = std::max(maxSize, size);
1287: }
1288: sieve->allocate();
1289: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1290: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1292: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1293: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1294: const Point& localPoint = *p_iter;
1295: const int rank = *points->begin();
1296: const Point& remotePoint = points->begin().color();
1297: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1298: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1300: for(int i = 0; i < size; ++i) {
1301: localValues[i] = values[i].first;
1302: localOrientation[i] = values[i].second;
1303: }
1304: sieve->setCone(localValues, localPoint);
1305: sieve->setConeOrientation(localOrientation, localPoint);
1306: }
1307: delete [] localValues;
1308: delete [] localOrientation;
1309: #else
1310: const Obj<typename RecvOverlap::capSequence> rRanks = recvOverlap->cap();
1311: const typename RecvOverlap::capSequence::iterator rEnd = rRanks->end();
1312: int maxSize = 0;
1314: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
1315: const int rank = *r_iter;
1316: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
1317: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
1318: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
1320: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1321: const Point& localPoint = *p_iter;
1322: const Point& remotePoint = p_iter.color();
1323: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1324: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1326: sieve->setConeSize(localPoint, size);
1327: for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
1328: maxSize = std::max(maxSize, size);
1329: }
1330: }
1331: sieve->allocate();
1332: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1333: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1335: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
1336: const int rank = *r_iter;
1337: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
1338: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
1339: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
1341: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1342: const Point& localPoint = *p_iter;
1343: const Point& remotePoint = p_iter.color();
1344: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1345: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1347: for(int i = 0; i < size; ++i) {
1348: localValues[i] = values[i].first;
1349: localOrientation[i] = values[i].second;
1350: }
1351: sieve->setCone(localValues, localPoint);
1352: sieve->setConeOrientation(localOrientation, localPoint);
1353: }
1354: }
1355: delete [] localValues;
1356: delete [] localOrientation;
1357: #endif
1358: }
1359: // Generic
1360: template<typename OverlapSection, typename RecvOverlap, typename Section, typename Bundle>
1361: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, const Obj<Bundle>& bundle) {
1362: typedef typename OverlapSection::point_type overlap_point_type;
1363: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1364: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1366: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1367: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1368: const typename Section::point_type& localPoint = *p_iter;
1369: const int rank = *points->begin();
1370: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1372: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1373: }
1374: bundle->allocate(section);
1375: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1376: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1377: const typename Section::point_type& localPoint = *p_iter;
1378: const int rank = *points->begin();
1379: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1381: section->update(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1382: }
1383: }
1384: };
1385: class InterpolateMultipleFusion {
1386: public:
1387: // Interpolate the overlapSection values into section along recvOverlap
1388: template<typename OverlapSection, typename RecvOverlap, typename Section>
1389: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1390: typedef typename Section::point_type point_type;
1391: typedef typename Section::value_type value_type;
1392: typedef typename OverlapSection::point_type overlap_point_type;
1393: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1394: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1395: int maxFiberDim = -1;
1397: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1398: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1399: const typename RecvOverlap::coneSequence::iterator rpEnd = points->end();
1400: const point_type& localPoint = *p_iter;
1401: bool inOverlap = false;
1402: int fiberDim = -1;
1404: for(typename RecvOverlap::coneSequence::iterator rp_iter = points->begin(); rp_iter != rpEnd; ++rp_iter) {
1405: const int rank = *rp_iter;
1406: const point_type& remotePoint = rp_iter.color();
1408: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1409: inOverlap = true;
1410: fiberDim = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1411: break;
1412: }
1413: }
1414: if (inOverlap) {
1415: if (!section->hasPoint(localPoint)) {
1416: std::cout <<"["<<section->commRank()<<"]: Destination section does not have local point " << localPoint << " remote point " << (points->begin().color()) << " fiber dim " << fiberDim << std::endl;
1417: }
1418: section->setFiberDimension(localPoint, fiberDim);
1419: maxFiberDim = std::max(fiberDim, maxFiberDim);
1420: }
1421: }
1422: if (rPoints->size()) {section->allocatePoint();}
1423: value_type *interpolant = new value_type[maxFiberDim];
1425: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1426: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1427: const typename RecvOverlap::coneSequence::iterator rpEnd = points->end();
1428: const point_type& localPoint = *p_iter;
1429: bool inOverlap = false;
1430: int numArgs = 0;
1432: for(int d = 0; d < maxFiberDim; ++d) {interpolant[d] = 0.0;}
1433: for(typename RecvOverlap::coneSequence::iterator rp_iter = points->begin(); rp_iter != rpEnd; ++rp_iter) {
1434: const int rank = *rp_iter;
1435: const point_type& remotePoint = rp_iter.color();
1436: const overlap_point_type opoint(rank, remotePoint);
1438: if (overlapSection->hasPoint(opoint)) {
1439: const int fiberDim = overlapSection->getFiberDimension(opoint);
1440: const value_type *values = overlapSection->restrictPoint(opoint);
1442: // TODO: Include interpolation weights (stored in overlap)
1443: for(int d = 0; d < fiberDim; ++d) {
1444: interpolant[d] += values[d];
1445: }
1446: inOverlap = true;
1447: ++numArgs;
1448: }
1449: }
1450: if (inOverlap) {
1451: for(int d = 0; d < maxFiberDim; ++d) {interpolant[d] /= numArgs;}
1452: section->updatePoint(localPoint, interpolant);
1453: }
1454: }
1455: delete [] interpolant;
1456: }
1457: };
1458: }
1459: }
1461: #endif