ParMetisPartitioner.cc 15.5 KB
Newer Older
1
#include <queue>
2
3
4
5
6
7
8
9
10
11
12
13
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "FixVec.h"
#include "PartitionElementData.h"
#include "DOFVector.h"
#include "mpi.h"

namespace AMDiS {

14
15
  ParMetisMesh::ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm, 
			     std::map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal)
Thomas Witkowski's avatar
Thomas Witkowski committed
16
17
    : dim(mesh->getDim()),
      nElements(0),
18
      mpiComm(comm)
19
20
  {
    FUNCNAME("ParMetisMesh::ParMetisMesh()");
21
22

    int mpiSize = mpiComm->Get_size();
23
24
25
26
27
    int nodeCounter = 0;
    int elementCounter = 0;
    int dow = Global::getGeo(WORLD);

    TraverseStack stack;
28
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
29
    while (elInfo) {
30
31
      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
Thomas Witkowski's avatar
Thomas Witkowski committed
32
	(elInfo->getElement()->getElementData(PARTITION_ED));
33

34
35
      if (partitionData && 
	  partitionData->getPartitionStatus() == IN &&
36
37
	  partitionData->getLevel() == 0)
	elementCounter++;      
38
39
40
41

      elInfo = stack.traverseNext(elInfo);
    }

42
    nElements = elementCounter;
43

44
    TEST_EXIT(nElements > 0)("No elements in ParMETIS mesh!\n");
45
46

    // allocate memory
Thomas Witkowski's avatar
Thomas Witkowski committed
47
    eptr = new int[nElements + 1];
Thomas Witkowski's avatar
huh    
Thomas Witkowski committed
48
    eind = new int[nElements * (dim + 1)];
Thomas Witkowski's avatar
Thomas Witkowski committed
49
50
    elmdist = new int[mpiSize + 1];
    elem_p2a = new int[nElements];
51

52
    if (dim == dow)
Thomas Witkowski's avatar
Thomas Witkowski committed
53
      xyz = new float[nElements * dim];
54
55
    else
      xyz = NULL;    
56

Thomas Witkowski's avatar
Thomas Witkowski committed
57
    eptr[0] = 0;
58

Thomas Witkowski's avatar
Thomas Witkowski committed
59
60
61
    int *ptr_eptr = eptr + 1;
    int *ptr_eind = eind;
    float *ptr_xyz = xyz;
62
63
    
    // gather element numbers and create elmdist
64
    mpiComm->Allgather(&nElements, 1, MPI_INT, elmdist + 1, 1, MPI_INT);
65

66
    elmdist[0] = 0;
67
    for (int i = 2; i < mpiSize + 1; i++)
68
      elmdist[i] += elmdist[i - 1];
69
70

    // traverse mesh and fill distributed ParMETIS data
Thomas Witkowski's avatar
Thomas Witkowski committed
71
    DimVec<double> bary(dim, DEFAULT_VALUE, 1.0 / (dim + 1));
72
73
74
75
76
    WorldVector<double> world;

    elementCounter = 0;

    elInfo = stack.traverseFirst(mesh, -1, 
77
				 Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_COORDS);
78
    while (elInfo) {
79
80
81
82
83
84
85
86
      Element *element = elInfo->getElement();
      int index = element->getIndex();

      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

      // if element in partition
87
88
89
      if (partitionData && 
	  partitionData->getPartitionStatus() == IN &&
	  partitionData->getLevel() == 0) {
90
91
92
93
94
	// remember index
	setParMetisIndex(index, elementCounter);
	setAMDiSIndex(elementCounter, index);

	// write eptr entry
Thomas Witkowski's avatar
Thomas Witkowski committed
95
	nodeCounter += dim + 1;
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
102
103
104
105
	  if (mapLocalGlobal)
	    *ptr_eind = (*mapLocalGlobal)[element->getDOF(i, 0)];
	  else
	    *ptr_eind = element->getDOF(i, 0);

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
114
115
116
117
118
119
120
121
122
123
	    *ptr_xyz = static_cast<float>(world[i]); 
	    ptr_xyz++;
	  }
	}

	elementCounter++;
      }
      elInfo = stack.traverseNext(elInfo);
    }
  }

124

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

143

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

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

154
155
    int numflag = 0;

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

159
    MPI_Comm tmpComm = MPI_Comm(*comm);
160
161
    int mpiRank = MPI::COMM_WORLD.Get_rank();
    int mpiSize = MPI::COMM_WORLD.Get_size();
