Actual source code: Completion.hh
1: #ifndef included_ALE_Completion_hh
2: #define included_ALE_Completion_hh
4: #ifndef included_ALE_Sections_hh
5: #include <Sections.hh>
6: #endif
8: #ifndef included_ALE_ParallelMapping_hh
9: #include <ParallelMapping.hh>
10: #endif
12: #include <iostream>
13: #include <fstream>
15: namespace ALE {
16: class Completion {
17: public:
18: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
19: static void completeSection(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
20: typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
21: //typedef ALE::Section<typename SendSection::point_type, typename SendSection::value_type> OverlapSection;
22: Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());
24: if (sendSection->debug()) {sendSection->view("Send Section");}
25: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
26: if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
27: ALE::Pullback::InsertionBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
28: if (recvSection->debug()) {recvSection->view("Receieve Section");}
29: };
30: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
31: static void completeSectionAdd(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
32: typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
33: Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());
35: if (sendSection->debug()) {sendSection->view("Send Section");}
36: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
37: if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
38: ALE::Pullback::AdditiveBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
39: if (recvSection->debug()) {recvSection->view("Receieve Section");}
40: };
41: };
42: namespace New {
43: template<typename Bundle_, typename Value_, typename Alloc_ = malloc_allocator<typename Bundle_::point_type> >
44: class Completion {
45: public:
46: typedef int point_type;
47: typedef Value_ value_type;
48: typedef Bundle_ bundle_type;
49: typedef Alloc_ alloc_type;
50: typedef typename alloc_type::template rebind<int>::other int_alloc_type;
51: typedef typename alloc_type::template rebind<value_type>::other value_alloc_type;
52: typedef typename bundle_type::sieve_type sieve_type;
53: typedef typename ALE::DiscreteSieve<point_type, alloc_type> dsieve_type;
54: typedef typename ALE::Topology<int, dsieve_type, alloc_type> topology_type;
55: typedef typename bundle_type::rank_type rank_type;
56: typedef typename ALE::Sifter<point_type, rank_type, point_type> send_overlap_type;
57: typedef typename ALE::Sifter<rank_type, point_type, point_type> recv_overlap_type;
58: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > send_sizer_type;
59: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > recv_sizer_type;
60: typedef typename ALE::New::ConeSizeSection<bundle_type, sieve_type> cone_size_section;
61: typedef typename ALE::New::ConeSection<sieve_type> cone_section;
62: typedef typename ALE::New::SectionCompletion<bundle_type, value_type, alloc_type> completion;
63: public:
64: template<typename PartitionType>
65: static void scatterSieve(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, const int dim, const Obj<sieve_type>& sieveNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height, const int numCells, const PartitionType assignment[]) {
66: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
67: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
68: int rank = sieve->commRank();
69: int debug = sieve->debug();
71: // Create local sieve
72: const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(height);
73: int e = 0;
75: if (sieve->debug()) {
76: int e2 = 0;
77: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
78: std::cout << "assignment["<<*e_iter<<"]" << assignment[e2++] << std::endl;
79: }
80: }
81: PetscTruth flg;
82: PetscOptionsHasName(PETSC_NULL, "-output_partition", &flg);
83: if (flg) {
84: ostringstream fname;
85: fname << "part." << sieve->commSize() << ".dat";
86: std::ofstream f(fname.str().c_str());
87: int e2 = 0;
88: f << sieve->commSize() << std::endl;
89: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
90: f << assignment[e2++] << std::endl;
91: }
92: }
93: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
94: if (assignment[e] == rank) {
95: Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
96: Obj<typename sieve_type::coneSet> next = new typename sieve_type::coneSet();
97: Obj<typename sieve_type::coneSet> tmp;
99: current->insert(*e_iter);
100: while(current->size()) {
101: for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
102: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
103:
104: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
105: sieveNew->addArrow(*c_iter, *p_iter, c_iter.color());
106: next->insert(*c_iter);
107: }
108: }
109: tmp = current; current = next; next = tmp;
110: next->clear();
111: }
112: if (height) {
113: current->insert(*e_iter);
114: while(current->size()) {
115: for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
116: const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
117:
118: for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
119: sieveNew->addArrow(*p_iter, *s_iter, s_iter.color());
120: next->insert(*s_iter);
121: }
122: }
123: tmp = current; current = next; next = tmp;
124: next->clear();
125: }
126: }
127: }
128: e++;
129: }
130: // Complete sizer section
131: typedef typename ALE::New::PartitionSizeSection<bundle_type, PartitionType> partition_size_section;
132: typedef typename ALE::New::PartitionSection<bundle_type, PartitionType> partition_section;
133: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
134: Obj<partition_size_section> partitionSizeSection = new partition_size_section(bundle, height, numCells, assignment);
135: Obj<partition_section> partitionSection = new partition_section(bundle, height, numCells, assignment);
136: Obj<send_section_type> sendSection = new send_section_type(sieve->comm(), sieve->debug());
137: Obj<recv_section_type> recvSection = new recv_section_type(sieve->comm(), sendSection->getTag(), sieve->debug());
139: completion::completeSection(sendOverlap, recvOverlap, partitionSizeSection, partitionSection, sendSection, recvSection);
140: // Unpack the section into the overlap
141: sendOverlap->clear();
142: recvOverlap->clear();
143: const typename send_section_type::sheaf_type& sendPatches = sendSection->getPatches();
145: for(typename send_section_type::sheaf_type::const_iterator p_iter = sendPatches.begin(); p_iter != sendPatches.end(); ++p_iter) {
146: const typename send_section_type::patch_type rank = p_iter->first;
147: const Obj<typename send_section_type::section_type>& section = p_iter->second;
148: const typename send_section_type::section_type::chart_type chart = section->getChart();
150: for(typename send_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
151: const typename send_section_type::value_type *points = section->restrictPoint(*c_iter);
152: int size = section->getFiberDimension(*c_iter);
154: for(int p = 0; p < size; p++) {
155: sendOverlap->addArrow(points[p], rank, points[p]);
156: }
157: }
158: }
159: const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();
161: for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
162: const typename send_section_type::patch_type rank = p_iter->first;
163: const Obj<typename send_section_type::section_type>& section = p_iter->second;
164: const typename send_section_type::section_type::chart_type chart = section->getChart();
166: for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
167: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
168: int size = section->getFiberDimension(*c_iter);
170: for(int p = 0; p < size; p++) {
171: recvOverlap->addArrow(rank, points[p], points[p]);
172: }
173: }
174: }
175: if (debug) {
176: sendOverlap->view(std::cout, "Send overlap for points");
177: recvOverlap->view(std::cout, "Receive overlap for points");
178: }
179: // Receive the point section
180: ALE::New::Completion<bundle_type, value_type>::scatterCones(sieve, sieveNew, sendOverlap, recvOverlap, bundle, height);
181: if (height) {
182: ALE::New::Completion<bundle_type, value_type>::scatterSupports(sieve, sieveNew, sendOverlap, recvOverlap, bundle, bundle->depth()-height);
183: }
184: };
185: template<typename SifterType>
186: static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumHeight = 0) {
187: typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
188: typedef typename ALE::New::ConeSection<SifterType> cone_section;
189: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
190: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
191: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
192: Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter, minimumHeight);
193: Obj<cone_section> coneSection = new cone_section(sifter);
194: Obj<send_section_type> sendSection = new send_section_type(sifter->comm(), sifter->debug());
195: Obj<recv_section_type> recvSection = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());
197: completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
198: // Unpack the section into the sieve
199: const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();
201: for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
202: const Obj<typename recv_section_type::section_type>& section = p_iter->second;
203: const typename recv_section_type::section_type::chart_type& chart = section->getChart();
205: for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
206: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
207: int size = section->getFiberDimension(*c_iter);
208: int c = 0;
210: for(int p = 0; p < size; p++) {
211: sifterNew->addArrow(points[p], *c_iter, c++);
212: }
213: }
214: }
215: };
216: template<typename SifterType, typename Renumbering>
217: static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, Renumbering& renumbering, const Obj<bundle_type>& bundle = NULL) {
218: PETSc::Log::Event("ScatterCones").begin();
219: typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
220: typedef typename ALE::New::ConeSection<SifterType> cone_section;
221: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
222: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
223: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
224: Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter);
225: Obj<cone_section> coneSection = new cone_section(sifter);
226: Obj<send_section_type> sendSection = new send_section_type(sifter->comm(), sifter->debug());
227: Obj<recv_section_type> recvSection = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());
229: PETSc::Log::Event("ScatterConesComplete").begin();
230: completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
231: PETSc::Log::Event("ScatterConesComplete").end();
232: PETSc::Log::Event("ScatterConesUpdate").begin();
233: // Unpack the section into the sieve
234: const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();
236: for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
237: const Obj<typename recv_section_type::section_type>& section = p_iter->second;
238: const typename recv_section_type::section_type::chart_type& chart = section->getChart();
240: for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
241: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
242: int size = section->getFiberDimension(*c_iter);
243: int c = 0;
245: for(int p = 0; p < size; p++) {
246: sifterNew->addArrow(points[p], renumbering[*c_iter], c++);
247: }
248: }
249: }
250: PETSc::Log::Event("ScatterConesUpdate").end();
251: PETSc::Log::Event("ScatterCones").end();
252: };
253: template<typename SifterType>
254: static void scatterSupports(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumDepth = 0) {
255: typedef typename ALE::New::SupportSizeSection<bundle_type, SifterType> support_size_section;
256: typedef typename ALE::New::SupportSection<SifterType> support_section;
257: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
258: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
259: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
260: Obj<support_size_section> supportSizeSection = new support_size_section(bundle, sifter, minimumDepth);
261: Obj<support_section> supportSection = new support_section(sifter);
262: Obj<send_section_type> sendSection = new send_section_type(sifter->comm(), sifter->debug());
263: Obj<recv_section_type> recvSection = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());
265: completion::completeSection(sendOverlap, recvOverlap, supportSizeSection, supportSection, sendSection, recvSection);
266: // Unpack the section into the sieve
267: const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();
269: for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
270: const Obj<typename send_section_type::section_type>& section = p_iter->second;
271: const typename send_section_type::section_type::chart_type chart = section->getChart();
273: for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
274: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
275: int size = section->getFiberDimension(*c_iter);
276: int c = 0;
278: for(int p = 0; p < size; p++) {
279: sifterNew->addArrow(*c_iter, points[p], c++);
280: }
281: }
282: }
283: };
284: };
286: template<typename Value_>
287: class ParallelFactory {
288: public:
289: typedef Value_ value_type;
290: protected:
291: int _debug;
292: MPI_Datatype _mpiType;
293: protected:
294: MPI_Datatype constructMPIType() {
295: if (sizeof(value_type) == 1) {
296: return MPI_BYTE;
297: } else if (sizeof(value_type) == 2) {
298: return MPI_SHORT;
299: } else if (sizeof(value_type) == 4) {
300: return MPI_INT;
301: } else if (sizeof(value_type) == 8) {
302: return MPI_DOUBLE;
303: } else if (sizeof(value_type) == 28) {
304: int blen[2];
305: MPI_Aint indices[2];
306: MPI_Datatype oldtypes[2], newtype;
307: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
308: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
309: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
310: MPI_Type_commit(&newtype);
311: return newtype;
312: } else if (sizeof(value_type) == 32) {
313: int blen[2];
314: MPI_Aint indices[2];
315: MPI_Datatype oldtypes[2], newtype;
316: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_DOUBLE;
317: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
318: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
319: MPI_Type_commit(&newtype);
320: return newtype;
321: }
322: throw ALE::Exception("Cannot determine MPI type for value type");
323: };
324: ParallelFactory(const int debug) : _debug(debug) {
325: this->_mpiType = this->constructMPIType();
326: };
327: public:
328: ~ParallelFactory() {};
329: public:
330: static const Obj<ParallelFactory>& singleton(const int debug, bool cleanup = false) {
331: static Obj<ParallelFactory> *_singleton = NULL;
333: if (cleanup) {
334: if (debug) {std::cout << "Destroying ParallelFactory" << std::endl;}
335: if (_singleton) {delete _singleton;}
336: _singleton = NULL;
337: } else if (_singleton == NULL) {
338: if (debug) {std::cout << "Creating new ParallelFactory" << std::endl;}
339: _singleton = new Obj<ParallelFactory>();
340: *_singleton = new ParallelFactory(debug);
341: }
342: return *_singleton;
343: };
344: public: // Accessors
345: int debug() const {return this->_debug;};
346: MPI_Datatype getMPIType() const {return this->_mpiType;};
347: };
348: }
349: }
350: #endif