ParMetisPartitioner.cc 14.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 "parallel/ParMetisPartitioner.h"
#include "parallel/MpiHelper.h"
16
17
18
19
20
21
22
23
24
25
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "mpi.h"

namespace AMDiS {

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

    int mpiSize = mpiComm->Get_size();
36
37
38
39
    int elementCounter = 0;
    int dow = Global::getGeo(WORLD);

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

      elInfo = stack.traverseNext(elInfo);
    }

48
    nElements = elementCounter;
49

50
    TEST_EXIT(nElements > 0)("No elements in ParMETIS mesh!\n");
51
52

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
63
    eptr[0] = 0;
64

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

72
    elmdist[0] = 0;
73
    for (int i = 2; i < mpiSize + 1; i++)
74
      elmdist[i] += elmdist[i - 1];
75
76

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

    elementCounter = 0;
81
    int nodeCounter = 0;
82

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

      // if element in partition
89
      if (elementInRank[index]) {
90
91
92
93
94
	// remember index
	setParMetisIndex(index, elementCounter);
	setAMDiSIndex(elementCounter, index);

	// write eptr entry
95
	nodeCounter += mesh->getGeo(VERTEX);
96
97
98
99
	*ptr_eptr = nodeCounter;
	ptr_eptr++;

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

106
107
108
109
	  ptr_eind++;
	}

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

	elementCounter++;
      }
120

121
122
123
124
      elInfo = stack.traverseNext(elInfo);
    }
  }

125

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

144

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

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

155
156
    int numflag = 0;

157
158
    if (ncommonnodes == -1) 
      ncommonnodes = parMetisMesh->getDim();
159

160
    MPI_Comm tmpComm = MPI_Comm(*comm);
161

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

172

173
174
  ParMetisGraph::~ParMetisGraph()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
175
176
    free(xadj);
    free(adjncy);
177
178
  }

179

180
181
182
183
184
185
186
187
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
  void ParMetisGraph::print()
  {
    FUNCNAME("ParMetisGraph::print()");

    std::stringstream oss;
    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;
  }


224
225
  void ParMetisPartitioner::createPartitionData() 
  {
226
227
    FUNCNAME("ParMetrisPartitioner::createPartitionData()");

228
229
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
230
231
    int nLeaves = mesh->getNumberOfLeaves();
    int elPerRank = nLeaves / mpiSize;
232

233
    // === Create initial partitioning of the AMDiS mesh. ===
234

235
236
    elementInRank.clear();

237
    TraverseStack stack;
238
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
239
    while (elInfo) {
240
241
      Element *element = elInfo->getElement();

242
243
244
245
      if ((element->getIndex() >= mpiRank * elPerRank &&
	   element->getIndex() < (mpiRank + 1) * elPerRank) ||
	  (element->getIndex() >= mpiSize * elPerRank &&
	   mpiRank == mpiSize - 1))
246
	elementInRank[element->getIndex()] = true;
247
      else
248
	elementInRank[element->getIndex()] = false;
249
      
250
251
252
253
      elInfo = stack.traverseNext(elInfo);
    }
  }

254

