ParMetisDistribution.hpp 15.5 KB
Newer Older
Pietro Incardona's avatar
Pietro Incardona committed
1 2 3 4
/*
 * ParMetisDistribution.hpp
 *
 *  Created on: Nov 19, 2015
5
 *      Author: Antonio Leo
Pietro Incardona's avatar
Pietro Incardona committed
6 7
 */

8 9 10 11 12

#ifndef SRC_DECOMPOSITION_PARMETISDISTRIBUTION_HPP_
#define SRC_DECOMPOSITION_PARMETISDISTRIBUTION_HPP_


13 14
#include "SubdomainGraphNodes.hpp"
#include "parmetis_util.hpp"
15
#include "Graph/ids.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
16

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
#define PARMETIS_DISTRIBUTION_ERROR 100002

/*! \brief Class that distribute sub-sub-domains across processors using ParMetis Library
 *
 * Given a graph and setting Computational cost, Communication cost (on the edge) and
 * Migration cost or total Communication costs, it produce the optimal balanced distribution
 *
 * In addition to Metis it provide the functionality to refine the previously computed
 * decomposition
 *
 * ### Initialize a Cartesian graph and decompose
 * \snippet Distribution_unit_tests.hpp Initialize a ParMetis Cartesian graph and decompose
 *
 * ### Refine the decomposition
 * \snippet Distribution_unit_tests.hpp refine with parmetis the decomposition
 *
 */
34
template<unsigned int dim, typename T>
Pietro Incardona's avatar
Pietro Incardona committed
35 36 37 38 39
class ParMetisDistribution
{
	//! Vcluster
	Vcluster & v_cl;

40 41
	//! Structure that store the cartesian grid information
	grid_sm<dim, void> gr;
Pietro Incardona's avatar
Pietro Incardona committed
42

43
	//! rectangular domain to decompose
44
	Box<dim, T> domain;
Pietro Incardona's avatar
Pietro Incardona committed
45 46 47 48

	//! Global sub-sub-domain graph
	Graph_CSR<nm_v, nm_e> gp;

49 50 51
	//! Convert the graph to parmetis format
	Parmetis<Graph_CSR<nm_v, nm_e>> parmetis_graph;

Pietro Incardona's avatar
Pietro Incardona committed
52
	//! Init vtxdist needed for Parmetis
53 54 55 56 57 58 59 60 61 62 63 64 65 66
	//
	// vtxdist is a common array across processor, it indicate how
	// vertex are distributed across processors
	//
	// Example we have 3 processors
	//
	// processor 0 has 3 vertices
	// processor 1 has 5 vertices
	// processor 2 has 4 vertices
	//
	// vtxdist contain, 0,3,8,12
	//
	// vtx dist is the unique global-id of the vertices
	//
67
	openfpm::vector<rid> vtxdist;
Pietro Incardona's avatar
Pietro Incardona committed
68 69 70 71 72

	//! partitions
	openfpm::vector<openfpm::vector<idx_t>> partitions;