162

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

173

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

180

181
182
183
  void ParMetisPartitioner::deletePartitionData() 
  {
    TraverseStack stack;
184
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
185
    while (elInfo) {
186
187
188
189
190
191
      Element *element = elInfo->getElement();
      element->deleteElementData(PARTITION_ED);
      elInfo = stack.traverseNext(elInfo);
    }
  }

192

193
194
  void ParMetisPartitioner::createPartitionData() 
  {
195
196
    FUNCNAME("ParMetrisPartitioner::createPartitionData()");

197
198
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
199
200
    int nLeaves = mesh->getNumberOfLeaves();
    int elPerRank = nLeaves / mpiSize;
201

202
    // === Create initial partitioning of the AMDiS mesh. ===
203

204
205
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
206
    while (elInfo) {
207
208
209
210
211
212
      Element *element = elInfo->getElement();

      TEST_EXIT(element->getElementData(PARTITION_ED) == NULL)
	("mesh already partitioned\n");

      PartitionElementData *elData = 
Thomas Witkowski's avatar
Thomas Witkowski committed
213
	new PartitionElementData(element->getElementData());
214
215
      element->setElementData(elData);

216
217
218
219
      if ((element->getIndex() >= mpiRank * elPerRank &&
	   element->getIndex() < (mpiRank + 1) * elPerRank) ||
	  (element->getIndex() >= mpiSize * elPerRank &&
	   mpiRank == mpiSize - 1))
220
	elData->setPartitionStatus(IN);
221
      else
222
	elData->setPartitionStatus(UNDEFINED);
223
      
224
225
226
227
      elInfo = stack.traverseNext(elInfo);
    }
  }

228

229
230
231
232
  void ParMetisPartitioner::partition(std::map<int, double> *elemWeights,
				      PartitionMode mode,
				      float itr) 
  {
233
    FUNCNAME("ParMetisPartitioner::partition()");
234

235
    int mpiSize = mpiComm->Get_size();
236
237

    // === create parmetis mesh ===
238
    if (parMetisMesh) 
Thomas Witkowski's avatar
Thomas Witkowski committed
239
      delete parMetisMesh;
240

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

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

    // === create weight array ===
Thomas Witkowski's avatar
Thomas Witkowski committed
246
247
    int *wgts = elemWeights ? new int[nElements] : NULL;
    float *floatWgts = elemWeights ? new float[nElements] : NULL;
248
249
250
    float maxWgt = 0.0;
    float *ptr_floatWgts = floatWgts;

251
252
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
253
    while (elInfo) {
254
255
256
257
258
259
      Element *element = elInfo->getElement();

      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

260
261
262
      if (partitionData &&
	  partitionData->getPartitionStatus() == IN &&
	  partitionData->getLevel() == 0) {
263

264
265
266
267
268
269
270
271
272
273
274
275
276
277
	int index = element->getIndex();

	// get weight 
	float wgt = static_cast<float>((*elemWeights)[index]);
	maxWgt = max(wgt, maxWgt);

	// write float weight
	*ptr_floatWgts = wgt;
	ptr_floatWgts++;
      }
      elInfo = stack.traverseNext(elInfo);
    }

    float tmp;
278
    mpiComm->Allreduce(&maxWgt, &tmp, 1, MPI_FLOAT, MPI_MAX);
279
280
281
    maxWgt = tmp;

    // === create dual graph ===
282
    ParMetisGraph parMetisGraph(parMetisMesh, mpiComm);
283
284
285
286
287
288

    // === partitioning of dual graph ===
    int wgtflag = elemWeights ? 2 : 0; // weights at vertices only!
    int numflag = 0; // c numbering style!
    int ncon = elemWeights ? 1 : 0; // one weight at each vertex!
    int nparts = mpiSize; // number of partitions
289

Thomas Witkowski's avatar
Thomas Witkowski committed
290
    float *tpwgts = elemWeights ? new float[mpiSize] : NULL;
291
292
293
    float ubvec = 1.05;
    int options[3] = {0, 0, 0}; // default options
    int edgecut = -1;
Thomas Witkowski's avatar
Thomas Witkowski committed
294
    int *part = new int[nElements];
295

296
    if (elemWeights) {
Thomas Witkowski's avatar
Thomas Witkowski committed
297
298
      // set tpwgts
      for (int i = 0; i < mpiSize; i++)
299
	tpwgts[i] = 1.0 / nparts;
300

301
      float scale = 10000.0 / maxWgt;
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
      // scale wgts
      for (int i = 0; i < nElements; i++)
304
305
306
	wgts[i] = static_cast<int>(floatWgts[i] * scale);
    }

307
308
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

309
310
    switch(mode) {
    case INITIAL:
311
      ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
312
313
			   parMetisGraph.getXAdj(),
			   parMetisGraph.getAdjncy(),
314
			   wgts,
315
			   NULL,
316
			   &wgtflag,
317
318
319
320
321
322
323
324
			   &numflag,
			   &ncon,
			   &nparts,
			   tpwgts,
			   &ubvec,
			   options,
			   &edgecut,
			   part,
325
			   &tmpComm);
326
327
328
      break;
    case ADAPTIVE_REPART:
      {
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
	int *vsize = new int[nElements];
	for (int i = 0; i < nElements; i++)
331
	  vsize[i] = 1;
332
	ParMETIS_V3_AdaptiveRepart(parMetisMesh->getElementDist(),
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
				   parMetisGraph.getXAdj(),
				   parMetisGraph.getAdjncy(),
				   wgts,
				   NULL,
				   vsize,
				   &wgtflag,
				   &numflag,
				   &ncon,
				   &nparts,
				   tpwgts,
				   &ubvec,
				   &itr,
				   options,
				   &edgecut,
				   part,
348
				   &tmpComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
349
	delete [] vsize;
350
351
352
      }
      break;
    case REFINE_PART:
353
      ParMETIS_V3_RefineKway(parMetisMesh->getElementDist(),
354
355
356
357
358
359
360
361
362
363
364
365
366
			     parMetisGraph.getXAdj(),
			     parMetisGraph.getAdjncy(),
			     wgts,
			     NULL,
			     &wgtflag,
			     &numflag,
			     &ncon,
			     &nparts,
			     tpwgts,
			     &ubvec,
			     options,
			     &edgecut,
			     part,
367
			     &tmpComm);
368
369
370
371
372
373
374
375
      break;
    default: 
      ERROR_EXIT("unknown partitioning mode\n");
    }

    // === distribute new partition data ===
    distributePartitioning(part);

376
    if (floatWgts) 
Thomas Witkowski's avatar
Thomas Witkowski committed
377
      delete [] floatWgts;
378
    if (wgts) 
Thomas Witkowski's avatar
Thomas Witkowski committed
379
      delete [] wgts;
380
    if (tpwgts) 
Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
383
      delete [] tpwgts;

    delete [] part;
384
385
  }

