1/*
2 * MetisDistribution.hpp
3 *
4 * Created on: Nov 19, 2015
5 * Author: Antonio Leo
6 */
7
8#ifndef SRC_DECOMPOSITION_BOXDISTRIBUTION_HPP_
9#define SRC_DECOMPOSITION_BOXDISTRIBUTION_HPP_
10
11#include "SubdomainGraphNodes.hpp"
12#include "Vector/map_vector.hpp"
13#include "VCluster/VCluster.hpp"
14#include "Graph/map_graph.hpp"
15#include "Graph/CartesianGraphFactory.hpp"
16#include "VTKWriter/VTKWriter.hpp"
17
18#define BOX_DISTRIBUTION_ERROR_OBJECT std::runtime_error("Box distribution runtime error");
19
20
21void getPrimeFactors(int n, openfpm::vector<int> & fact )
22{
23 while (n%2 == 0)
24 {
25 fact.add(2);
26 n = n/2;
27 }
28 for (int i = 3; i <= sqrt(n); i = i+2)
29 {
30 while (n%i == 0)
31 {
32 fact.add(i);
33 n = n/i;
34 }
35 }
36 if (n > 2)
37 {fact.add(n);}
38}
39
40
41/*! \brief Class that distribute sub-sub-domains across processors using Metis Library
42 *
43 * Given a graph and setting Computational cost, Communication cost (on the edge) and
44 * Migration cost or total Communication costs, it produce the optimal distribution
45 *
46 * ### Initialize a Cartesian graph and decompose
47 * \snippet Distribution_unit_tests.hpp Initialize a Metis Cartesian graph and decompose
48 *
49 * ### Set Computation Communication and Migration cost
50 * \snippet Distribution_unit_tests.hpp Decomposition Metis with weights
51 *
52 */
53
54template<unsigned int dim, typename T>
55class BoxDistribution
56{
57 //! Vcluster
58 Vcluster<> & v_cl;
59
60 //! Structure that store the cartesian grid information
61 grid_sm<dim, void> gr;
62
63 //! rectangular domain to decompose
64 Box<dim, T> domain;
65
66 //! Global sub-sub-domain graph
67 Graph_CSR<nm_v<dim>, nm_e> gp;
68
69 //! Flag that indicate if we are doing a test (In general it fix the seed)
70 bool testing = false;
71
72 openfpm::vector<int> subsub_own;
73
74 /*! \brief Check that the sub-sub-domain id exist
75 *
76 * \param id sub-sub-domain id
77 *
78 */
79 inline void check_overflow(size_t id)
80 {
81#ifdef SE_CLASS1
82 if (id >= gp.getNVertex())
83 {
84 std::cerr << "Error " << __FILE__ ":" << __LINE__ << " such sub-sub-domain doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
85 ACTION_ON_ERROR(BOX_DISTRIBUTION_ERROR_OBJECT)
86 }
87#endif
88 }
89
90 /*! \brief Check that the sub-sub-domain id exist
91 *
92 * \param id sub-sub-domain id
93 *
94 */
95 inline void check_overflowe(size_t id, size_t e)
96 {
97#ifdef SE_CLASS1
98 if (e >= gp.getNChilds(id))
99 {
100 std::cerr << "Error " << __FILE__ ":" << __LINE__ << " for the sub-sub-domain " << id << " such neighborhood doesn't exist (e = " << e << ", " << "total size = " << gp.getNChilds(id) << ")\n";
101 ACTION_ON_ERROR(BOX_DISTRIBUTION_ERROR_OBJECT)
102 }
103#endif
104 }
105
106public:
107
108 static constexpr unsigned int computation = nm_v_computation;
109
110 /*! \brief constructor
111 *
112 * \param v_cl vcluster
113 *
114 */
115 BoxDistribution(Vcluster<> & v_cl)
116 :v_cl(v_cl)
117 {
118#ifdef SE_CLASS2
119 check_new(this,8,VECTOR_EVENT,1);
120#endif
121 }
122
123 /*! \brief Copy constructor
124 *
125 * \param mt distribution to copy
126 *
127 */
128 BoxDistribution(const BoxDistribution & mt)
129 :v_cl(mt.v_cl)
130 {
131#ifdef SE_CLASS2
132 check_valid(mt);
133 check_new(this,8,VECTOR_EVENT,1);
134#endif
135 this->operator=(mt);
136 }
137
138 /*! \brief Copy constructor
139 *
140 * \param mt distribution to copy
141 *
142 */
143 BoxDistribution(BoxDistribution && mt)
144 {
145#ifdef SE_CLASS2
146 check_valid(mt);
147 check_new(this,8,VECTOR_EVENT,1);
148#endif
149 this->operator=(mt);
150 }
151
152 /*! \brief Destructor
153 *
154 *
155 */
156 ~BoxDistribution()
157 {
158#ifdef SE_CLASS2
159 check_delete(this);
160#endif
161 }
162
163
164 /*! \brief create a Cartesian distribution graph
165 *
166 * \param grid grid info (sub-sub somains on each dimension)
167 * \param dom domain (domain where the sub-sub-domains are defined)
168 *
169 */
170 void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
171 {
172#ifdef SE_CLASS2
173 check_valid(this,8);
174#endif
175 // NON periodic boundary conditions
176 size_t bc[dim];
177
178 for (size_t i = 0 ; i < dim ; i++)
179 bc[i] = NON_PERIODIC;
180
181 // Set grid and domain
182 gr = grid;
183 domain = dom;
184
185 // Create a cartesian grid graph
186 CartesianGraphFactory<dim, Graph_CSR<nm_v<dim>, nm_e>> g_factory_part;
187 gp = g_factory_part.template construct<NO_EDGE, nm_v_id, T, dim - 1, 0>(gr.getSize(), domain, bc);
188
189 // Init to 0.0 axis z (to fix in graphFactory)
190 if (dim < 3)
191 {
192 for (size_t i = 0; i < gp.getNVertex(); i++)
193 {
194 gp.vertex(i).template get<nm_v_x>()[2] = 0.0;
195 }
196 }
197
198 for (size_t i = 0; i < gp.getNVertex(); i++)
199 gp.vertex(i).template get<nm_v_global_id>() = i;
200 }
201
202 /*! \brief Get the current graph (main)
203 *
204 * \return the current sub-sub domain Graph
205 *
206 */
207 Graph_CSR<nm_v<dim>, nm_e> & getGraph()
208 {
209#ifdef SE_CLASS2
210 check_valid(this,8);
211#endif
212 return gp;
213 }
214
215 /*! \brief Distribute the sub-sub-domains
216 *
217 */
218 void decompose()
219 {
220#ifdef SE_CLASS2
221 check_valid(this,8);
222#endif
223
224 subsub_own.clear();
225 grid_key_dx_iterator<dim> it(gr);
226
227 openfpm::vector<int> facts;
228 getPrimeFactors(v_cl.size(),facts);
229
230 size_t div[dim];
231 size_t ln[dim];
232 double ln_d[dim];
233
234 for (int i = 0 ; i < dim ; i++)
235 {div[i] = 1;}
236
237 for (int i = 0 ; i < facts.size() ; i++)
238 {div[i % dim] *= facts.get(i);}
239
240 for (int i = 0 ; i < dim ; i++)
241 {
242 ln[i] = gr.size(i) / div[i];
243 ln_d[i] = (double)gr.size(i) / div[i];
244 }
245
246 grid_sm<dim,void> gr_proc(div);
247
248 while (it.isNext())
249 {
250 auto key = it.get();
251
252 grid_key_dx<dim> key_prc;
253
254 for (int i = 0 ; i < dim ; i++)
255 {
256 key_prc.set_d(i,key.get(i)/ln_d[i]);
257 if (key_prc.get(i) >= div[i])
258 {key_prc.set_d(i,div[i]-1);}
259 }
260
261 size_t i = gr.LinId(key);
262
263 gp.vertex(i).template get<nm_v_proc_id>() = gr_proc.LinId(key_prc);
264
265 if (gp.vertex(i).template get<nm_v_proc_id>() == v_cl.rank())
266 {subsub_own.add(i);}
267
268 ++it;
269 }
270 }
271
272 /*! \brief Refine current decomposition
273 *
274 * In metis case it just re-decompose
275 *
276 */
277 void refine()
278 {
279#ifdef SE_CLASS2
280 check_valid(this,8);
281#endif
282
283 decompose();
284 }
285
286 /*! \brief Redecompose current decomposition
287 *
288 */
289 void redecompose()
290 {
291#ifdef SE_CLASS2
292 check_valid(this,8);
293#endif
294 decompose();
295 }
296
297
298 /*! \brief Function that return the position (point P1) of the sub-sub domain box in the space
299 *
300 * \param id vertex id
301 * \param pos vector that contain x, y, z
302 *
303 */
304 void getSSDomainPos(size_t id, T (&pos)[dim])
305 {
306#ifdef SE_CLASS2
307 check_valid(this,8);
308#endif
309 check_overflow(id);
310
311 // Copy the geometrical informations inside the pos vector
312 pos[0] = gp.vertex(id).template get<nm_v_x>()[0];
313 pos[1] = gp.vertex(id).template get<nm_v_x>()[1];
314 if (dim == 3)
315 {pos[2] = gp.vertex(id).template get<nm_v_x>()[2];}
316 }
317
318 /*! \brief function that get the computational cost of the sub-sub-domain
319 *
320 * \param id sub-sub-domain
321 *
322 * \return the comutational cost
323 *
324 */
325 size_t getComputationalCost(size_t id)
326 {
327#ifdef SE_CLASS2
328 check_valid(this,8);
329#endif
330 check_overflow(id);
331 return gp.vertex(id).template get<nm_v_computation>();
332 }
333
334
335 /*! \brief Set computation cost on a sub-sub domain
336 *
337 * \param id sub-sub domain id
338 * \param cost
339 *
340 */
341 void setComputationCost(size_t id, size_t cost)
342 {
343#ifdef SE_CLASS2
344 check_valid(this,8);
345#endif
346#ifdef SE_CLASS1
347 check_overflow(id);
348#endif
349 }
350
351 /*! \brief Set migration cost on a sub-sub domain
352 *
353 * \param id of the sub-sub domain
354 * \param cost
355 */
356 void setMigrationCost(size_t id, size_t cost)
357 {
358#ifdef SE_CLASS2
359 check_valid(this,8);
360#endif
361#ifdef SE_CLASS1
362 check_overflow(id);
363#endif
364 }
365
366 /*! \brief Set communication cost between neighborhood sub-sub-domains (weight on the edge)
367 *
368 * \param id sub-sub domain
369 * \param e id in the neighborhood list (id in the adjacency list)
370 * \param cost
371 */
372 void setCommunicationCost(size_t id, size_t e, size_t cost)
373 {
374#ifdef SE_CLASS2
375 check_valid(this,8);
376#endif
377#ifdef SE_CLASS1
378 check_overflow(id);
379 check_overflowe(id,e);
380#endif
381 }
382
383 /*! \brief Returns total number of sub-sub-domains
384 *
385 * \return sub-sub domain numbers
386 *
387 */
388 size_t getNSubSubDomains()
389 {
390#ifdef SE_CLASS2
391 check_valid(this,8);
392#endif
393 return gp.getNVertex();
394 }
395
396 /*! \brief Returns total number of neighbors of one sub-sub-domain
397 *
398 * \param id of the sub-sub-domain
399 */
400 size_t getNSubSubDomainNeighbors(size_t id)
401 {
402#ifdef SE_CLASS2
403 check_valid(this,8);
404#endif
405#ifdef SE_CLASS1
406 check_overflow(id);
407#endif
408
409 return gp.getNChilds(id);
410 }
411
412 /*! \brief Compute the unbalance of the processor compared to the optimal balance
413 *
414 * \warning all processor must call this function
415 *
416 * \return the unbalance from the optimal one 0.01 mean 1%
417 */
418 float getUnbalance()
419 {
420#ifdef SE_CLASS2
421 check_valid(this,8);
422#endif
423 return 0.0;
424 }
425
426 /*! \brief Return the total number of sub-sub-domains in the distribution graph
427 *
428 * \return the total number of sub-sub-domains set
429 *
430 */
431 size_t getNOwnerSubSubDomains() const
432 {
433 return subsub_own.size();
434 }
435
436 /*! \brief Return the id of the set sub-sub-domain
437 *
438 * \param id id in the list of the set sub-sub-domains
439 *
440 * \return the id
441 *
442 */
443 size_t getOwnerSubSubDomain(size_t id) const
444 {
445 return subsub_own.get(id);
446 }
447
448 /*! \brief It set the Classs on test mode
449 *
450 * At the moment it fix the seed to have reproducible results
451 *
452 */
453 void onTest()
454 {
455#ifdef SE_CLASS2
456 check_valid(this,8);
457#endif
458 testing = true;
459 }
460
461 /*! \brief Write the distribution graph into file
462 *
463 * \param out output filename
464 *
465 */
466 void write(std::string out)
467 {
468#ifdef SE_CLASS2
469 check_valid(this,8);
470#endif
471
472 VTKWriter<Graph_CSR<nm_v<dim>, nm_e>, VTK_GRAPH> gv2(gp);
473 gv2.write(std::to_string(v_cl.getProcessUnitID()) + "_" + out + ".vtk");
474
475 }
476
477 /*! \brief Compute the processor load
478 *
479 * \warning all processors must call this function
480 *
481 * \return the total computation cost
482 */
483 size_t getProcessorLoad()
484 {
485#ifdef SE_CLASS2
486 check_valid(this,8);
487#endif
488 openfpm::vector<size_t> loads(v_cl.getProcessingUnits());
489
490 size_t load = 0;
491
492 if (v_cl.getProcessUnitID() == 0)
493 {
494 for (size_t i = 0; i < gp.getNVertex(); i++)
495 {loads.get(gp.template vertex_p<nm_v_proc_id>(i)) += gp.template vertex_p<nm_v_computation>(i);}
496
497 for (size_t i = 0 ; i < v_cl.getProcessingUnits() ; i++)
498 {
499 v_cl.send(i,1234,&loads.get(i),sizeof(size_t));
500 }
501 }
502 v_cl.recv(0,1234,&load,sizeof(size_t));
503 v_cl.execute();
504
505 return load;
506 }
507
508 /*! \brief operator=
509 *
510 * \param mt object to copy
511 *
512 * \return itself
513 *
514 */
515 BoxDistribution & operator=(const BoxDistribution & mt)
516 {
517#ifdef SE_CLASS2
518 check_valid(&mt,8);
519 check_valid(this,8);
520#endif
521
522 this->gr = mt.gr;
523 this->domain = mt.domain;
524 this->gp = mt.gp;
525 this->subsub_own = mt.subsub_own;
526 return *this;
527 }
528
529 /*! \brief operator=
530 *
531 * \param mt object to copy
532 *
533 * \return itself
534 *
535 */
536 BoxDistribution & operator=(BoxDistribution && mt)
537 {
538#ifdef SE_CLASS2
539 check_valid(mt);
540 check_valid(this,8);
541#endif
542 this->gr = mt.gr;
543 this->domain = mt.domain;
544 this->gp.swap(mt.gp);
545 this->subsub_own.swap(mt.subsub_own);
546 return *this;
547 }
548
549 /*! \brief operator==
550 *
551 * \param mt Metis distribution to compare with
552 *
553 * \return true if the distribution match
554 *
555 */
556 inline bool operator==(const BoxDistribution & mt)
557 {
558#ifdef SE_CLASS2
559 check_valid(&mt,8);
560 check_valid(this,8);
561#endif
562 bool ret = true;
563
564 ret &= (this->gr == mt.gr);
565 ret &= (this->domain == mt.domain);
566 ret &= (this->gp == mt.gp);
567
568 ret &= this->subsub_own == mt.subsub_own;
569
570 return ret;
571 }
572
573 /*! \brief Set the tolerance for each partition
574 *
575 * \param tol tolerance
576 *
577 */
578 void setDistTol(double tol)
579 {
580 }
581
582
583
584 /*! \brief function that get the weight of the vertex
585 *
586 * \param id vertex id
587 *
588 */
589 size_t getSubSubDomainComputationCost(size_t id)
590 {
591#ifdef SE_CLASS1
592 if (id >= gp.getNVertex())
593 std::cerr << __FILE__ << ":" << __LINE__ << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
594#endif
595 return 1;
596 }
597
598 /*! \brief Get the decomposition counter
599 *
600 * \return the decomposition counter
601 *
602 */
603 size_t get_ndec()
604 {
605 return 0;
606 }
607};
608
609#endif /* SRC_DECOMPOSITION_METISDISTRIBUTION_HPP_ */
610