1/*
2 * parmetis_util.hpp
3 *
4 * Created on: Oct 07, 2015
5 * Author: Antonio Leo
6 */
7
8#ifndef DISTPARMETIS_UTIL_HPP
9#define DISTPARMETIS_UTIL_HPP
10
11#include <iostream>
12#include "parmetis.h"
13#include "VTKWriter/VTKWriter.hpp"
14#include "VCluster/VCluster.hpp"
15
16/*! \brief Metis graph structure
17 *
18 * Metis graph structure
19 *
20 */
21struct Parmetis_dist_graph
22{
23 //! The number of vertices in the graph
24 idx_t * nvtxs;
25
26 //! number of balancing constrains
27 //! more practical, are the number of weights for each vertex
28 //! PS even we you specify vwgt == NULL ncon must be set at leat to
29 //! one
30 idx_t * ncon;
31
32 //! For each vertex it store the adjacency lost start for the vertex i
33 idx_t * xadj;
34
35 //! For each vertex it store a list of all neighborhood vertex
36 idx_t * adjncy;
37
38 //! Array that store the weight for each vertex
39 idx_t * vwgt;
40
41 //! Array of the vertex size, basically is the total communication amount
42 idx_t * vsize;
43
44 //! The weight of the edge
45 idx_t * adjwgt;
46
47 //! number of part to partition the graph
48 idx_t * nparts;
49
50 //! Desired weight for each partition (one for each constrain)
51 real_t * tpwgts;
52
53 //! For each partition load imbalance tollerated
54 real_t * ubvec;
55
56 //! Additional option for the graph partitioning
57 idx_t * options;
58
59 //! return the total comunication cost for each partition
60 idx_t * objval;
61
62 //! Is a output vector containing the partition for each vertex
63 idx_t * part;
64
65 //! Upon successful completion, the number of edges that are cut by the partitioning is written to this parameter.
66 idx_t * edgecut;
67
68 //! 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.
69 real_t * itr;
70
71 //! 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)
72 idx_t * numflag;
73
74 //! This is used to indicate if the graph is weighted. wgtflag can take one of four values:
75 //! 0 No weights (vwgt and adjwgt are both NULL).
76 //! 1 Weights on the edges only (vwgt is NULL).
77 //! 2 Weights on the vertices only (adjwgt is NULL).
78 //! 3 Weights on both the vertices and edges.
79 idx_t * wgtflag;
80};
81
82//! Balance communication and computation
83#define BALANCE_CC 1
84//! Balance communication computation and memory
85#define BALANCE_CCM 2
86//! Balance computation and comunication and others
87#define BALANCE_CC_O(c) c+1
88
89/*! \brief Helper class to define Metis graph
90 *
91 * TODO Transform pointer to openfpm vector
92 *
93 * \tparam graph structure that store the graph
94 *
95 */
96template<typename Graph>
97class DistParmetis
98{
99 //! Graph in parmetis reppresentation
100 Parmetis_dist_graph Mg;
101
102 //! Communticator for OpenMPI
103 MPI_Comm comm = (MPI_Comm)NULL;
104
105 //! VCluster
106 Vcluster<> & v_cl;
107
108 //! Process rank information
109 int p_id = 0;
110
111 //! nc Number of partition
112 size_t nc = 0;
113
114 /*! \brief Construct Adjacency list
115 *
116 * \param sub_g graph in which we construct the adjacency list
117 *
118 */
119 void constructAdjList(Graph & sub_g)
120 {
121
122 // init basic graph informations and part vector
123 Mg.nvtxs[0] = sub_g.getNVertex();
124 Mg.part = new idx_t[sub_g.getNVertex()];
125 for (size_t i = 0; i < sub_g.getNVertex(); i++)
126 Mg.part[i] = p_id;
127
128 // create xadj, adjlist, vwgt, adjwgt and vsize
129 Mg.xadj = new idx_t[sub_g.getNVertex() + 1];
130 Mg.adjncy = new idx_t[sub_g.getNEdge()];
131 Mg.vwgt = new idx_t[sub_g.getNVertex()];
132 Mg.adjwgt = new idx_t[sub_g.getNEdge()];
133 Mg.vsize = new idx_t[sub_g.getNVertex()];
134
135 //! starting point in the adjacency list
136 size_t prev = 0;
137
138 // actual position
139 size_t id = 0;
140
141 for (size_t i = 0, j = sub_g.firstId(); i < sub_g.getNVertex() && j <= sub_g.lastId(); i++, j++)
142 {
143 size_t idx = sub_g.nodeById(j);
144
145 // Add weight to vertex and migration cost
146 Mg.vwgt[i] = sub_g.vertex(idx).template get<nm_v_computation>();
147 Mg.vsize[i] = sub_g.vertex(idx).template get<nm_v_migration>();
148
149 // Calculate the starting point in the adjacency list
150 Mg.xadj[id] = prev;
151
152 // Create the adjacency list and the weights for edges
153 for (size_t s = 0; s < sub_g.getNChilds(idx); s++)
154 {
155 Mg.adjncy[prev + s] = sub_g.getChild(idx, s);
156
157 Mg.adjwgt[prev + s] = sub_g.getChildEdge(idx, s).template get<nm_e::communication>();
158 }
159
160 // update the position for the next vertex
161 prev += sub_g.getNChilds(idx);
162
163 id++;
164 }
165
166 // Fill the last point
167 Mg.xadj[id] = prev;
168
169 }
170
171public:
172
173 /*! \brief Constructor
174 *
175 * Construct a metis graph from a Graph_CSR
176 *
177 * \param v_cl Vcluster
178 * \param nc number of partitions
179 *
180 */
181 DistParmetis(Vcluster<> & v_cl, size_t nc) :
182 v_cl(v_cl), nc(nc)
183 {
184 // TODO Move into VCluster
185 MPI_Comm_dup(MPI_COMM_WORLD, &comm);
186
187 // Nullify Mg
188 Mg.nvtxs = NULL;
189 Mg.ncon = NULL;
190 Mg.xadj = NULL;
191 Mg.adjncy = NULL;
192 Mg.vwgt = NULL;
193 Mg.vsize = NULL;
194 Mg.adjwgt = NULL;
195 Mg.nparts = NULL;
196 Mg.tpwgts = NULL;
197 Mg.ubvec = NULL;
198 Mg.options = NULL;
199 Mg.objval = NULL;
200 Mg.part = NULL;
201 Mg.edgecut = NULL;
202 Mg.itr = NULL;
203 Mg.numflag = NULL;
204 Mg.wgtflag = NULL;
205 }
206
207 //TODO deconstruct new variables
208 /*! \brief destructor
209 *
210 * Destructor, It destroy all the memory allocated
211 *
212 */
213 ~DistParmetis()
214 {
215 // Deallocate the Mg structure
216 if (Mg.nvtxs != NULL)
217 {
218 delete[] Mg.nvtxs;
219 }
220
221 if (Mg.ncon != NULL)
222 {
223 delete[] Mg.ncon;
224 }
225
226 if (Mg.xadj != NULL)
227 {
228 delete[] Mg.xadj;
229 }
230
231 if (Mg.adjncy != NULL)
232 {
233 delete[] Mg.adjncy;
234 }
235
236 if (Mg.vwgt != NULL)
237 {
238 delete[] Mg.vwgt;
239 }
240
241 if (Mg.adjwgt != NULL)
242 {
243 delete[] Mg.adjwgt;
244 }
245
246 if (Mg.nparts != NULL)
247 {
248 delete[] Mg.nparts;
249 }
250
251 if (Mg.tpwgts != NULL)
252 {
253 delete[] Mg.tpwgts;
254 }
255
256 if (Mg.ubvec != NULL)
257 {
258 delete[] Mg.ubvec;
259 }
260
261 if (Mg.options != NULL)
262 {
263 delete[] Mg.options;
264 }
265
266 if (Mg.part != NULL)
267 {
268 delete[] Mg.part;
269 }
270
271 if (Mg.edgecut != NULL)
272 {
273 delete[] Mg.edgecut;
274 }
275
276 if (Mg.numflag != NULL)
277 {
278 delete[] Mg.numflag;
279 }
280
281 if (Mg.wgtflag != NULL)
282 {
283 delete[] Mg.wgtflag;
284 }
285
286 if (Mg.itr != NULL)
287 {
288 delete[] Mg.itr;
289 }
290
291 if (Mg.vsize != NULL)
292 {
293 delete[] Mg.vsize;
294 }
295 }
296
297 /*! \brief Set the Sub-graph
298 *
299 * \param sub_g Sub-graph to set
300 *
301 */
302 void initSubGraph(Graph & sub_g)
303 {
304 p_id = v_cl.getProcessUnitID();
305
306 // Get the number of vertex
307
308 if (Mg.nvtxs != NULL)
309 {delete[] Mg.nvtxs;}
310 Mg.nvtxs = new idx_t[1];
311 Mg.nvtxs[0] = sub_g.getNVertex();
312
313 // Set the number of constrains
314
315 if (Mg.ncon != NULL)
316 {delete[] Mg.ncon;}
317 Mg.ncon = new idx_t[1];
318 Mg.ncon[0] = 1;
319
320 // Set to null the weight of the vertex (init after in constructAdjList) (can be removed)
321
322 if (Mg.vwgt != NULL)
323 {delete[] Mg.vwgt;}
324 Mg.vwgt = NULL;
325
326 // Set to null the weight of the edge (init after in constructAdjList) (can be removed)
327
328 if (Mg.adjwgt != NULL)
329 {delete[] Mg.adjwgt;}
330 Mg.adjwgt = NULL;
331
332 // construct the adjacency list
333
334 constructAdjList(sub_g);
335
336 // Set the total number of partitions
337
338 if (Mg.nparts != NULL)
339 {delete[] Mg.nparts;}
340 Mg.nparts = new idx_t[1];
341 Mg.nparts[0] = nc;
342
343 //! Set option for the graph partitioning (set as default)
344
345 if (Mg.options != NULL)
346 {delete[] Mg.options;}
347 Mg.options = new idx_t[4];
348 Mg.options[0] = 0;
349 Mg.options[1] = 0;
350 Mg.options[2] = 0;
351 Mg.options[3] = 0;
352
353 //! adaptiveRepart itr value
354 if (Mg.itr != NULL)
355 {delete[] Mg.itr;}
356 Mg.itr = new real_t[1];
357 Mg.itr[0] = 1000.0;
358
359 //! init tpwgts to have balanced vertices and ubvec
360
361 if (Mg.tpwgts != NULL)
362 {delete[] Mg.tpwgts;}
363 if (Mg.ubvec != NULL)
364 {delete[] Mg.ubvec;}
365 Mg.tpwgts = new real_t[Mg.nparts[0]];
366 Mg.ubvec = new real_t[Mg.nparts[0]];
367
368 for (int s = 0; s < Mg.nparts[0]; s++)
369 {
370 Mg.tpwgts[s] = 1.0 / Mg.nparts[0];
371 Mg.ubvec[s] = 1.05;
372 }
373
374 if (Mg.edgecut != NULL)
375 {delete[] Mg.edgecut;}
376 Mg.edgecut = new idx_t[1];
377 Mg.edgecut[0] = 0;
378
379 //! 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)
380 if (Mg.numflag != NULL)
381 {delete[] Mg.numflag;}
382 Mg.numflag = new idx_t[1];
383 Mg.numflag[0] = 0;
384
385 //! This is used to indicate if the graph is weighted.
386 if (Mg.wgtflag != NULL)
387 {delete[] Mg.wgtflag;}
388 Mg.wgtflag = new idx_t[1];
389 Mg.wgtflag[0] = 3;
390 }
391
392 /*! \brief Decompose the graph
393 *
394 * \tparam i which property store the decomposition
395 *
396 * \param sub_g graph to decompose
397 *
398 */
399 template<unsigned int i>
400 void decompose(Graph & sub_g)
401 {
402
403 // Decompose
404
405 ParMETIS_V3_PartKway((idx_t *) sub_g.getVtxdist()->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);
406 /*
407 ParMETIS_V3_AdaptiveRepart( (idx_t *) vtxdist.getPointer(), Mg.xadj,Mg.adjncy,Mg.vwgt,Mg.vsize,Mg.adjwgt, Mg.wgtflag, Mg.numflag,
408 Mg.ncon, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.itr, Mg.options, Mg.edgecut,
409 Mg.part, &comm );
410 */
411
412 // For each vertex store the processor that contain the data
413 for (size_t id = 0, j = sub_g.firstId(); id < sub_g.getNVertex() && j <= sub_g.lastId(); id++, j++)
414 {
415 sub_g.vertex(sub_g.nodeById(j)).template get<i>() = Mg.part[id];
416 }
417 }
418
419 /*! \brief Refine the graph
420 *
421 * \tparam i which property store the refined decomposition
422 *
423 * \param sub_g graph to decompose
424 *
425 */
426
427 template<unsigned int i>
428 void refine(Graph & sub_g)
429 {
430 // Refine
431 //ParMETIS_V3_PartKway((idx_t *) sub_g.getVtxdist()->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);
432 ParMETIS_V3_AdaptiveRepart((idx_t *) sub_g.getVtxdist()->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);
433
434 // For each vertex store the processor that contain the data
435 for (size_t id = 0, j = sub_g.firstId(); id < sub_g.getNVertex() && j <= sub_g.lastId(); id++, j++)
436 {
437 sub_g.vertex(sub_g.nodeById(j)).template get<i>() = Mg.part[id];
438 }
439 }
440
441 /*! \brief Get graph partition vector
442 *
443 * \return the partition or the assignment of each sub-sub-domain
444 *
445 */
446 idx_t * getPartition()
447 {
448 return Mg.part;
449 }
450
451 /*! \brief Reset graph and reconstruct it
452 *
453 * \param sub_g graph to decompose
454 *
455 */
456 void reset(Graph & sub_g)
457 {
458 // Deallocate the graph structures
459
460 if (Mg.xadj != NULL)
461 {
462 delete[] Mg.xadj;
463 }
464
465 if (Mg.adjncy != NULL)
466 {
467 delete[] Mg.adjncy;
468 }
469
470 if (Mg.vwgt != NULL)
471 {
472 delete[] Mg.vwgt;
473 }
474
475 if (Mg.adjwgt != NULL)
476 {
477 delete[] Mg.adjwgt;
478 }
479
480 if (Mg.part != NULL)
481 {
482 delete[] Mg.part;
483 }
484
485 if (Mg.vsize != NULL)
486 {
487 delete[] Mg.vsize;
488 }
489
490 sub_g.deleteGhosts();
491
492 constructAdjList(sub_g);
493 }
494}
495;
496
497#endif
498