386

387
388
389
390
391
392
393
  void ParMetisPartitioner::fillCoarsePartitionVec(std::map<int, int> *partitionVec)
  {
    TEST_EXIT(partitionVec)("no partition vector\n");

    partitionVec->clear();

    // update ParMETIS mesh to new partitioning
394
    if (!parMetisMesh)
395
      parMetisMesh = new ParMetisMesh(mesh, mpiComm, mapLocalGlobal);
396

397
398
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
399
    std::vector<int> nPartitionElements(mpiSize);
400
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
401

402
    for (int i = 0; i < mpiSize; i++)
403
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
404
405

    // === count number of elements ===
406
407
408
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
409

410
    std::vector<int> partitionElements(nElements);
411
412

    // distribute partition elements
413
414
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
415
			MPI_INT, 
416
417
			&(partitionElements[0]), 
			&(nPartitionElements[0]), 
418
419
			elmdist, 
			MPI_INT);
420
421

    // fill partitionVec
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
    for (int i = 0; i < mpiSize; i++)
      for (int j = 0; j < nPartitionElements[i]; j++)
424
425
426
	(*partitionVec)[partitionElements[elmdist[i] + j]] = i;
  }

427

428
429
  void ParMetisPartitioner::distributePartitioning(int *part) 
  {
430
431
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

432
433
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
434
    int nElements = parMetisMesh->getNumElements();
435

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

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

    // sum up partition elements over all ranks
Thomas Witkowski's avatar
Thomas Witkowski committed
448
    int *sumPartitionElements = new int[mpiSize];
449
    mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize,
450
		       MPI_INT, MPI_SUM);
451
452
    
    // prepare distribution (fill partitionElements with AMDiS indices)
Thomas Witkowski's avatar
Thomas Witkowski committed
453
    int *bufferOffset = new int[mpiSize];
