Actual source code: ParallelMapping.hh
petsc-3.3-p7 2013-05-11
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: sendOverlap->assemblePoints();
344: recvOverlap->assemblePoints();
345: }
346: };
347: namespace Pullback {
348: class SimpleCopy {
349: public:
350: // Copy the overlap section to the related processes
351: // This version is for Constant sections, meaning the same, single value over all points
352: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
353: static void copyConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
354: MPIMover<char> pMover(sendSection->comm(), sendSection->debug());
355: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
356: std::map<int, char *> sendPoints;
357: std::map<int, char *> recvPoints;
358: typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
359: typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;
361: const typename SendOverlap::baseSequence::iterator sBegin = sendOverlap->baseBegin();
362: const typename SendOverlap::baseSequence::iterator sEnd = sendOverlap->baseEnd();
363: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
365: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
366: const int pSize = sendOverlap->getConeSize(*r_iter);
367: const typename SendOverlap::coneSequence::iterator pBegin = sendOverlap->coneBegin(*r_iter);
368: const typename SendOverlap::coneSequence::iterator pEnd = sendOverlap->coneEnd(*r_iter);
369: char *v = sendAllocator.allocate(pSize);
370: int k = 0;
372: for(int i = 0; i < pSize; ++i) {sendAllocator.construct(v+i, 0);}
373: for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter, ++k) {
374: v[k] = (char) sendSection->hasPoint(*p_iter);
375: }
376: sendPoints[*r_iter] = v;
377: pMover.send(*r_iter, pSize, sendPoints[*r_iter]);
378: vMover.send(*r_iter, 2, sValues);
379: }
380: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
381: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
382: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
384: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
385: const int pSize = recvOverlap->getSupportSize(*r_iter);
386: char *v = recvAllocator.allocate(pSize);
388: for(int i = 0; i < pSize; ++i) {recvAllocator.construct(v+i, 0);}
389: recvPoints[*r_iter] = v;
390: pMover.recv(*r_iter, pSize, recvPoints[*r_iter]);
391: vMover.recv(*r_iter, 2, rValues);
392: }
393: pMover.start();
394: pMover.end();
395: vMover.start();
396: vMover.end();
397: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
398: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
399: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
400: const char *v = recvPoints[*r_iter];
401: int k = 0;
403: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter, ++k) {
404: if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(*r_iter, s_iter.color()));}
405: }
406: }
407: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
408: sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->getConeSize(*r_iter));
409: }
410: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
411: recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->getSupportSize(*r_iter));
412: }
413: }
414: // Specialize to ArrowSection
415: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
416: static void copyConstantArrow(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
417: MPIMover<char> pMover(sendSection->comm(), sendSection->debug());
418: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
419: std::map<int, char *> sendPoints;
420: std::map<int, char *> recvPoints;
421: typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
422: typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;
424: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
425: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
426: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
428: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
429: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
430: const int pSize = sendOverlap->getConeSize(*r_iter);
431: const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
432: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
433: char *v = sendAllocator.allocate(pSize*pSize);
434: int k = 0;
436: for(size_t i = 0; i < pSize*pSize; ++i) {sendAllocator.construct(v+i, 0);}
437: for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
438: for(typename SendOverlap::coneSequence::iterator q_iter = pBegin; q_iter != pEnd; ++q_iter, ++k) {
439: v[k] = (char) sendSection->hasPoint(typename SendSection::point_type(*p_iter,*q_iter));
440: }
441: }
442: sendPoints[*r_iter] = v;
443: pMover.send(*r_iter, pSize*pSize, sendPoints[*r_iter]);
444: vMover.send(*r_iter, 2, sValues);
445: }
446: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
447: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
448: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
450: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
451: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
452: const int pSize = recvOverlap->getSupportSize(*r_iter);
453: char *v = recvAllocator.allocate(pSize*pSize);
455: for(size_t i = 0; i < pSize*pSize; ++i) {recvAllocator.construct(v+i, 0);}
456: recvPoints[*r_iter] = v;
457: pMover.recv(*r_iter, pSize*pSize, recvPoints[*r_iter]);
458: vMover.recv(*r_iter, 2, rValues);
459: }
460: pMover.start();
461: pMover.end();
462: vMover.start();
463: vMover.end();
464: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
465: const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
466: const typename RecvOverlap::traits::supportSequence::iterator pBegin = points->begin();
467: const typename RecvOverlap::traits::supportSequence::iterator pEnd = points->end();
468: const char *v = recvPoints[*r_iter];
469: int k = 0;
471: for(typename RecvOverlap::traits::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
472: for(typename RecvOverlap::traits::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter, ++k) {
473: if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(s_iter.color(),t_iter.color()));}
474: }
475: }
476: }
477: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
478: sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->getConeSize(*r_iter)*sendOverlap->getConeSize(*r_iter));
479: }
480: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
481: recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->getSupportSize(*r_iter)*recvOverlap->getSupportSize(*r_iter));
482: }
483: }
484: // Copy the overlap section to the related processes
485: // This version is for IConstant sections, meaning the same, single value over all points
486: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
487: static void copyIConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
488: MPIMover<typename SendSection::point_type> pMover(sendSection->comm(), sendSection->debug());
489: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
490: std::map<int, typename SendSection::point_type *> sendPoints;
491: std::map<int, typename SendSection::point_type *> recvPoints;
492: typename SendSection::alloc_type::template rebind<typename SendSection::point_type>::other sendAllocator;
493: typename RecvSection::alloc_type::template rebind<typename SendSection::point_type>::other recvAllocator;
495: const Obj<typename SendOverlap::baseSequence> sRanks = sendOverlap->base();
496: const typename SendOverlap::baseSequence::iterator sEnd = sRanks->end();
497: const typename SendSection::value_type *sValues = sendSection->restrictSpace();
499: for(typename SendOverlap::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
500: typename SendSection::point_type *v = sendAllocator.allocate(2);
502: for(size_t i = 0; i < 2; ++i) {sendAllocator.construct(v+i, 0);}
503: v[0] = sendSection->getChart().min();
504: v[1] = sendSection->getChart().max();
505: sendPoints[*r_iter] = v;
506: pMover.send(*r_iter, 2, sendPoints[*r_iter]);
507: vMover.send(*r_iter, 2, sValues);
508: ///std::cout << "["<<sendOverlap->commRank()<<"]Sending chart (" << v[0] << ", " << v[1] << ") with values (" << sValues[0] << ", " << sValues[1] << ") to process " << *r_iter << std::endl;
509: }
510: const Obj<typename RecvOverlap::capSequence> rRanks = recvOverlap->cap();
511: const typename RecvOverlap::capSequence::iterator rEnd = rRanks->end();
512: const typename RecvSection::value_type *rValues = recvSection->restrictSpace();
514: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
515: typename SendSection::point_type *v = recvAllocator.allocate(2);
517: for(size_t i = 0; i < 2; ++i) {recvAllocator.construct(v+i, 0);}
518: recvPoints[*r_iter] = v;
519: pMover.recv(*r_iter, 2, recvPoints[*r_iter]);
520: vMover.recv(*r_iter, 2, rValues);
521: }
522: pMover.start();
523: pMover.end();
524: vMover.start();
525: vMover.end();
527: typename SendSection::point_type min = -1;
528: typename SendSection::point_type max = -1;
530: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
531: const typename RecvSection::point_type *v = recvPoints[*r_iter];
532: typename SendSection::point_type newMin = v[0];
533: typename SendSection::point_type newMax = v[1]-1;
534: //int pSize = 0;
536: ///std::cout << "["<<recvOverlap->commRank()<<"]Received chart (" << v[0] << ", " << v[1] << ") from process " << *r_iter << std::endl;
537: #if 0
538: // Translate to local numbering
539: if (recvOverlap->support(*r_iter)->size()) {
540: while(!pSize) {
541: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMin);
542: pSize = points->size();
543: if (pSize) {
544: newMin = *points->begin();
545: } else {
546: newMin++;
547: }
548: }
549: pSize = 0;
550: while(!pSize) {
551: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMax);
552: pSize = points->size();
553: if (pSize) {
554: newMax = *points->begin();
555: } else {
556: newMax--;
557: }
558: }
559: }
560: std::cout << "["<<recvOverlap->commRank()<<"]Translated to chart (" << newMin << ", " << newMax+1 << ") from process " << *r_iter << std::endl;
561: #endif
562: // Update chart
563: if (min < 0) {
564: min = newMin;
565: max = newMax+1;
566: } else {
567: min = std::min(min, newMin);
568: max = std::max(max, (typename SendSection::point_type) (newMax+1));
569: }
570: }
571: if (!rRanks->size()) {min = max = 0;}
572: recvSection->setChart(typename RecvSection::chart_type(min, max));
573: for(typename SendOverlap::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
574: sendAllocator.deallocate(sendPoints[*r_iter], 2);
575: }
576: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
577: recvAllocator.deallocate(recvPoints[*r_iter], 2);
578: }
579: }
580: // Copy the overlap section to the related processes
581: // This version is for different sections, possibly with different data types
582: // TODO: Can cache MPIMover objects (like a VecScatter)
583: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
584: 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) {
585: typedef std::pair<int, typename SendSection::value_type *> allocPair;
586: typedef typename RecvSection::point_type recv_point_type;
587: const Obj<typename SendSection::atlas_type>& sendAtlas = sendSection->getAtlas();
588: const Obj<typename RecvSection::atlas_type>& recvAtlas = recvSection->getAtlas();
589: MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
590: std::map<int, allocPair> sendValues;
591: std::map<int, allocPair> recvValues;
592: typename SendSection::alloc_type sendAllocator;
593: typename RecvSection::alloc_type recvAllocator;
595: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
596: const typename SendOverlap::baseSequence::iterator sBegin = sendOverlap->baseBegin();
597: const typename SendOverlap::baseSequence::iterator sEnd = sendOverlap->baseEnd();
599: // TODO: This should be const_iterator, but Sifter sucks
600: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
601: const typename SendOverlap::coneSequence::iterator pBegin = sendOverlap->coneBegin(*r_iter);
602: const typename SendOverlap::coneSequence::iterator pEnd = sendOverlap->coneEnd(*r_iter);
603: const int numPoints = sendOverlap->getConeSize(*r_iter);
604: std::valarray<typename SendOverlap::source_type> sortedPoints(numPoints);
605: int numVals = 0, p = 0;
607: // TODO: This should be const_iterator, but Sifter sucks
608: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter, ++p) {
609: numVals += sendSection->getFiberDimension(*c_iter);
610: sortedPoints[p] = *c_iter;
611: }
612: typename SendSection::value_type *v = sendAllocator.allocate(numVals);
613: int k = 0;
615: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
616: for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
617: for(p = 0; p < numPoints; ++p) {
618: const typename SendSection::value_type *vals = sendSection->restrictPoint(sortedPoints[p]);
620: for(int i = 0; i < sendSection->getFiberDimension(sortedPoints[p]); ++i, ++k) v[k] = vals[i];
621: }
622: sendValues[*r_iter] = allocPair(numVals, v);
623: vMover.send(*r_iter, numVals, v);
624: }
625: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
626: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
628: recvSection->allocatePoint();
629: // TODO: This should be const_iterator, but Sifter sucks
630: int maxVals = 0;
631: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
632: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
633: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
634: int numVals = 0;
636: // TODO: This should be const_iterator, but Sifter sucks
637: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
638: numVals += recvSection->getFiberDimension(recv_point_type(*r_iter, s_iter.color()));
639: }
640: typename SendSection::value_type *v = sendAllocator.allocate(numVals);
642: for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
643: recvValues[*r_iter] = allocPair(numVals, v);
644: vMover.recv(*r_iter, numVals, v);
645: maxVals = std::max(maxVals, numVals);
646: }
647: vMover.start();
648: vMover.end();
649: typename RecvSection::value_type *convertedValues = recvAllocator.allocate(maxVals);
650: for(int i = 0; i < maxVals; ++i) {recvAllocator.construct(convertedValues+i, 0);}
651: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
652: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
653: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
654: const int numPoints = recvOverlap->getSupportSize(*r_iter);
655: std::valarray<typename RecvOverlap::color_type> sortedPoints(numPoints);
656: typename SendSection::value_type *v = recvValues[*r_iter].second;
657: const int numVals = recvValues[*r_iter].first;
658: int k = 0, p = 0;
660: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter, ++p) {
661: sortedPoints[p] = s_iter.color();
662: }
663: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
665: //for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
666: for(p = 0; p < numPoints; ++p) {
667: const int size = recvSection->getFiberDimension(recv_point_type(*r_iter, sortedPoints[p]));
669: for(int i = 0; i < size; ++i) {convertedValues[i] = (typename RecvSection::value_type) v[k+i];}
670: if (size) {recvSection->updatePoint(recv_point_type(*r_iter, sortedPoints[p]), convertedValues);}
671: k += size;
672: }
673: assert(k == numVals);
674: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
675: }
676: for(int i = 0; i < maxVals; ++i) {recvAllocator.destroy(convertedValues+i);}
677: recvAllocator.deallocate(convertedValues, maxVals);
678: for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
679: typename SendSection::value_type *v = sendValues[*r_iter].second;
680: const int numVals = sendValues[*r_iter].first;
682: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
683: sendAllocator.deallocate(v, numVals);
684: }
685: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
686: typename SendSection::value_type *v = recvValues[*r_iter].second;
687: const int numVals = recvValues[*r_iter].first;
689: for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
690: sendAllocator.deallocate(v, numVals);
691: }
692: //recvSection->view("After copy");
693: }
694: // Copy the overlap section to the related processes
695: // This version is for sections with the same type
696: template<typename SendOverlap, typename RecvOverlap, typename Section>
697: 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) {
698: typedef std::pair<int, typename Section::value_type *> allocPair;
699: const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
700: const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
701: MPIMover<typename Section::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
702: std::map<int, allocPair> sendValues;
703: std::map<int, allocPair> recvValues;
704: typename Section::alloc_type allocator;
706: ///sendAtlas->view("Send Atlas in same type copy()");
707: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
708: ///recvAtlas->view("Recv Atlas after same type copy()");
709: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
710: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
712: // TODO: This should be const_iterator, but Sifter sucks
713: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
714: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
715: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
716: const int numPoints = sendOverlap->getConeSize(*r_iter);
717: std::valarray<typename SendOverlap::source_type> sortedPoints(numPoints);
718: int numVals = 0, p = 0;
720: // TODO: This should be const_iterator, but Sifter sucks
721: for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter, ++p) {
722: numVals += sendSection->getFiberDimension(*c_iter);
723: sortedPoints[p] = *c_iter;
724: }
725: typename Section::value_type *v = allocator.allocate(numVals);
726: int k = 0;
728: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
729: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
730: //for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
731: for(p = 0; p < numPoints; ++p) {
732: const typename Section::value_type *vals = sendSection->restrictPoint(sortedPoints[p]);
733: const int dim = sendSection->getFiberDimension(sortedPoints[p]);
735: for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
736: }
737: sendValues[*r_iter] = allocPair(numVals, v);
738: vMover.send(*r_iter, numVals, v);
739: }
740: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
741: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
743: recvSection->allocatePoint();
744: ///recvSection->view("Recv Section after same type copy() allocation");
745: // TODO: This should be const_iterator, but Sifter sucks
746: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
747: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
748: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
749: int numVals = 0;
751: // TODO: This should be const_iterator, but Sifter sucks
752: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
753: numVals += recvSection->getFiberDimension(s_iter.color());
754: }
755: typename Section::value_type *v = allocator.allocate(numVals);
757: recvValues[*r_iter] = allocPair(numVals, v);
758: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
759: vMover.recv(*r_iter, numVals, v);
760: }
761: vMover.start();
762: vMover.end();
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: const int numPoints = recvOverlap->getSupportSize(*r_iter);
767: std::valarray<typename RecvOverlap::color_type> sortedPoints(numPoints);
768: typename Section::value_type *v = recvValues[*r_iter].second;
769: const int numVals = recvValues[*r_iter].first;
770: int k = 0, p = 0;
772: for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter, ++p) {
773: sortedPoints[p] = s_iter.color();
774: }
775: std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
777: //for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
778: for(p = 0; p < numPoints; ++p) {
779: const int size = recvSection->getFiberDimension(sortedPoints[p]);
781: if (size) {recvSection->updatePoint(sortedPoints[p], &v[k]);}
782: k += size;
783: }
784: assert(k == numVals);
785: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
786: allocator.deallocate(v, numVals);
787: }
788: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
789: typename Section::value_type *v = sendValues[*r_iter].second;
790: const int numVals = sendValues[*r_iter].first;
792: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
793: allocator.deallocate(v, numVals);
794: }
795: //recvSection->view("After copy");
796: }
797: // Specialize to ArrowSection
798: template<typename SendOverlap, typename RecvOverlap>
799: 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) {
800: typedef UniformSection<MinimalArrow<int,int>,int> Section;
801: typedef std::pair<int, typename Section::value_type *> allocPair;
802: const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
803: const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
804: MPIMover<typename Section::value_type> vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
805: std::map<int, allocPair> sendValues;
806: std::map<int, allocPair> recvValues;
807: typename Section::alloc_type allocator;
809: copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
810: const Obj<typename SendOverlap::traits::baseSequence> sRanks = sendOverlap->base();
811: const typename SendOverlap::traits::baseSequence::iterator sEnd = sRanks->end();
813: // TODO: This should be const_iterator, but Sifter sucks
814: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
815: const Obj<typename SendOverlap::coneSequence>& points = sendOverlap->cone(*r_iter);
816: const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
817: const typename SendOverlap::coneSequence::iterator pEnd = points->end();
818: int numVals = 0;
820: // TODO: This should be const_iterator, but Sifter sucks
821: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
822: for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
823: numVals += sendSection->getFiberDimension(typename Section::point_type(*c_iter, *d_iter));
824: }
825: }
826: typename Section::value_type *v = allocator.allocate(numVals);
827: int k = 0;
829: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
830: for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
831: for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
832: const typename Section::point_type arrow(*c_iter, *d_iter);
833: const typename Section::value_type *vals = sendSection->restrictPoint(arrow);
834: const int dim = sendSection->getFiberDimension(arrow);
836: for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
837: }
838: }
839: sendValues[*r_iter] = allocPair(numVals, v);
840: vMover.send(*r_iter, numVals, v);
841: }
842: const Obj<typename RecvOverlap::traits::capSequence> rRanks = recvOverlap->cap();
843: const typename RecvOverlap::traits::capSequence::iterator rEnd = rRanks->end();
845: recvSection->allocatePoint();
846: // TODO: This should be const_iterator, but Sifter sucks
847: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
848: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
849: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
850: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
851: int numVals = 0;
853: // TODO: This should be const_iterator, but Sifter sucks
854: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
855: for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
856: numVals += recvSection->getFiberDimension(typename Section::point_type(s_iter.color(), t_iter.color()));
857: }
858: }
859: typename Section::value_type *v = allocator.allocate(numVals);
861: recvValues[*r_iter] = allocPair(numVals, v);
862: for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
863: vMover.recv(*r_iter, numVals, v);
864: }
865: vMover.start();
866: vMover.end();
867: for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
868: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
869: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
870: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
871: typename Section::value_type *v = recvValues[*r_iter].second;
872: const int numVals = recvValues[*r_iter].first;
873: int k = 0;
875: for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
876: for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
877: const typename Section::point_type arrow(s_iter.color(), t_iter.color());
878: const int size = recvSection->getFiberDimension(arrow);
880: if (size) {recvSection->updatePoint(arrow, &v[k]);}
881: k += size;
882: }
883: }
884: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
885: allocator.deallocate(v, numVals);
886: }
887: for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
888: typename Section::value_type *v = sendValues[*r_iter].second;
889: const int numVals = sendValues[*r_iter].first;
891: for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
892: allocator.deallocate(v, numVals);
893: }
894: //recvSection->view("After copy");
895: }
896: // Specialize to a ConstantSection
897: #if 0
898: template<typename SendOverlap, typename RecvOverlap, typename Value>
899: 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) {
900: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
901: };
902: template<typename SendOverlap, typename RecvOverlap, typename Value>
903: 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) {
904: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
905: };
906: #else
907: template<typename SendOverlap, typename RecvOverlap, typename Value>
908: 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) {
909: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
910: }
911: template<typename SendOverlap, typename RecvOverlap, typename Value>
912: 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) {
913: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
914: }
915: #endif
916: // Specialize to an IConstantSection
917: template<typename SendOverlap, typename RecvOverlap, typename Value>
918: 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) {
919: // Why doesn't this work?
920: // This supposed to be a copy, BUT filtered through the sendOverlap
921: // However, an IConstant section does not allow filtration of its
922: // chart. Therefore, you end up with either
923: //
924: // a) Too many items in the chart, copied from the remote sendSection
925: // b) A chart mapped to the local numbering, which we do not want
926: copyIConstant(sendOverlap, recvOverlap, sendSection, recvSection);
927: }
928: // Specialize to an BaseSection/ConstantSection pair
929: #if 0
930: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
931: 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) {
932: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
933: };
934: #else
935: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
936: 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) {
937: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
938: }
939: #endif
940: // Specialize to an BaseSectionV/ConstantSection pair
941: #if 0
942: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
943: 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) {
944: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
945: };
946: #else
947: template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
948: 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) {
949: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
950: }
951: #endif
952: // Specialize to an LabelBaseSection/ConstantSection pair
953: #if 0
954: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
955: 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) {
956: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
957: };
958: #else
959: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
960: 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) {
961: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
962: }
963: #endif
964: template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
965: 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) {
966: copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
967: }
968: // Specialize to a ConstantSection for ArrowSection
969: template<typename SendOverlap, typename RecvOverlap, typename Value>
970: 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) {
971: copyConstantArrow(sendOverlap, recvOverlap, sendSection, recvSection);
972: }
973: };
974: class BinaryFusion {
975: public:
976: template<typename OverlapSection, typename RecvOverlap, typename Section, typename Function>
977: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, Function binaryOp) {
978: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
979: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
981: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
982: // TODO: This must become a more general iterator over colors
983: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
984: // Just taking the first value
985: const typename Section::point_type& localPoint = *p_iter;
986: const typename OverlapSection::point_type& remotePoint = points->begin().color();
987: const typename OverlapSection::value_type *overlapValues = overlapSection->restrictPoint(remotePoint);
988: const typename Section::value_type *localValues = section->restrictPoint(localPoint);
989: const int dim = section->getFiberDimension(localPoint);
990: // TODO: optimize allocation
991: typename Section::value_type *values = new typename Section::value_type[dim];
993: for(int d = 0; d < dim; ++d) {
994: values[d] = binaryOp(overlapValues[d], localValues[d]);
995: }
996: section->updatePoint(localPoint, values);
997: delete [] values;
998: }
999: }
1000: };
1001: class ReplacementBinaryFusion {
1002: public:
1003: template<typename OverlapSection, typename RecvOverlap, typename Section>
1004: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1005: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1006: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1008: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1009: // TODO: This must become a more general iterator over colors
1010: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1011: // Just taking the first value
1012: const typename Section::point_type& localPoint = *p_iter;
1013: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1015: section->update(localPoint, overlapSection->restrictPoint(remotePoint));
1016: }
1017: }
1018: };
1019: class AdditiveBinaryFusion {
1020: public:
1021: template<typename OverlapSection, typename RecvOverlap, typename Section>
1022: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1023: typedef typename Section::point_type point_type;
1024: typedef typename OverlapSection::point_type overlap_point_type;
1025: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
1026: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
1028: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1029: const int rank = *r_iter;
1030: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1031: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1033: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1034: const point_type& localPoint = *p_iter;
1035: const point_type& remotePoint = p_iter.color();
1037: section->updateAddPoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1038: }
1039: }
1040: }
1041: };
1042: class InsertionBinaryFusion {
1043: public:
1044: // Insert the overlapSection values into section along recvOverlap
1045: template<typename OverlapSection, typename RecvOverlap, typename Section>
1046: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1047: typedef typename Section::point_type point_type;
1048: typedef typename OverlapSection::point_type overlap_point_type;
1049: #if 0
1050: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1051: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1053: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1054: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1055: const point_type& localPoint = *p_iter;
1056: const int rank = *points->begin();
1057: const point_type& remotePoint = points->begin().color();
1059: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1060: if (!section->hasPoint(localPoint)) {
1061: 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;
1062: }
1063: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1064: }
1065: }
1066: if (rPoints->size()) {section->allocatePoint();}
1067: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1068: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1069: const point_type& localPoint = *p_iter;
1070: const int rank = *points->begin();
1071: const point_type& remotePoint = points->begin().color();
1073: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1074: section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1075: }
1076: }
1077: #else
1078: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
1079: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
1080: bool hasPoints = false;
1082: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1083: const int rank = *r_iter;
1084: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1085: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1087: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1088: const point_type& localPoint = *p_iter;
1089: const point_type& remotePoint = p_iter.color();
1091: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1092: if (!section->hasPoint(localPoint)) {
1093: 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;
1094: }
1095: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1096: }
1097: hasPoints = true;
1098: }
1099: }
1100: if (hasPoints) {section->allocatePoint();}
1101: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1102: const int rank = *r_iter;
1103: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1104: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1106: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1107: const point_type& localPoint = *p_iter;
1108: const point_type& remotePoint = p_iter.color();
1110: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1111: section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1112: }
1113: }
1114: }
1115: #endif
1116: }
1117: // Specialize to ArrowSection
1118: template<typename OverlapSection, typename RecvOverlap>
1119: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& section) {
1120: typedef UniformSection<MinimalArrow<int,int>,int> Section;
1121: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1122: const typename RecvOverlap::traits::baseSequence::iterator rBegin = rPoints->begin();
1123: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1125: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1126: const Obj<typename RecvOverlap::coneSequence>& sources = recvOverlap->cone(*p_iter);
1127: const typename RecvOverlap::target_type localSource = *p_iter;
1128: const typename RecvOverlap::target_type remoteSource = sources->begin().color();
1130: for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1131: const Obj<typename RecvOverlap::coneSequence>& targets = recvOverlap->cone(*q_iter);
1132: const typename RecvOverlap::target_type localTarget = *q_iter;
1133: const typename RecvOverlap::target_type remoteTarget = targets->begin().color();
1134: const typename Section::point_type localPoint(localSource, localTarget);
1135: const typename OverlapSection::point_type remotePoint(remoteSource, remoteTarget);
1137: if (overlapSection->hasPoint(remotePoint)) {section->setFiberDimension(localPoint, overlapSection->getFiberDimension(remotePoint));}
1138: }
1139: }
1140: if (rPoints->size()) {section->allocatePoint();}
1141: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1142: const Obj<typename RecvOverlap::coneSequence>& sources = recvOverlap->cone(*p_iter);
1143: const typename RecvOverlap::target_type localSource = *p_iter;
1144: const typename RecvOverlap::target_type remoteSource = sources->begin().color();
1146: for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1147: const Obj<typename RecvOverlap::coneSequence>& targets = recvOverlap->cone(*q_iter);
1148: const typename RecvOverlap::target_type localTarget = *q_iter;
1149: const typename RecvOverlap::target_type remoteTarget = targets->begin().color();
1150: const typename Section::point_type localPoint(localSource, localTarget);
1151: const typename OverlapSection::point_type remotePoint(remoteSource, remoteTarget);
1153: if (overlapSection->hasPoint(remotePoint)) {section->updatePoint(localPoint, overlapSection->restrictPoint(remotePoint));}
1154: }
1155: }
1156: }
1157: // Specialize to the Sieve
1158: template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1159: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<Sieve<Point,Point,int> >& sieve) {
1160: typedef typename OverlapSection::point_type overlap_point_type;
1161: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1162: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1164: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1165: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1166: const Point& localPoint = *p_iter;
1167: const int rank = *points->begin();
1168: const Point& remotePoint = points->begin().color();
1169: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1170: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1171: int c = 0;
1173: sieve->clearCone(localPoint);
1174: for(int i = 0; i < size; ++i, ++c) {sieve->addCone(renumbering[values[i]], localPoint, c);}
1175: }
1176: }
1177: // Specialize to the ISieve
1178: template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1179: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<IFSieve<Point> >& sieve) {
1180: typedef typename OverlapSection::point_type overlap_point_type;
1181: #if 0
1182: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1183: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1184: int maxSize = 0;
1186: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1187: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1188: const Point& localPoint = *p_iter;
1189: const int rank = *points->begin();
1190: const Point& remotePoint = points->begin().color();
1191: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1192: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1194: sieve->setConeSize(localPoint, size);
1195: ///for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i]], 1);}
1196: for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
1197: maxSize = std::max(maxSize, size);
1198: }
1199: sieve->allocate();
1200: ///typename OverlapSection::value_type *localValues = new typename OverlapSection::value_type[maxSize];
1201: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1202: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1204: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1205: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1206: const Point& localPoint = *p_iter;
1207: const int rank = *points->begin();
1208: const Point& remotePoint = points->begin().color();
1209: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1210: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1212: ///for(int i = 0; i < size; ++i) {localValues[i] = renumbering[values[i]];}
1213: for(int i = 0; i < size; ++i) {
1214: localValues[i] = renumbering[values[i].first];
1215: localOrientation[i] = values[i].second;
1216: }
1217: sieve->setCone(localValues, localPoint);
1218: sieve->setConeOrientation(localOrientation, localPoint);
1219: }
1220: delete [] localValues;
1221: delete [] localOrientation;
1222: #else
1223: const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
1224: const typename RecvOverlap::capSequence::iterator rEnd = recvOverlap->capEnd();
1225: int maxSize = 0;
1227: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1228: const int rank = *r_iter;
1229: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1230: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1232: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1233: const Point& localPoint = *p_iter;
1234: const Point& remotePoint = p_iter.color();
1235: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1236: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1238: sieve->setConeSize(localPoint, size);
1239: for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
1240: maxSize = std::max(maxSize, size);
1241: }
1242: }
1243: sieve->allocate();
1244: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1245: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1247: for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1248: const int rank = *r_iter;
1249: const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1250: const typename RecvOverlap::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1252: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1253: const Point& localPoint = *p_iter;
1254: const Point& remotePoint = p_iter.color();
1255: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1256: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1258: for(int i = 0; i < size; ++i) {
1259: localValues[i] = renumbering[values[i].first];
1260: localOrientation[i] = values[i].second;
1261: }
1262: sieve->setCone(localValues, localPoint);
1263: sieve->setConeOrientation(localOrientation, localPoint);
1264: }
1265: }
1266: delete [] localValues;
1267: delete [] localOrientation;
1268: #endif
1269: }
1270: template<typename OverlapSection, typename RecvOverlap, typename Point>
1271: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<IFSieve<Point> >& sieve) {
1272: typedef typename OverlapSection::point_type overlap_point_type;
1273: #if 0
1274: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1275: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1276: int maxSize = 0;
1278: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1279: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1280: const Point& localPoint = *p_iter;
1281: const int rank = *points->begin();
1282: const Point& remotePoint = points->begin().color();
1283: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1284: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1286: sieve->setConeSize(localPoint, size);
1287: for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
1288: maxSize = std::max(maxSize, size);
1289: }
1290: sieve->allocate();
1291: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1292: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1294: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1295: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1296: const Point& localPoint = *p_iter;
1297: const int rank = *points->begin();
1298: const Point& remotePoint = points->begin().color();
1299: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1300: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1302: for(int i = 0; i < size; ++i) {
1303: localValues[i] = values[i].first;
1304: localOrientation[i] = values[i].second;
1305: }
1306: sieve->setCone(localValues, localPoint);
1307: sieve->setConeOrientation(localOrientation, localPoint);
1308: }
1309: delete [] localValues;
1310: delete [] localOrientation;
1311: #else
1312: const Obj<typename RecvOverlap::capSequence> rRanks = recvOverlap->cap();
1313: const typename RecvOverlap::capSequence::iterator rEnd = rRanks->end();
1314: int maxSize = 0;
1316: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
1317: const int rank = *r_iter;
1318: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
1319: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
1320: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
1322: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1323: const Point& localPoint = *p_iter;
1324: const Point& remotePoint = p_iter.color();
1325: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1326: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1328: sieve->setConeSize(localPoint, size);
1329: for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
1330: maxSize = std::max(maxSize, size);
1331: }
1332: }
1333: sieve->allocate();
1334: typename OverlapSection::value_type::first_type *localValues = new typename OverlapSection::value_type::first_type[maxSize];
1335: typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];
1337: for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
1338: const int rank = *r_iter;
1339: const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter);
1340: const typename RecvOverlap::supportSequence::iterator pBegin = points->begin();
1341: const typename RecvOverlap::supportSequence::iterator pEnd = points->end();
1343: for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1344: const Point& localPoint = *p_iter;
1345: const Point& remotePoint = p_iter.color();
1346: const int size = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1347: const typename OverlapSection::value_type *values = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1349: for(int i = 0; i < size; ++i) {
1350: localValues[i] = values[i].first;
1351: localOrientation[i] = values[i].second;
1352: }
1353: sieve->setCone(localValues, localPoint);
1354: sieve->setConeOrientation(localOrientation, localPoint);
1355: }
1356: }
1357: delete [] localValues;
1358: delete [] localOrientation;
1359: #endif
1360: }
1361: // Generic
1362: template<typename OverlapSection, typename RecvOverlap, typename Section, typename Bundle>
1363: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, const Obj<Bundle>& bundle) {
1364: typedef typename OverlapSection::point_type overlap_point_type;
1365: const Obj<typename RecvOverlap::baseSequence> rPoints = recvOverlap->base();
1366: const typename RecvOverlap::baseSequence::iterator rEnd = rPoints->end();
1368: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1369: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1370: const typename Section::point_type& localPoint = *p_iter;
1371: const int rank = *points->begin();
1372: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1374: section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1375: }
1376: bundle->allocate(section);
1377: for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1378: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1379: const typename Section::point_type& localPoint = *p_iter;
1380: const int rank = *points->begin();
1381: const typename OverlapSection::point_type& remotePoint = points->begin().color();
1383: section->update(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1384: }
1385: }
1386: };
1387: class InterpolateMultipleFusion {
1388: public:
1389: // Interpolate the overlapSection values into section along recvOverlap
1390: template<typename OverlapSection, typename RecvOverlap, typename Section>
1391: static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1392: typedef typename Section::point_type point_type;
1393: typedef typename Section::value_type value_type;
1394: typedef typename OverlapSection::point_type overlap_point_type;
1395: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
1396: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
1397: int maxFiberDim = -1;
1399: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1400: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1401: const typename RecvOverlap::coneSequence::iterator rpEnd = points->end();
1402: const point_type& localPoint = *p_iter;
1403: bool inOverlap = false;
1404: int fiberDim = -1;
1406: for(typename RecvOverlap::coneSequence::iterator rp_iter = points->begin(); rp_iter != rpEnd; ++rp_iter) {
1407: const int rank = *rp_iter;
1408: const point_type& remotePoint = rp_iter.color();
1410: if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1411: inOverlap = true;
1412: fiberDim = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1413: break;
1414: }
1415: }
1416: if (inOverlap) {
1417: if (!section->hasPoint(localPoint)) {
1418: std::cout <<"["<<section->commRank()<<"]: Destination section does not have local point " << localPoint << " remote point " << (points->begin().color()) << " fiber dim " << fiberDim << std::endl;
1419: }
1420: section->setFiberDimension(localPoint, fiberDim);
1421: maxFiberDim = std::max(fiberDim, maxFiberDim);
1422: }
1423: }
1424: if (rPoints->size()) {section->allocatePoint();}
1425: value_type *interpolant = new value_type[maxFiberDim];
1427: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1428: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1429: const typename RecvOverlap::coneSequence::iterator rpEnd = points->end();
1430: const point_type& localPoint = *p_iter;
1431: bool inOverlap = false;
1432: int numArgs = 0;
1434: for(int d = 0; d < maxFiberDim; ++d) {interpolant[d] = 0.0;}
1435: for(typename RecvOverlap::coneSequence::iterator rp_iter = points->begin(); rp_iter != rpEnd; ++rp_iter) {
1436: const int rank = *rp_iter;
1437: const point_type& remotePoint = rp_iter.color();
1438: const overlap_point_type opoint(rank, remotePoint);
1440: if (overlapSection->hasPoint(opoint)) {
1441: const int fiberDim = overlapSection->getFiberDimension(opoint);
1442: const value_type *values = overlapSection->restrictPoint(opoint);
1444: // TODO: Include interpolation weights (stored in overlap)
1445: for(int d = 0; d < fiberDim; ++d) {
1446: interpolant[d] += values[d];
1447: }
1448: inOverlap = true;
1449: ++numArgs;
1450: }
1451: }
1452: if (inOverlap) {
1453: for(int d = 0; d < maxFiberDim; ++d) {interpolant[d] /= numArgs;}
1454: section->updatePoint(localPoint, interpolant);
1455: }
1456: }
1457: delete [] interpolant;
1458: }
1459: };
1460: }
1461: }
1463: #endif