255
  void ParMetisPartitioner::partition(std::map<int, double> &elemWeights,
256
257
258
				      PartitionMode mode,
				      float itr) 
  {
259
    FUNCNAME("ParMetisPartitioner::partition()");
260

261
    int mpiSize = mpiComm->Get_size();
262

263
264
265

    // === Create parmetis mesh ===

266
    if (parMetisMesh) 
Thomas Witkowski's avatar
Thomas Witkowski committed
267
      delete parMetisMesh;
268

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

271
    parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
272

273
    int nElements = parMetisMesh->getNumElements();
274

275
276
277

    // === Create weight array ===

278
279
280
    std::vector<int> wgts(nElements);
    std::vector<float> floatWgts(nElements);
    unsigned int floatWgtsPos = 0;
281
282
    float maxWgt = 0.0;

283
    TraverseStack stack;
284
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
285
    while (elInfo) {
286
      int index = elInfo->getElement()->getIndex();
287

288
      if (elementInRank[index]) {
289
	// get weight 
290
	float wgt = static_cast<float>(elemWeights[index]);
291
292
293
	maxWgt = max(wgt, maxWgt);

	// write float weight
294
295
	TEST_EXIT_DBG(floatWgtsPos < floatWgts.size())("Should not happen!\n");
	floatWgts[floatWgtsPos++] = wgt;
296
297
298
299
      }
      elInfo = stack.traverseNext(elInfo);
    }

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

302
    float tmp;
303
    mpiComm->Allreduce(&maxWgt, &tmp, 1, MPI_FLOAT, MPI_MAX);
304
305
    maxWgt = tmp;

306
307
308

    // === Create dual graph ===

309
    ParMetisGraph parMetisGraph(parMetisMesh, mpiComm);
310

311
312
313

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

314
    int wgtflag = 2; // weights at vertices only!
315
    int numflag = 0; // c numbering style!
316
    int ncon = 1; // one weight at each vertex!
317
    int nparts = mpiSize; // number of partitions
318

319
    std::vector<float> tpwgts(mpiSize);
320
    float ubvec = 1.05;
321
    int options[4] = {0, 0, 15, 1}; // default options
322
    int edgecut = -1;
323
    std::vector<int> part(nElements);
324

325
326
327
    // set tpwgts
    for (int i = 0; i < mpiSize; i++)
      tpwgts[i] = 1.0 / nparts;
328
   
329
    float scale = 10000.0 / maxWgt;
330
331
332
333
334
335
336


    // === Scale element weights. ===

    double weightSum = 0.0;

    for (int i = 0; i < nElements; i++) {
337
      wgts[i] = static_cast<int>(floatWgts[i] * scale);
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
      weightSum += wgts[i];
    }

    mpi::globalAdd(weightSum);
    weightSum /= mpiSize;

    for (int i = 0; i < nElements; i++) {
      if (wgts[i] > weightSum) {
	MSG("ICH HABE EINEN %d > %f\n", wgts[i], weightSum);
	wgts[i] = static_cast<int>(weightSum);
      }
    }


    // === Start ParMETIS. ===
353

354
355
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

356
    switch (mode) {
357
    case INITIAL:
358
      ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
359
360
			   parMetisGraph.getXAdj(),
			   parMetisGraph.getAdjncy(),
361
			   &(wgts[0]),
362
			   NULL,
363
			   &wgtflag,
364
365
366
			   &numflag,
			   &ncon,
			   &nparts,
367
			   &(tpwgts[0]),
368
369
370
			   &ubvec,
			   options,
			   &edgecut,
371
			   &(part[0]),
372
			   &tmpComm);
373
374
375
      break;
    case ADAPTIVE_REPART:
      {
376
	std::vector<int> vsize(nElements);
Thomas Witkowski's avatar
Thomas Witkowski committed
377
	for (int i = 0; i < nElements; i++)
378
379
	  vsize[i] = static_cast<int>(floatWgts[i]);

380
	ParMETIS_V3_AdaptiveRepart(parMetisMesh->getElementDist(),
381
382
				   parMetisGraph.getXAdj(),
				   parMetisGraph.getAdjncy(),
383
				   &(wgts[0]),
384
				   NULL,
385
				   &(vsize[0]),
386
387
388
389
				   &wgtflag,
				   &numflag,
				   &ncon,
				   &nparts,
390
				   &(tpwgts[0]),
391
392
393
394
				   &ubvec,
				   &itr,
				   options,
				   &edgecut,
395
				   &(part[0]),
396
				   &tmpComm);
397
398
399
      }
      break;
    case REFINE_PART:
400
      ParMETIS_V3_RefineKway(parMetisMesh->getElementDist(),
401
402
			     parMetisGraph.getXAdj(),
			     parMetisGraph.getAdjncy(),
403
			     &(wgts[0]),
404
405
406
407
408
			     NULL,
			     &wgtflag,
			     &numflag,
			     &ncon,
			     &nparts,
409
			     &(tpwgts[0]),
410
411
412
			     &ubvec,
			     options,
			     &edgecut,
413
			     &(part[0]),
414
			     &tmpComm);
415

416
417
418
419
420
      break;
    default: 
      ERROR_EXIT("unknown partitioning mode\n");
    }

421
422
423

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

424
    distributePartitioning(&(part[0]));
425
426
  }

427

