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