1/*
2 * parmetis_util.hpp
3 *
4 * Created on: Oct 07, 2015
5 * Author: Antonio Leo
6 */
7
8#ifndef PARMETIS_UTIL_HPP
9#define PARMETIS_UTIL_HPP
10
11#include <iostream>
12#include "parmetis.h"
13#include "VTKWriter/VTKWriter.hpp"
14#include "VCluster/VCluster.hpp"
15#include "Graph/ids.hpp"
16
17#define PARMETIS_ERROR_OBJECT std::runtime_error("Runtime Parmetis error");
18
19
20/*! \brief Metis graph structure
21 *
22 * Metis graph structure
23 *
24 */
25struct Parmetis_graph
26{
27 //! The number of vertices in the graph
28 idx_t * nvtxs;
29
30 //! number of balancing constrains
31 //! more practical, are the number of weights for each vertex
32 //! PS even we you specify vwgt == NULL ncon must be set at leat to
33 //! one
34 idx_t * ncon;
35
36 //! For each vertex it store the adjacency lost start for the vertex i
37 idx_t * xadj;
38
39 //! For each vertex it store a list of all neighborhood vertex
40 idx_t * adjncy;
41
42 //! Array that store the weight for each vertex
43 idx_t * vwgt;
44
45 //! Array of the vertex size, basically is the total communication amount
46 idx_t * vsize;
47
48 //! The weight of the edge
49 idx_t * adjwgt;
50
51 //! number of part to partition the graph
52 idx_t * nparts;
53
54 //! Desired weight for each partition (one for each constrain)
55 real_t * tpwgts;
56
57 //! For each partition load imbalance tollerated
58 real_t * ubvec;
59
60 //! Additional option for the graph partitioning
61 idx_t * options;
62
63 //! return the total comunication cost for each partition
64 idx_t * objval;
65
66 //! Is a output vector containing the partition for each vertex
67 idx_t * part;
68
69 //! Upon successful completion, the number of edges that are cut by the partitioning is written to this parameter.
70 idx_t * edgecut;
71
72 //! This parameter describes the ratio of inter-processor communication time compared to data redistri- bution time. It should be set between 0.000001 and 1000000.0. If ITR is set high, a repartitioning with a low edge-cut will be computed. If it is set low, a repartitioning that requires little data redistri- bution will be computed. Good values for this parameter can be obtained by dividing inter-processor communication time by data redistribution time. Otherwise, a value of 1000.0 is recommended.
73 real_t * itr;
74
75 //! This is used to indicate the numbering scheme that is used for the vtxdist, xadj, adjncy, and part arrays. (0 for C-style, start from 0 index)
76 idx_t * numflag;
77
78 //! This is used to indicate if the graph is weighted. wgtflag can take one of four values:
79 // 0 No weights (vwgt and adjwgt are both NULL).
80 // 1 Weights on the edges only (vwgt is NULL).
81 // 2 Weights on the vertices only (adjwgt is NULL).
82 // 3 Weights on both the vertices and edges.
83 idx_t * wgtflag;
84};
85
86//! Balance communication and computation
87#define BALANCE_CC 1
88//! Balance communication computation and memory
89#define BALANCE_CCM 2
90//! Balance computation and comunication and others
91#define BALANCE_CC_O(c) c+1
92
93/*! \brief Helper class to define Metis graph
94 *
95 * \tparam graph structure that store the graph
96 *
97 */
98template<typename Graph>
99class Parmetis
100{
101 //! Graph in metis reppresentation
102 Parmetis_graph Mg;
103
104 // Original graph
105 // Graph & g;
106
107 //! Communticator for OpenMPI
108 MPI_Comm comm = (MPI_Comm)NULL;
109
110 //! VCluster
111 Vcluster<> & v_cl;
112
113 //! Process rank information
114 int p_id = 0;
115
116 //! nc Number of partition
117 size_t nc = 0;
118
119 //! first re-mapped id
120 rid first;
121
122 //! last re-mapped id
123 rid last;
124
125 //! number of vertices that the processor has
126 size_t nvertex;
127
128 //! indicate how many time decompose/refine/re-decompose has been called
129 size_t n_dec;
130
131 //! Distribution tolerance
132 real_t dist_tol = 1.05;
133
134 /*! \brief Construct Adjacency list
135 *
136 * \param g Global graph
137 * \param m2g map from local index to global index
138 *
139 */
140 void constructAdjList(Graph &g, const std::unordered_map<rid, gid> & m2g)
141 {
142 // init basic graph informations and part vector
143 // Put the total communication size to NULL
144
145 Mg.nvtxs[0] = nvertex;
146 if (Mg.part != NULL)
147 {delete[] Mg.part;}
148 Mg.part = new idx_t[nvertex];
149
150 size_t nedge = 0;
151 size_t i = 0;
152 for (rid j = first; i < nvertex; i++, ++j)
153 {
154 Mg.part[i] = p_id;
155 nedge += g.getNChilds(m2g.find(j)->second.id);
156 }
157
158 // create xadj, adjlist, vwgt, adjwgt and vsize
159 if (Mg.xadj != NULL)
160 {delete[] Mg.xadj;}
161 if (Mg.adjncy != NULL)
162 {delete[] Mg.adjncy;}
163 if (Mg.vwgt != NULL)
164 {delete[] Mg.vwgt;}
165 if (Mg.adjwgt != NULL)
166 {delete[] Mg.adjwgt;}
167 if (Mg.vsize != NULL)
168 {delete[] Mg.vsize;}
169 Mg.xadj = new idx_t[nvertex + 1];
170 Mg.adjncy = new idx_t[nedge];
171 Mg.vwgt = new idx_t[nvertex];
172 Mg.adjwgt = new idx_t[nedge];
173 Mg.vsize = new idx_t[nvertex];
174
175 //! starting point in the adjacency list
176 size_t prev = 0;
177
178 // actual position
179 size_t id = 0;
180
181 size_t j = 0;
182
183 // for each vertex calculate the position of the starting point in the adjacency list
184 for (rid i = first; i <= last; ++i, j++)
185 {
186 gid idx = m2g.find(i)->second;
187
188 // Add weight to vertex and migration cost
189 Mg.vwgt[j] = g.vertex(idx.id).template get<nm_v_computation>();
190 Mg.vsize[j] = g.vertex(idx.id).template get<nm_v_migration>();
191
192 // Calculate the starting point in the adjacency list
193 Mg.xadj[id] = prev;
194
195 // Create the adjacency list and the weights for edges
196 for (size_t s = 0; s < g.getNChilds(idx.id); s++)
197 {
198
199 size_t child = g.getChild(idx.id, s);
200
201 Mg.adjncy[prev + s] = g.vertex(child).template get<nm_v_id>();
202 Mg.adjwgt[prev + s] = g.getChildEdge(idx.id, s).template get<nm_e::communication>();
203 }
204
205 // update the position for the next vertex
206 prev += g.getNChilds(idx.id);
207
208 id++;
209 }
210
211 // Fill the last point
212 Mg.xadj[id] = prev;
213 }
214
215public:
216
217 /*! \brief Constructor
218 *
219 * Construct a metis graph from Graph_CSR
220 *
221 * \param v_cl Vcluster object
222 * \param nc number of partitions
223 *
224 */
225 Parmetis(Vcluster<> & v_cl, size_t nc)
226 :v_cl(v_cl), nc(nc),n_dec(0)
227 {
228#ifdef SE_CLASS1
229
230 if (sizeof(idx_t) != 8)
231 {
232 std::cerr << __FILE__ << ":" << __LINE__ << " Error detected invalid installation of Parmetis. OpenFPM support Parmetis/Metis version with 64 bit idx_t" << std::endl;
233 ACTION_ON_ERROR(PARMETIS_ERROR_OBJECT);
234 }
235
236#endif
237
238 // TODO Move into VCluster
239 MPI_Comm_dup(MPI_COMM_WORLD, &comm);
240
241 // Nullify Mg
242 Mg.nvtxs = NULL;
243 Mg.ncon = NULL;
244 Mg.xadj = NULL;
245 Mg.adjncy = NULL;
246 Mg.vwgt = NULL;
247 Mg.vsize = NULL;
248 Mg.adjwgt = NULL;
249 Mg.nparts = NULL;
250 Mg.tpwgts = NULL;
251 Mg.ubvec = NULL;
252 Mg.options = NULL;
253 Mg.objval = NULL;
254 Mg.part = NULL;
255 Mg.edgecut = NULL;
256 Mg.itr = NULL;
257 Mg.numflag = NULL;
258 Mg.wgtflag = NULL;
259 first.id = 0;
260 last.id = 0;
261 nvertex = 0;
262 }
263
264 /*! \brief destructor
265 *
266 * Destructor, It destroy all the memory allocated
267 *
268 */
269 ~Parmetis()
270 {
271 // Deallocate the Mg structure
272 if (Mg.nvtxs != NULL)
273 {
274 delete[] Mg.nvtxs;
275 }
276
277 if (Mg.ncon != NULL)
278 {
279 delete[] Mg.ncon;
280 }
281
282 if (Mg.xadj != NULL)
283 {
284 delete[] Mg.xadj;
285 }
286
287 if (Mg.adjncy != NULL)
288 {
289 delete[] Mg.adjncy;
290 }
291
292 if (Mg.vwgt != NULL)
293 {
294 delete[] Mg.vwgt;
295 }
296
297 if (Mg.adjwgt != NULL)
298 {
299 delete[] Mg.adjwgt;
300 }
301
302 if (Mg.nparts != NULL)
303 {
304 delete[] Mg.nparts;
305 }
306
307 if (Mg.tpwgts != NULL)
308 {
309 delete[] Mg.tpwgts;
310 }
311
312 if (Mg.ubvec != NULL)
313 {
314 delete[] Mg.ubvec;
315 }
316
317 if (Mg.options != NULL)
318 {
319 delete[] Mg.options;
320 }
321
322 if (Mg.part != NULL)
323 {
324 delete[] Mg.part;
325 }
326
327 if (Mg.edgecut != NULL)
328 {
329 delete[] Mg.edgecut;
330 }
331
332 if (Mg.numflag != NULL)
333 {
334 delete[] Mg.numflag;
335 }
336
337 if (Mg.wgtflag != NULL)
338 {
339 delete[] Mg.wgtflag;
340 }
341
342 if (Mg.itr != NULL)
343 {
344 delete[] Mg.itr;
345 }
346
347 if (Mg.vsize != NULL)
348 {
349 delete[] Mg.vsize;
350 }
351
352 if (Mg.objval != NULL)
353 {
354 delete[] Mg.objval;
355 }
356
357 if (is_openfpm_init() == true)
358 {MPI_Comm_free(&comm);}
359 }
360
361 /*! \brief Set the Sub-graph
362 *
363 * \param g Global graph to set
364 * \param vtxdist indicate how the vertex of the graph are distrubuted across
365 * processors.
366 * \param m2g map the local ids of the vertex into global-ids
367 * \param w true if vertices have weights
368 */
369 void initSubGraph(Graph & g,
370 const openfpm::vector<rid> & vtxdist,
371 const std::unordered_map<rid, gid> & m2g,
372 bool w)
373 {
374 p_id = v_cl.getProcessUnitID();
375
376 first = vtxdist.get(p_id);
377 last = vtxdist.get(p_id + 1) - 1;
378 nvertex = last.id - first.id + 1;
379
380 setDefaultParameters(w);
381
382 // construct the adjacency list
383 constructAdjList(g, m2g);
384 }
385
386 /*! \brief Decompose the graph
387 *
388 * \tparam i which property store the decomposition
389 *
390 */
391 void decompose(const openfpm::vector<rid> & vtxdist)
392 {
393 // Decompose
394
395 ParMETIS_V3_PartKway((idx_t *) vtxdist.getPointer(), Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.adjwgt, Mg.wgtflag, Mg.numflag, Mg.ncon, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.options, Mg.edgecut, Mg.part, &comm);
396
397 n_dec++;
398 }
399
400 /*! \brief Refine the graph
401 *
402 * \tparam i which property store the refined decomposition
403 *
404 */
405 void refine(openfpm::vector<rid> & vtxdist)
406 {
407 // Refine
408
409 ParMETIS_V3_RefineKway((idx_t *) vtxdist.getPointer(), Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.adjwgt, Mg.wgtflag, Mg.numflag, Mg.ncon, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.options, Mg.edgecut, Mg.part, &comm);
410
411 n_dec++;
412 }
413
414 /*! \brief Redecompose the graph
415 *
416 * \tparam i which property
417 *
418 */
419 void redecompose(openfpm::vector<rid> & vtxdist)
420 {
421 ParMETIS_V3_AdaptiveRepart((idx_t *)vtxdist.getPointer(), Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.vsize, Mg.adjwgt, Mg.wgtflag, Mg.numflag, Mg.ncon, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.itr, Mg.options, Mg.edgecut, Mg.part, &comm);
422
423 n_dec++;
424 }
425
426 /*! \brief Get graph partition vector
427 *
428 */
429 idx_t* getPartition()
430 {
431 return Mg.part;
432 }
433
434 /*! \brief Reset graph and reconstruct it
435 *
436 * \param g Global graph
437 * \param vtxdist Distribution vector
438 * \param m2g Mapped id to global id map
439 * \param vgw Using weights on vertices
440 */
441 void reset(Graph & g, const openfpm::vector<rid> & vtxdist, const std::unordered_map<rid, gid> & m2g, bool vgw)
442 {
443 first = vtxdist.get(p_id);
444 last = vtxdist.get(p_id + 1) - 1;
445 nvertex = last.id - first.id + 1;
446
447 // Deallocate the graph structures
448
449 setDefaultParameters(vgw);
450
451 // construct the adjacency list
452 constructAdjList(g, m2g);
453 }
454
455 /*! \brief Seth the default parameters for parmetis
456 *
457 *
458 */
459 void setDefaultParameters(bool w)
460 {
461 if (Mg.nvtxs != NULL)
462 {delete[] Mg.nvtxs;}
463 Mg.nvtxs = new idx_t[1];
464
465 // Set the number of constrains
466 if (Mg.ncon != NULL)
467 {delete[] Mg.ncon;}
468 Mg.ncon = new idx_t[1];
469 Mg.ncon[0] = 1;
470
471 // Set to null the weight of the vertex (init after in constructAdjList) (can be removed)
472 if (Mg.vwgt != NULL)
473 {delete[] Mg.vwgt;}
474 Mg.vwgt = NULL;
475
476 // Set to null the weight of the edge (init after in constructAdjList) (can be removed)
477 if (Mg.adjwgt != NULL)
478 {delete[] Mg.adjwgt;}
479 Mg.adjwgt = NULL;
480
481 // Set the total number of partitions
482 if (Mg.nparts != NULL)
483 {delete[] Mg.nparts;}
484 Mg.nparts = new idx_t[1];
485 Mg.nparts[0] = nc;
486
487 //! Set option for the graph partitioning (set as default)
488
489 if (Mg.options != NULL)
490 {delete[] Mg.options;}
491 Mg.options = new idx_t[4];
492 Mg.options[0] = 0;
493 Mg.options[1] = 0;
494 Mg.options[2] = 0;
495 Mg.options[3] = 0;
496
497 //! is an output vector containing the partition for each vertex
498
499 //! adaptiveRepart itr value
500 if(Mg.itr != NULL)
501 {delete[] Mg.itr;}
502 Mg.itr = new real_t[1];
503 Mg.itr[0] = 1000.0;
504
505 if (Mg.objval != NULL)
506 {delete[] Mg.objval;}
507 Mg.objval = new idx_t[1];
508
509 //! init tpwgts to have balanced vertices and ubvec
510
511 if (Mg.tpwgts != NULL)
512 {delete[] Mg.tpwgts;}
513 if (Mg.ubvec != NULL)
514 {delete[] Mg.ubvec;}
515 Mg.tpwgts = new real_t[Mg.nparts[0]];
516 Mg.ubvec = new real_t[Mg.nparts[0]];
517
518 for (size_t s = 0; s < (size_t) Mg.nparts[0]; s++)
519 {
520 Mg.tpwgts[s] = 1.0 / Mg.nparts[0];
521 Mg.ubvec[s] = dist_tol;
522 }
523
524 if (Mg.edgecut != NULL)
525 {delete[] Mg.edgecut;}
526 Mg.edgecut = new idx_t[1];
527 Mg.edgecut[0] = 0;
528
529 //! This is used to indicate the numbering scheme that is used for the vtxdist, xadj, adjncy, and part arrays. (0 for C-style, start from 0 index)
530 if (Mg.numflag != NULL)
531 {delete[] Mg.numflag;}
532 Mg.numflag = new idx_t[1];
533 Mg.numflag[0] = 0;
534
535 //! This is used to indicate if the graph is weighted. wgtflag can take one of four values:
536 if (Mg.wgtflag != NULL)
537 {delete[] Mg.wgtflag;}
538 Mg.wgtflag = new idx_t[1];
539
540 if (w)
541 Mg.wgtflag[0] = 3;
542 else
543 Mg.wgtflag[0] = 0;
544 }
545
546 /*! \brief Copy the object
547 *
548 * \param pm object to copy
549 *
550 * \return itself
551 *
552 */
553 const Parmetis<Graph> & operator=(const Parmetis<Graph> & pm)
554 {
555 MPI_Comm_dup(pm.comm, &comm);
556 p_id = pm.p_id;
557 nc = pm.nc;
558 n_dec = pm.n_dec;
559 dist_tol = pm.dist_tol;
560
561 setDefaultParameters(pm.Mg.wgtflag[0] == 3);
562
563 return *this;
564 }
565
566 /*! \brief Copy the object
567 *
568 * \param pm object to copy
569 *
570 * \return itself
571 *
572 */
573 const Parmetis<Graph> & operator=(Parmetis<Graph> && pm)
574 {
575 // TODO Move into VCluster
576 MPI_Comm_dup(pm.comm, &comm);
577 p_id = pm.p_id;
578 nc = pm.nc;
579 n_dec = pm.n_dec;
580 dist_tol = pm.dist_tol;
581
582 setDefaultParameters(pm.Mg.wgtflag[0] == 3);
583
584 return *this;
585 }
586
587 /*! \brief Get the decomposition counter
588 *
589 * \return the decomposition counter
590 *
591 */
592 size_t get_ndec()
593 {
594 return n_dec;
595 }
596
597 /*! \brief Distribution tolerance
598 *
599 * \param tol tolerance
600 *
601 */
602 void setDistTol(real_t tol)
603 {
604 dist_tol = tol;
605 }
606};
607
608#endif
609