	//! Init data structure to keep trace of new vertices distribution in processors (needed to update main graph)
73
	openfpm::vector<openfpm::vector<gid>> v_per_proc;
Pietro Incardona's avatar
Pietro Incardona committed
74

75
	//! Hashmap to access to the global position given the re-mapped one (needed for access the map)
76
	std::unordered_map<rid, gid> m2g;
Pietro Incardona's avatar
Pietro Incardona committed
77

78 79
	//! Flag to check if weights are used on vertices
	bool verticesGotWeights = false;
Pietro Incardona's avatar
Pietro Incardona committed
80

81
	/*! \brief Update main graph ad subgraph with the received data of the partitions from the other processors
Pietro Incardona's avatar
Pietro Incardona committed
82 83 84 85
	 *
	 */
	void updateGraphs()
	{
86
		size_t Np = v_cl.getProcessingUnits();
Pietro Incardona's avatar
Pietro Incardona committed
87 88

		// Init n_vtxdist to gather informations about the new decomposition
89
		openfpm::vector<rid> n_vtxdist(Np + 1);
90
		for (size_t i = 0; i <= Np; i++)
91
			n_vtxdist.get(i).id = 0;
Pietro Incardona's avatar
Pietro Incardona committed
92

93 94
		// Update the main graph with received data from processor i
		for (size_t i = 0; i < Np; i++)
Pietro Incardona's avatar
Pietro Incardona committed
95
		{
96
			size_t ndata = partitions.get(i).size();
97
			size_t k = 0;
Pietro Incardona's avatar
Pietro Incardona committed
98

99
			// Update the main graph with the received informations
100
			for (rid l = vtxdist.get(i); k < ndata && l < vtxdist.get(i + 1); k++, ++l)
Pietro Incardona's avatar
Pietro Incardona committed
101
			{
102
				// Create new n_vtxdist (just count processors vertices)
103
				++n_vtxdist.get(partitions.get(i).get(k) + 1);
Pietro Incardona's avatar
Pietro Incardona committed
104

105
				// Update proc id in the vertex (using the old map)
106
				vertexByMapId(l).template get<nm_v::proc_id>() = partitions.get(i).get(k);
Pietro Incardona's avatar
Pietro Incardona committed
107 108

				// Add vertex to temporary structure of distribution (needed to update main graph)
109
				v_per_proc.get(partitions.get(i).get(k)).add(getVertexGlobalId(l));
Pietro Incardona's avatar
Pietro Incardona committed
110 111 112
			}
		}

113 114
		// Create new n_vtxdist (accumulate the counters)
		for (size_t i = 2; i <= Np; i++)
Pietro Incardona's avatar
Pietro Incardona committed
115 116 117
			n_vtxdist.get(i) += n_vtxdist.get(i - 1);

		// Copy the new decomposition in the main vtxdist
118
		for (size_t i = 0; i <= Np; i++)
Pietro Incardona's avatar
Pietro Incardona committed
119 120
			vtxdist.get(i) = n_vtxdist.get(i);

121
		// Renumber the main graph and re-create the map
122
		for (size_t p = 0; p < (size_t)Np; p++)
Pietro Incardona's avatar
Pietro Incardona committed
123
		{
124 125
			size_t i = 0;
			for (rid j = vtxdist.get(p); j < vtxdist.get(p + 1); ++j, i++)
Pietro Incardona's avatar
Pietro Incardona committed
126
			{
127
				setMapId(j, v_per_proc.get(p).get(i));
128
				gp.vertex(v_per_proc.get(p).get(i).id).template get<nm_v::id>() = j.id;
Pietro Incardona's avatar
Pietro Incardona committed
129 130 131 132
			}
		}
	}

133 134 135 136 137 138 139 140 141 142 143
	void createMapsFromGlobalGraph(openfpm::vector<size_t> & vtxdist)
	{
/*		openfpm::vector<size_t> cnt_np;

		for (size_t i = 0 ; i < gp.getNVertex() ; i++)
		{
			cnt_np(gp.template vertex<nm_v::proc_id>)++;

			gp.setMapId()
		}*/
	}
Pietro Incardona's avatar
Pietro Incardona committed
144

145 146 147 148 149 150 151
	/*! \brief operator to access the vertex by mapped position
	 *
	 * operator to access the vertex
	 *
	 * \param id re-mapped id of the vertex to access
	 *
	 */
152
	inline auto vertexByMapId(rid id) -> decltype( gp.vertex(m2g.find(id)->second.id) )
153
	{
154
		return gp.vertex(m2g.find(id)->second.id);
155
	}
Pietro Incardona's avatar
Pietro Incardona committed
156

157 158 159 160 161 162
	/*! \brief operator to remap vertex to a new position
	 *
	 * \param n re-mapped position
	 * \param g global position
	 *
	 */
163
	inline void setMapId(rid n, gid g)
164 165 166 167 168 169 170 171 172 173
	{
		m2g[n] = g;
	}

	/*! \brief Get the global id of the vertex given the re-mapped one
	 *
	 * \param remapped id
	 * \return global id
	 *
	 */
174
	gid getVertexGlobalId(rid n)
175 176 177 178 179 180 181 182 183 184 185
	{
		return m2g.find(n)->second;
	}

