ParMetisPartitioner.cc 13.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include <queue>
14
15
#include <mpi.h>

16
17
#include "parallel/ParMetisPartitioner.h"
#include "parallel/MpiHelper.h"
18
#include "Serializer.h"
19
20
21
22
23
24
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "FixVec.h"
#include "DOFVector.h"
25

26
27
28

namespace AMDiS {

29
  ParMetisMesh::ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm, 
30
			     std::map<int, bool>& elementInRank,
31
			     map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal)
Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
    : dim(mesh->getDim()),
      nElements(0),
34
      mpiComm(comm)
35
36
  {
    FUNCNAME("ParMetisMesh::ParMetisMesh()");
37
38

    int mpiSize = mpiComm->Get_size();
39
40
41
42
    int elementCounter = 0;
    int dow = Global::getGeo(WORLD);

    TraverseStack stack;
43
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
44
    while (elInfo) {
45
      if (elementInRank[elInfo->getElement()->getIndex()])
46
	elementCounter++;      
47
48
49
50

      elInfo = stack.traverseNext(elInfo);
    }

51
    nElements = elementCounter;
52

53
    TEST_EXIT(nElements > 0)("No elements in ParMETIS mesh!\n");
54
55

    // allocate memory
Thomas Witkowski's avatar
Thomas Witkowski committed
56
    eptr = new int[nElements + 1];
57
    eind = new int[nElements * (mesh->getGeo(VERTEX))];
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
    elmdist = new int[mpiSize + 1];
    elem_p2a = new int[nElements];
60

61
    if (dim == dow)
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      xyz = new float[nElements * dim];
63
64
    else
      xyz = NULL;    
65

Thomas Witkowski's avatar
Thomas Witkowski committed
66
    eptr[0] = 0;
67

Thomas Witkowski's avatar
Thomas Witkowski committed
68
69
70
    int *ptr_eptr = eptr + 1;
    int *ptr_eind = eind;
    float *ptr_xyz = xyz;
71
72
    
    // gather element numbers and create elmdist
73
    mpiComm->Allgather(&nElements, 1, MPI_INT, elmdist + 1, 1, MPI_INT);
74

75
    elmdist[0] = 0;
76
    for (int i = 2; i < mpiSize + 1; i++)
77
      elmdist[i] += elmdist[i - 1];
78
79

    // traverse mesh and fill distributed ParMETIS data
80
    DimVec<double> bary(dim, DEFAULT_VALUE, 1.0 / mesh->getGeo(VERTEX));
81
82
83
    WorldVector<double> world;

    elementCounter = 0;
84
    int nodeCounter = 0;
85

86
    elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL | Mesh::FILL_COORDS);
87
    while (elInfo) {
88
89
90
91
      Element *element = elInfo->getElement();
      int index = element->getIndex();

      // if element in partition
92
      if (elementInRank[index]) {
93
94
95
96
97
	// remember index
	setParMetisIndex(index, elementCounter);
	setAMDiSIndex(elementCounter, index);

	// write eptr entry
98
	nodeCounter += mesh->getGeo(VERTEX);
99
100
101
102
	*ptr_eptr = nodeCounter;
	ptr_eptr++;

	// write eind entries (element nodes)
Thomas Witkowski's avatar
Thomas Witkowski committed
103
	for (int i = 0; i < dim + 1; i++) {
104
	  if (mapLocalGlobal)
105
	    *ptr_eind = (*mapLocalGlobal)[element->getDof(i, 0)];
106
107
	  else
	    *ptr_eind = element->getDof(i, 0);	  
108

109
110
111
112
	  ptr_eind++;
	}

	// write xyz element coordinates
113
	if (ptr_xyz) {
114
	  elInfo->coordToWorld(bary, world);
Thomas Witkowski's avatar
Thomas Witkowski committed
115
	  for (int i = 0; i < dim; i++) {
116
	    *ptr_xyz = static_cast<float>(world[i]);
117
118
119
120
121
122
	    ptr_xyz++;
	  }
	}

	elementCounter++;
      }
123

124
125
126
127
      elInfo = stack.traverseNext(elInfo);
    }
  }

128

129
130
  ParMetisMesh::~ParMetisMesh()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
131
    if (eptr)
Thomas Witkowski's avatar
Thomas Witkowski committed
132
      delete [] eptr;
133
    
Thomas Witkowski's avatar
Thomas Witkowski committed
134
    if (eind)     
Thomas Witkowski's avatar
Thomas Witkowski committed
135
      delete [] eind;