454
    bufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
455
    for (int i = 1; i < mpiSize; i++)
456
      bufferOffset[i] = bufferOffset[i - 1] + nPartitionElements[i - 1];
457

Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
460

Thomas Witkowski's avatar
Thomas Witkowski committed
461
    for (int i = 0; i < mpiSize; i++)
462
463
      partitionPtr[i] = partitionElements + bufferOffset[i];

464
    for (int i = 0; i < nElements; i++) {
465
      int partition = part[i];
466
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
467
468
469
470
471
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

    // all to all: partition elements to rank elements
Thomas Witkowski's avatar
Thomas Witkowski committed
472
473
    int *rankElements = new int[sumPartitionElements[mpiRank]];
    int *recvBufferOffset = new int[mpiSize];
474
    recvBufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
475
    for (int i = 1; i < mpiSize; i++)
476
      recvBufferOffset[i] = recvBufferOffset[i - 1] + nRankElements[i - 1];
477

478
    mpiComm->Alltoallv(partitionElements, 
479
		       nPartitionElements,
480
481
482
		       bufferOffset,
		       MPI_INT,
		       rankElements,
483
		       nRankElements,
484
485
		       recvBufferOffset,
		       MPI_INT);
486
487
    

488
489
    // Create map which stores for each element index on ther partitioning level
    // if the element is in the partition of this rank.
490
    std::map<int, bool> elementInPartition;
491
    for (int i = 0; i < mpiSize; i++) {
492
      int *rankStart = rankElements + recvBufferOffset[i];
493
      int *rankEnd = rankStart + nRankElements[i];
494
      for (int *rankPtr = rankStart; rankPtr < rankEnd; ++rankPtr)
495
496
497
498
	elementInPartition[*rankPtr] = true;
    }

    TraverseStack stack;
499
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
500
    while (elInfo) {
501
502
503
504
505
506
      Element *element = elInfo->getElement();

      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

507
      if (partitionData && partitionData->getLevel() == 0) {
508
	int amdisIndex = element->getIndex();
509
	if (elementInPartition[amdisIndex])
510
	  partitionData->setPartitionStatus(IN);
511
	else
512
	  partitionData->setPartitionStatus(OUT);
513
	
514
515
516
517
518
519
	descendPartitionData(element);
      }

      elInfo = stack.traverseNext(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
520
    delete parMetisMesh;
521
    parMetisMesh = NULL;
522

Thomas Witkowski's avatar
Thomas Witkowski committed
523
524
525
526
527
528
529
    delete [] rankElements;
    delete [] nPartitionElements;
    delete [] nRankElements;
    delete [] sumPartitionElements;
    delete [] partitionElements;
    delete [] partitionPtr;
    delete [] bufferOffset;
Thomas Witkowski's avatar
huh    
Thomas Witkowski committed
530
    delete [] recvBufferOffset;
531
532
  }

533

534
535
  void ParMetisPartitioner::descendPartitionData(Element *element) 
  {
536
537
    FUNCNAME("ParMetisPartitioner::descendPartitionData()");

538
    if (!element->isLeaf()) {
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
      Element *child0 = element->getChild(0);
      Element *child1 = element->getChild(1);

      // get partition data
      PartitionElementData *parentData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));
      PartitionElementData *child0Data = dynamic_cast<PartitionElementData*>
	(child0->getElementData(PARTITION_ED));
      PartitionElementData *child1Data = dynamic_cast<PartitionElementData*>
	(child1->getElementData(PARTITION_ED));
      
      TEST_EXIT(parentData && child0Data && child1Data)("no partition data\n");

      child0Data->setPartitionStatus(parentData->getPartitionStatus());
      child1Data->setPartitionStatus(parentData->getPartitionStatus());

      descendPartitionData(child0);
      descendPartitionData(child1);
    }
  }


  void ParMetisPartitioner::fillLeafPartitionVec(std::map<int, int> *coarseVec,
						 std::map<int, int> *fineVec)
  {
    int partition = -1;
    TraverseStack stack;
566
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
567
    while (elInfo) {
568
569
570
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));
571
      if (partitionData) {
572
	if (partitionData->getLevel() == 0)
573
	  partition = (*(coarseVec))[element->getIndex()];
574
	if (element->isLeaf())
575
576
577
578
579
580
	  (*(fineVec))[element->getIndex()] = partition;
      }
      elInfo = stack.traverseNext(elInfo);
    }
  }
}