	/*! \brief operator to init ids vector
	 *
	 * operator to init ids vector
	 *
	 */
	void initLocalToGlobalMap()
	{
186 187 188 189
		gid g;
		rid i;
		i.id = 0;

190
		m2g.clear();
191
		for ( ; (size_t)i.id < gp.getNVertex(); ++i)
192
		{
193 194 195
			g.id = i.id;

			m2g.insert( { i, g });
196
		}
Pietro Incardona's avatar
Pietro Incardona committed
197 198
	}

199 200 201 202 203 204 205 206 207
	/*! \brief Callback of the sendrecv to set the size of the array received
	 *
	 * \param msg_i Index of the message
	 * \param total_msg Total numeber of messages
	 * \param total_p Total number of processors to comunicate with
	 * \param i Processor id
	 * \param ri Request id
	 * \param ptr Void pointer parameter for additional data to pass to the call-back
	 */
Pietro Incardona's avatar
Pietro Incardona committed
208 209 210 211 212 213 214 215 216 217 218
	static void * message_receive(size_t msg_i, size_t total_msg, size_t total_p, size_t i, size_t ri, void * ptr)
	{
		openfpm::vector < openfpm::vector < idx_t >> *v = static_cast<openfpm::vector<openfpm::vector<idx_t>> *>(ptr);

		v->get(i).resize(msg_i / sizeof(idx_t));

		return &(v->get(i).get(0));
	}

public:

219 220
	/*! Constructor for the ParMetis class
	 *
221
	 * \param v_cl Vcluster to use as communication object in this class
222
	 */
223 224 225 226
	ParMetisDistribution(Vcluster & v_cl)
	:v_cl(v_cl), parmetis_graph(v_cl, v_cl.getProcessingUnits()), vtxdist(v_cl.getProcessingUnits() + 1), partitions(v_cl.getProcessingUnits()), v_per_proc(v_cl.getProcessingUnits())
	{
	}
227

228 229 230 231
	/*! Copy constructor
	 *
	 * \param pm Distribution to copy
	 *
232
	 */
233 234
	ParMetisDistribution(const ParMetisDistribution<dim,T> & pm)
	:v_cl(pm.v_cl),parmetis_graph(v_cl, v_cl.getProcessingUnits())
235
	{
236
		this->operator=(pm);
237 238
	}

239 240 241
	/*! Copy constructor
	 *
	 * \param pm Distribution to copy
242 243
	 *
	 */
244
	ParMetisDistribution(ParMetisDistribution<dim,T> && pm)
245
	{
246
		this->operator=(pm);
247 248
	}

249
	/*! \brief Create the Cartesian graph
250
	 *
251 252
	 * \param grid info
	 * \param dom domain
253
	 */
254
	void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
255
	{
256 257 258 259 260
		size_t bc[dim];

		for (size_t i = 0 ; i < dim ; i++)
			bc[i] = NON_PERIODIC;

261 262 263 264 265 266
		// Set grid and domain
		gr = grid;
		domain = dom;

		// Create a cartesian grid graph
		CartesianGraphFactory<dim, Graph_CSR<nm_v, nm_e>> g_factory_part;
267
		gp = g_factory_part.template construct<NO_EDGE, nm_v::id, T, dim - 1, 0>(gr.getSize(), domain, bc);
268
		initLocalToGlobalMap();
269

270 271 272 273 274 275 276 277
		//! Get the number of processing units
		size_t Np = v_cl.getProcessingUnits();

		//! Division of vertices in Np graphs
		//! Put (div+1) vertices in mod graphs
		//! Put div vertices in the rest of the graphs
		size_t mod_v = gr.size() % Np;
		size_t div_v = gr.size() / Np;
278

279 280 281
		for (size_t i = 0; i <= Np; i++)
		{
			if (i < mod_v)
282
				vtxdist.get(i).id = (div_v + 1) * i;
283
			else
284
				vtxdist.get(i).id = (div_v) * i + mod_v;
285
		}
286 287 288 289 290 291

		// Init to 0.0 axis z (to fix in graphFactory)
		if (dim < 3)
		{
			for (size_t i = 0; i < gp.getNVertex(); i++)
			{
292
				gp.vertex(i).template get<nm_v::x>()[2] = 0.0;
293 294
			}
		}
295 296 297 298 299
		for (size_t i = 0; i < gp.getNVertex(); i++)
		{
			gp.vertex(i).template get<nm_v::global_id>() = i;
		}

300 301 302 303 304 305 306 307 308 309
	}