136
    
137
    if (elmdist)
Thomas Witkowski's avatar
Thomas Witkowski committed
138
      delete [] elmdist;
139
    
Thomas Witkowski's avatar
Thomas Witkowski committed
140
    if (xyz)
Thomas Witkowski's avatar
Thomas Witkowski committed
141
      delete [] xyz;
142
    
Thomas Witkowski's avatar
Thomas Witkowski committed
143
    if (elem_p2a) 
Thomas Witkowski's avatar
Thomas Witkowski committed
144
      delete [] elem_p2a;
145
146
  }

147

148
  ParMetisGraph::ParMetisGraph(ParMetisMesh *parMesh,
149
			       MPI::Intracomm *comm,
150
			       int ncommonnodes)
151
    : parMetisMesh(parMesh)
152
  {
153
154
155
156
157
    FUNCNAME("ParMetisGraph::ParMetisGraph()");

    TEST_EXIT(parMesh)("No ParMetisMesh defined!\n");
    TEST_EXIT(comm)("No MPI communicator defined!\n");

158
159
    int numflag = 0;

160
161
    if (ncommonnodes == -1) 
      ncommonnodes = parMetisMesh->getDim();
162

163
    MPI_Comm tmpComm = MPI_Comm(*comm);
164

165
166
167
    ParMETIS_V3_Mesh2Dual(parMetisMesh->getElementDist(),
			  parMetisMesh->getElementPtr(),
			  parMetisMesh->getElementInd(),
168
169
			  &numflag,
			  &ncommonnodes,
Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
			  &xadj,
			  &adjncy,
172
			  &tmpComm);
173
174
  }

175

176
177
  ParMetisGraph::~ParMetisGraph()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
    free(xadj);
    free(adjncy);
180
181
  }

182

183
184
185
186
  void ParMetisGraph::print()
  {
    FUNCNAME("ParMetisGraph::print()");

187
    stringstream oss;
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
    for (int i = 0; i <= MPI::COMM_WORLD.Get_size(); i++)
      oss << parMetisMesh->getElementDist()[i] << " ";
    
    MSG("Element dist = %s\n", oss.str().c_str());

    int mpiRank = MPI::COMM_WORLD.Get_rank();
    int nElements = parMetisMesh->getElementDist()[mpiRank + 1] - 
      parMetisMesh->getElementDist()[mpiRank];

    MSG("nElements = %d   in index range %d - %d\n", 
	nElements, 
	parMetisMesh->getElementDist()[mpiRank],
	parMetisMesh->getElementDist()[mpiRank + 1]);

    oss.str("");
    oss.clear();
    
    for (int i = 0; i <= nElements; i++)
      oss << xadj[i] << ", ";

    MSG("xadj = {%s}\n", oss.str().c_str());

    oss.str("");
    oss.clear();

    for (int i = 0; i <= xadj[nElements] - 1; i++) 
      oss << adjncy[i] << ", ";

    MSG("adjncy = {%s}\n", oss.str().c_str());
  }


  ParMetisPartitioner::~ParMetisPartitioner()
  {
    if (parMetisMesh) 
      delete parMetisMesh;
  }


227
  bool ParMetisPartitioner::partition(map<int, double> &elemWeights,
228
				      PartitionMode mode) 