428
429
  void ParMetisPartitioner::fillCoarsePartitionVec(std::map<int, int> *partitionVec)
  {
430
431
    FUNCNAME("ParMetisPartitioner::fillCoarsePartitionVec()");

432
    TEST_EXIT_DBG(partitionVec)("No partition vector!\n");
433
434
435
436

    partitionVec->clear();

    // update ParMETIS mesh to new partitioning
437
    if (!parMetisMesh)
438
      parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
439

440
441
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
442
    std::vector<int> nPartitionElements(mpiSize);
443
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
444

445
    for (int i = 0; i < mpiSize; i++)
446
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
447
448

    // === count number of elements ===
449
450
451
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
452

453
    std::vector<int> partitionElements(nElements);
454
455

    // distribute partition elements
456
457
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
458
			MPI_INT, 
459
460
			&(partitionElements[0]), 
			&(nPartitionElements[0]), 
461
462
			elmdist, 
			MPI_INT);
463
464

    // fill partitionVec
Thomas Witkowski's avatar
Thomas Witkowski committed
465
466
    for (int i = 0; i < mpiSize; i++)
      for (int j = 0; j < nPartitionElements[i]; j++)
467
468
469
	(*partitionVec)[partitionElements[elmdist[i] + j]] = i;
  }

470

471
472
  void ParMetisPartitioner::distributePartitioning(int *part) 
  {
473
474
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

475
476
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
477
    int nElements = parMetisMesh->getNumElements();
478

479
    // nPartitionElements[i] is the number of elements for the i-th partition
Thomas Witkowski's avatar
Thomas Witkowski committed
480
    int *nPartitionElements = new int[mpiSize];
481
    for (int i = 0; i < mpiSize; i++) 
482
      nPartitionElements[i] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    for (int i = 0; i < nElements; i++)
484
      nPartitionElements[part[i]]++;    
485
486

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

490

491
    // sum up partition elements over all ranks
Thomas Witkowski's avatar
Thomas Witkowski committed
492
    int *sumPartitionElements = new int[mpiSize];
493
    mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize,
494
		       MPI_INT, MPI_SUM);
495
496
    
    // prepare distribution (fill partitionElements with AMDiS indices)
Thomas Witkowski's avatar
Thomas Witkowski committed
497
    int *bufferOffset = new int[mpiSize];
498
    bufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
499
    for (int i = 1; i < mpiSize; i++)
500
      bufferOffset[i] = bufferOffset[i - 1] + nPartitionElements[i - 1];
501

Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
504

505
    for (int i = 0; i < mpiSize; i++)
506
507
      partitionPtr[i] = partitionElements + bufferOffset[i];

508
    sendElements.clear();
509
    for (int i = 0; i < nElements; i++) {
510
      int partition = part[i];
511
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
512

513
514
515
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

516
517
518
519
520
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

    // all to all: partition elements to rank elements
Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
    int *rankElements = new int[sumPartitionElements[mpiRank]];
    int *recvBufferOffset = new int[mpiSize];
523
    recvBufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
524
    for (int i = 1; i < mpiSize; i++)
525
      recvBufferOffset[i] = recvBufferOffset[i - 1] + nRankElements[i - 1];
526

527
    mpiComm->Alltoallv(partitionElements, 
528
		       nPartitionElements,
529
530
531
		       bufferOffset,
		       MPI_INT,
		       rankElements,
532
		       nRankElements,
533
534
		       recvBufferOffset,
		       MPI_INT);
535
    
536
537
538
539
    TEST_EXIT(elementInRank.size() != 0)("Should not happen!\n");
    for (std::map<int, bool>::iterator it = elementInRank.begin();
	 it != elementInRank.end(); ++it)
      elementInRank[it->first] = false;
540

541
542
    // Create map which stores for each element index on ther partitioning level
    // if the element is in the partition of this rank.
543
    recvElements.clear();
544
    for (int i = 0; i < mpiSize; i++) {
545
      int *rankStart = rankElements + recvBufferOffset[i];
546
      int *rankEnd = rankStart + nRankElements[i];
547
548
      for (int *rankPtr = rankStart; rankPtr < rankEnd; ++rankPtr) {
	elementInRank[*rankPtr] = true;
549
550
	if (i != mpiRank)
	  recvElements[i].push_back(*rankPtr);
551
552
553
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
554
    delete parMetisMesh;
555
    parMetisMesh = NULL;
556

Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
559
560
561
562
563
    delete [] rankElements;
    delete [] nPartitionElements;
    delete [] nRankElements;
    delete [] sumPartitionElements;
    delete [] partitionElements;
    delete [] partitionPtr;
    delete [] bufferOffset;
Thomas Witkowski's avatar
huh    
Thomas Witkowski committed
564
    delete [] recvBufferOffset;
565
566
567
  }

}