	/*! \brief Get the current graph (main)
	 *
	 */
	Graph_CSR<nm_v, nm_e> & getGraph()
	{
		return gp;
	}

310
	/*! \brief Create the decomposition
311 312 313 314
	 *
	 */
	void decompose()
	{
315

316 317 318 319 320 321
		//! Get the processor id
		size_t p_id = v_cl.getProcessUnitID();

		//! Get the number of processing units
		size_t Np = v_cl.getProcessingUnits();

322
		// Number of local vertex
323
		size_t nl_vertex = vtxdist.get(p_id+1).id - vtxdist.get(p_id).id;
324

325
		parmetis_graph.initSubGraph(gp, vtxdist, m2g, verticesGotWeights);
326 327

		//! Decompose
328
		parmetis_graph.decompose<nm_v::proc_id>(vtxdist);
329 330 331 332 333

		//! Get result partition for this processors
		idx_t *partition = parmetis_graph.getPartition();

		//! Prepare vector of arrays to contain all partitions
334 335
		partitions.get(p_id).resize(nl_vertex);
		std::copy(partition, partition + nl_vertex, &partitions.get(p_id).get(0));
336

337 338
		// Communicate the local distribution to the other processors
		// to reconstruct individually the global graph
339 340 341 342 343 344 345 346 347
		openfpm::vector<size_t> prc;
		openfpm::vector<size_t> sz;
		openfpm::vector<void *> ptr;

		for (size_t i = 0; i < Np; i++)
		{
			if (i != v_cl.getProcessUnitID())
			{
				prc.add(i);
348
				sz.add(nl_vertex * sizeof(idx_t));
349 350 351
				ptr.add(partitions.get(p_id).getPointer());
			}
		}
Pietro Incardona's avatar
Pietro Incardona committed
352

353 354 355
		v_cl.sendrecvMultipleMessagesNBX(prc.size(), &sz.get(0), &prc.get(0), &ptr.get(0), message_receive, &partitions,
		NONE);

356
		// Update graphs with the received data
357 358 359 360
		updateGraphs();
	}