229
  {
230
    FUNCNAME("ParMetisPartitioner::partition()");
231

232
    int mpiSize = mpiComm->Get_size();
233

234
235
236

    // === Create parmetis mesh ===

237
    if (parMetisMesh) 
Thomas Witkowski's avatar
Thomas Witkowski committed
238
      delete parMetisMesh;
239

240
241
    TEST_EXIT_DBG(elementInRank.size() != 0)("Should not happen!\n");

242
    parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
243

244
    int nElements = parMetisMesh->getNumElements();
245

246
247
248

    // === Create weight array ===

249
250
    vector<int> wgts(nElements);
    vector<float> floatWgts(nElements);
251
    unsigned int floatWgtsPos = 0;
252
253
    float maxWgt = 0.0;

254
    TraverseStack stack;
255
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
256
    while (elInfo) {
257
      int index = elInfo->getElement()->getIndex();
258

259
      if (elementInRank[index]) {
260
	// get weight 
261
	float wgt = static_cast<float>(elemWeights[index]);
262
	maxWgt = std::max(wgt, maxWgt);
263
264

	// write float weight
265
266
	TEST_EXIT_DBG(floatWgtsPos < floatWgts.size())("Should not happen!\n");
	floatWgts[floatWgtsPos++] = wgt;
267
268
269
270
      }
      elInfo = stack.traverseNext(elInfo);
    }

271
272
    TEST_EXIT_DBG(floatWgtsPos == floatWgts.size())("Should not happen!\n");

273
    float tmp;
274
    mpiComm->Allreduce(&maxWgt, &tmp, 1, MPI_FLOAT, MPI_MAX);
275
276
    maxWgt = tmp;

277
278
279

    // === Create dual graph ===

280
    ParMetisGraph parMetisGraph(parMetisMesh, mpiComm);
281

282
283
284

    // === Partitioning of dual graph ===

285
    int wgtflag = 2; // weights at vertices only!
286
    int numflag = 0; // c numbering style!
287
    int ncon = 1; // one weight at each vertex!
288
    int nparts = mpiSize; // number of partitions
289

290
    vector<float> tpwgts(mpiSize);
291
    float ubvec = 1.05;
292
    int options[4] = {0, 0, 15, 1}; // default options
293
    int edgecut = -1;
294
    vector<int> part(nElements);
295

296
297
298
    // set tpwgts
    for (int i = 0; i < mpiSize; i++)
      tpwgts[i] = 1.0 / nparts;
299
   
300
    float scale = 10000.0 / maxWgt;
301
    for (int i = 0; i < nElements; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
302
      wgts[i] = static_cast<int>(floatWgts[i] * scale);
303
304
305


    // === Start ParMETIS. ===
306

307
308
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

309
    switch (mode) {
310
    case INITIAL:
311
      ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
312
313
			   parMetisGraph.getXAdj(),
			   parMetisGraph.getAdjncy(),
314
			   &(wgts[0]),
315
			   NULL,
316
			   &wgtflag,
317
318
319
			   &numflag,
			   &ncon,
			   &nparts,
320
			   &(tpwgts[0]),
321
322
323
			   &ubvec,
			   options,
			   &edgecut,
324
			   &(part[0]),
325
			   &tmpComm);
326
327
328
      break;
    case ADAPTIVE_REPART:
      {
329
	vector<int> vsize(nElements);
Thomas Witkowski's avatar
Thomas Witkowski committed
330
	for (int i = 0; i < nElements; i++)
331
332
	  vsize[i] = static_cast<int>(floatWgts[i]);

333
	ParMETIS_V3_AdaptiveRepart(parMetisMesh->getElementDist(),
334
335
				   parMetisGraph.getXAdj(),
				   parMetisGraph.getAdjncy(),
336
				   &(wgts[0]),
337
				   NULL,
338
				   &(vsize[0]),
339
340
341
342
				   &wgtflag,
				   &numflag,
				   &ncon,
				   &nparts,
343
				   &(tpwgts[0]),
344
345
346
347
				   &ubvec,
				   &itr,
				   options,
				   &edgecut,
348
				   &(part[0]),
349
				   &tmpComm);
350
351
352
      }
      break;
    case REFINE_PART:
353
      ParMETIS_V3_RefineKway(parMetisMesh->getElementDist(),
354
355
			     parMetisGraph.getXAdj(),
			     parMetisGraph.getAdjncy(),
356
			     &(wgts[0]),
357
358
359
360
361
			     NULL,
			     &wgtflag,
			     &numflag,
			     &ncon,
			     &nparts,
362
			     &(tpwgts[0]),
363
364
365
			     &ubvec,
			     options,
			     &edgecut,
366
			     &(part[0]),
367
			     &tmpComm);
368

369
370
371
372
373
      break;
    default: 
      ERROR_EXIT("unknown partitioning mode\n");
    }

374
375
376

    // === Distribute new partition data. ===

377
    return distributePartitioning(&(part[0]));
378
379
  }

380

381
  void ParMetisPartitioner::createPartitionMap(map<int, int>& pMap)
382
  {
383
    FUNCNAME("ParMetisPartitioner::createPartitionMap()");
384

385
    partitionMap.clear();
386
387

    // update ParMETIS mesh to new partitioning
388
    if (!parMetisMesh)
389
      parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
390

391
392
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
393
    vector<int> nPartitionElements(mpiSize);
394
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
395

396
    for (int i = 0; i < mpiSize; i++)
397
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
398
399

    // === count number of elements ===
400
401
402
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
403

404
    vector<int> partitionElements(nElements);
405
406

    // distribute partition elements
407
408
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
409
			MPI_INT, 
410
411
			&(partitionElements[0]), 
			&(nPartitionElements[0]), 
412
413
			elmdist, 
			MPI_INT);