	/*! \brief Refine current decomposition
Pietro Incardona's avatar
Pietro Incardona committed
361 362 363 364 365 366 367 368 369 370
	 *
	 * It makes a refinement of the current decomposition using Parmetis function RefineKWay
	 * After that it also does the remapping of the graph
	 *
	 */
	void refine()
	{
		size_t Np = v_cl.getProcessingUnits();
		size_t p_id = v_cl.getProcessUnitID();

371
		// Number of local vertex
372
		rid nl_vertex = vtxdist.get(p_id+1) - vtxdist.get(p_id);
373

Pietro Incardona's avatar
Pietro Incardona committed
374
		// Reset parmetis graph and reconstruct it
375
		parmetis_graph.reset(gp, vtxdist, m2g, verticesGotWeights);
Pietro Incardona's avatar
Pietro Incardona committed
376 377

		// Refine
378
		parmetis_graph.refine<nm_v::proc_id>(vtxdist);
Pietro Incardona's avatar
Pietro Incardona committed
379 380 381 382

		// Get result partition for this processor
		idx_t * partition = parmetis_graph.getPartition();

383 384
		partitions.get(p_id).resize(nl_vertex.id);
		std::copy(partition, partition + nl_vertex.id, &partitions.get(p_id).get(0));
Pietro Incardona's avatar
Pietro Incardona committed
385 386

		// Reset data structure to keep trace of new vertices distribution in processors (needed to update main graph)
387
		for (size_t i = 0; i < Np; ++i)
Pietro Incardona's avatar
Pietro Incardona committed
388 389 390 391
		{
			v_per_proc.get(i).clear();
		}

392 393
		openfpm::vector<size_t> prc;
		openfpm::vector<size_t> sz;
Pietro Incardona's avatar
Pietro Incardona committed
394 395 396 397 398 399 400 401
		openfpm::vector<void *> ptr;

		for (size_t i = 0; i < Np; i++)
		{
			if (i != v_cl.getProcessUnitID())
			{
				partitions.get(i).clear();
				prc.add(i);
402
				sz.add(nl_vertex.id * sizeof(idx_t));
Pietro Incardona's avatar
Pietro Incardona committed
403 404 405 406 407 408
				ptr.add(partitions.get(p_id).getPointer());
			}
		}

		// Exchange informations through processors
		v_cl.sendrecvMultipleMessagesNBX(prc.size(), &sz.get(0), &prc.get(0), &ptr.get(0), message_receive, &partitions,
409
		NONE);
Pietro Incardona's avatar
Pietro Incardona committed
410 411 412 413 414

		// Update graphs with the new distributions
		updateGraphs();
	}

415
	/*! \brief Compute the unbalance of the processor compared to the optimal balance
416
	 *
417
	 * \return the unbalance from the optimal one 0.01 mean 1%
418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436
	 */
	float getUnbalance()
	{
		long t_cost = 0;

		long min, max, sum;
		float unbalance;

		t_cost = getProcessorLoad();

		min = t_cost;
		max = t_cost;
		sum = t_cost;

		v_cl.min(min);
		v_cl.max(max);
		v_cl.sum(sum);
		v_cl.execute();

437
		unbalance = ((float) (max - min)) / (float) (sum / v_cl.getProcessingUnits());
438 439 440 441

		return unbalance * 100;
	}

442
	/*! \brief function that return the position of the vertex in the space
Pietro Incardona's avatar
Pietro Incardona committed
443 444 445 446 447
	 *
	 * \param id vertex id
	 * \param pos vector that will contain x, y, z
	 *
	 */
448
	void getSubSubDomainPosition(size_t id, T (&pos)[dim])
Pietro Incardona's avatar
Pietro Incardona committed
449
	{
450
		if (id >= gp.getNVertex())
451 452 453 454 455
			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";

		// Copy the geometrical informations inside the pos vector
		pos[0] = gp.vertex(id).template get<nm_v::x>()[0];
		pos[1] = gp.vertex(id).template get<nm_v::x>()[1];
Pietro Incardona's avatar
Pietro Incardona committed
456
		if (dim == 3)
457
			pos[2] = gp.vertex(id).template get<nm_v::x>()[2];
Pietro Incardona's avatar
Pietro Incardona committed
458 459
	}

460
	/*! \brief Function that set the weight of the vertex
Pietro Incardona's avatar
Pietro Incardona committed
461 462
	 *
	 * \param id vertex id
463
	 * \param weight to give to the vertex
Pietro Incardona's avatar
Pietro Incardona committed
464 465
	 *
	 */
466
	inline void setComputationCost(size_t id, size_t weight)
Pietro Incardona's avatar
Pietro Incardona committed
467
	{
468
		if (!verticesGotWeights)
469 470 471
			verticesGotWeights = true;

		if (id >= gp.getNVertex())
472 473 474
			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";

		// Update vertex in main graph
Pietro Incardona's avatar
Pietro Incardona committed
475 476 477
		gp.vertex(id).template get<nm_v::computation>() = weight;
	}

478 479
	/*! \brief Checks if weights are used on the vertices
	 *
480
	 * \return true if weights are used in the decomposition
481 482 483 484 485 486 487 488 489 490 491
	 */
	bool weightsAreUsed()
	{
		return verticesGotWeights;
	}