414

415
    // fill partitionMap
Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
    for (int i = 0; i < mpiSize; i++)
      for (int j = 0; j < nPartitionElements[i]; j++)
418
	partitionMap[partitionElements[elmdist[i] + j]] = i;
419
420

    pMap = partitionMap;
421
422
  }

423

424
  bool ParMetisPartitioner::distributePartitioning(int *part) 
425
  {
426
427
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

428
429
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
430
    int nElements = parMetisMesh->getNumElements();
431

432
    // nPartitionElements[i] is the number of elements for the i-th partition
Thomas Witkowski's avatar
Thomas Witkowski committed
433
    int *nPartitionElements = new int[mpiSize];
434
    for (int i = 0; i < mpiSize; i++) 
435
      nPartitionElements[i] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
436
    for (int i = 0; i < nElements; i++)
437
      nPartitionElements[part[i]]++;    
438
439

    // collect number of partition elements from all ranks for this rank
Thomas Witkowski's avatar
Thomas Witkowski committed
440
    int *nRankElements = new int[mpiSize];
441
    mpiComm->Alltoall(nPartitionElements, 1, MPI_INT, nRankElements, 1, MPI_INT);
442

443

444
    // sum up partition elements over all ranks
Thomas Witkowski's avatar
Thomas Witkowski committed
445
    int *sumPartitionElements = new int[mpiSize];
446
    mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize,
447
		       MPI_INT, MPI_SUM);
448
449
450
451
452
453
454
455
456

    // Test if there exists an empty partition
    bool emptyPartition = false;
    for (int i = 0; i < mpiSize; i++) 
      emptyPartition |= (sumPartitionElements[i] == 0);

    if (emptyPartition)
      return false;    

457
    // prepare distribution (fill partitionElements with AMDiS indices)
Thomas Witkowski's avatar
Thomas Witkowski committed
458
    int *bufferOffset = new int[mpiSize];
459
    bufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
460
    for (int i = 1; i < mpiSize; i++)
461
      bufferOffset[i] = bufferOffset[i - 1] + nPartitionElements[i - 1];
462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
465

466
    for (int i = 0; i < mpiSize; i++)
467
468
      partitionPtr[i] = partitionElements + bufferOffset[i];

469
    sendElements.clear();
470
    for (int i = 0; i < nElements; i++) {
471
      int partition = part[i];
472
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
473

474
475
476
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

477
478
479
480
481
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

    // all to all: partition elements to rank elements
Thomas Witkowski's avatar
Thomas Witkowski committed
482
483
    int *rankElements = new int[sumPartitionElements[mpiRank]];
    int *recvBufferOffset = new int[mpiSize];
484
    recvBufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
485
    for (int i = 1; i < mpiSize; i++)
486
      recvBufferOffset[i] = recvBufferOffset[i - 1] + nRankElements[i - 1];
487

488
    mpiComm->Alltoallv(partitionElements, 
489
		       nPartitionElements,
490
491
492
		       bufferOffset,
		       MPI_INT,
		       rankElements,
493
		       nRankElements,
494
495
		       recvBufferOffset,
		       MPI_INT);
496
    
497
    TEST_EXIT(elementInRank.size() != 0)("Should not happen!\n");
498
    for (map<int, bool>::iterator it = elementInRank.begin();
499
500
	 it != elementInRank.end(); ++it)
      elementInRank[it->first] = false;
501

502
    // Create map which stores for each element index on macro level
503
    // if the element is in the partition of this rank.
504
    recvElements.clear();
505
    for (int i = 0; i < mpiSize; i++) {
506
      int *rankStart = rankElements + recvBufferOffset[i];
507
      int *rankEnd = rankStart + nRankElements[i];
508
509
      for (int *rankPtr = rankStart; rankPtr < rankEnd; ++rankPtr) {
	elementInRank[*rankPtr] = true;
510
511
	if (i != mpiRank)
	  recvElements[i].push_back(*rankPtr);
512
513
514
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
515
    delete parMetisMesh;
516
    parMetisMesh = NULL;
517

Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
520
521
522
523
524
    delete [] rankElements;
    delete [] nPartitionElements;
    delete [] nRankElements;
    delete [] sumPartitionElements;
    delete [] partitionElements;
    delete [] partitionPtr;
    delete [] bufferOffset;
Thomas Witkowski's avatar
huh    
Thomas Witkowski committed
525
    delete [] recvBufferOffset;
526
527

    return true;
528
  }
529
 
530
}