	/*! \brief function that get the weight of the vertex
	 *
	 * \param id vertex id
	 *
	 */
492
	size_t getSubSubDomainComputationCost(size_t id)
493 494
	{
		if (id >= gp.getNVertex())
495
			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
496 497 498 499

		return gp.vertex(id).template get<nm_v::computation>();
	}

500 501 502 503 504 505 506 507
	/*! \brief Compute the processor load counting the total weights of its vertices
	 *
	 * \return the computational load of the processor graph
	 */
	size_t getProcessorLoad()
	{
		size_t load = 0;

508 509 510 511
		// Processor id
		size_t p_id = v_cl.getProcessUnitID();


512
		for (rid i = vtxdist.get(p_id); i < vtxdist.get(p_id+1) ; ++i)
513
		{
514
			load += gp.vertex(m2g.find(i)->second.id).template get<nm_v::computation>();
515
		}
516
		//std::cout << v_cl.getProcessUnitID() << " weight " << load << " size " << sub_g.getNVertex() << "\n";
517 518 519
		return load;
	}

520 521 522 523 524 525 526 527
	/*! \brief Set migration cost of the vertex id
	 *
	 * \param id of the vertex to update
	 * \param migration cost of the migration
	 */
	void setMigrationCost(size_t id, size_t migration)
	{
		if (id >= gp.getNVertex())
528 529
			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";

530 531 532 533 534
		gp.vertex(id).template get<nm_v::migration>() = migration;
	}

	/*! \brief Set communication cost of the edge id
	 *
535 536 537
	 * \param v_id Id of the source vertex of the edge
	 * \param e i child of the vertex
	 * \param communication Communication value
538
	 */
539
	void setCommunicationCost(size_t v_id, size_t e, size_t communication)
540
	{
541 542 543 544
		size_t e_id = v_id + e;

		if (e_id >= gp.getNEdge())
			std::cerr << "Such edge doesn't exist (id = " << e_id << ", " << "total size = " << gp.getNEdge() << ")\n";
545

546
		gp.getChildEdge(v_id, e).template get<nm_e::communication>() = communication;
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563
	}

	/*! \brief Returns total number of sub-sub-domains in the distribution graph
	 *
	 */
	size_t getNSubSubDomains()
	{
		return gp.getNVertex();
	}

	/*! \brief Returns total number of neighbors of the sub-sub-domain id
	 *
	 * \param i id of the sub-sub-domain
	 */
	size_t getNSubSubDomainNeighbors(size_t id)
	{
		if (id >= gp.getNVertex())
564
			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
565 566 567 568

		return gp.getNChilds(id);
	}

569
	/*! \brief Print the current distribution and save it to VTK file
570
	 *
571
	 * \param file filename
572 573
	 *
	 */
574
	void write(const std::string & file)
575
	{
576 577
		//f (v_cl.getProcessUnitID() == 0)
		//{
578
			VTKWriter<Graph_CSR<nm_v, nm_e>, VTK_GRAPH> gv2(gp);
579 580
			gv2.write(std::to_string(v_cl.getProcessUnitID()) + file);
		//}
Pietro Incardona's avatar
Pietro Incardona committed
581
	}
582 583 584 585 586 587 588 589 590 591 592 593 594 595 596

	const ParMetisDistribution<dim,T> & operator=(const ParMetisDistribution<dim,T> & dist)
	{
		v_cl = dist.v_cl;
		gr = dist.gr;
		domain = dist.domain;
		gp = dist.gp;
		vtxdist = dist.vtxdist;
		partitions = dist.partitions;
		v_per_proc = dist.v_per_proc;
		verticesGotWeights = dist.verticesGotWeights;

		return *this;
	}

597
	const ParMetisDistribution<dim,T> & operator=(ParMetisDistribution<dim,T> && dist)
598 599 600 601 602 603 604 605 606 607 608 609
	{
		v_cl = dist.v_cl;
		gr = dist.gr;
		domain = dist.domain;
		gp.swap(dist.gp);
		vtxdist.swap(dist.vtxdist);
		partitions.swap(dist.partitions);
		v_per_proc.swap(dist.v_per_proc);
		verticesGotWeights = dist.verticesGotWeights;

		return *this;
	}
Pietro Incardona's avatar
Pietro Incardona committed
610 611 612
};

#endif /* SRC_DECOMPOSITION_PARMETISDISTRIBUTION_HPP_ */