ParMetisPartitioner.cc 15.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#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"

#include <queue>

namespace AMDiS {

15
  ParMetisMesh::ParMetisMesh(Mesh *mesh, MPI::Comm *comm)
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
28
29
    int nodeCounter = 0;
    int elementCounter = 0;
    int dow = Global::getGeo(WORLD);

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

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

      elInfo = stack.traverseNext(elInfo);
    }

44
    nElements = elementCounter;
45

46
    TEST_EXIT(nElements > 0)("no elements in ParMETIS mesh\n");
47
48

    // allocate memory
Thomas Witkowski's avatar
Thomas Witkowski committed
49
50
    eptr = GET_MEMORY(int, nElements + 1);
    eind = GET_MEMORY(int, nElements * (dim + 1));
51
    elmdist = GET_MEMORY(int, mpiSize + 1);
52

Thomas Witkowski's avatar
Thomas Witkowski committed
53
    elem_p2a = GET_MEMORY(int, nElements);
54

Thomas Witkowski's avatar
Thomas Witkowski committed
55
56
    if (dim == dow) {
      xyz = GET_MEMORY(float, nElements * dim);
57
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
58
      xyz = NULL;
59
60
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
61
    eptr[0] = 0;
62

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

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

    // traverse mesh and fill distributed ParMETIS data
Thomas Witkowski's avatar
Thomas Witkowski committed
76
    DimVec<double> bary(dim, DEFAULT_VALUE, 1.0 / (dim + 1));
77
78
79
80
81
82
83
    WorldVector<double> world;

    elementCounter = 0;

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

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

      // if element in partition
93
94
95
      if (partitionData && 
	  partitionData->getPartitionStatus() == IN &&
	  partitionData->getLevel() == 0) {
96
97
98
99
100
	// remember index
	setParMetisIndex(index, elementCounter);
	setAMDiSIndex(elementCounter, index);

	// write eptr entry
Thomas Witkowski's avatar
Thomas Witkowski committed
101
	nodeCounter += dim + 1;
102
103
104
105
	*ptr_eptr = nodeCounter;
	ptr_eptr++;

	// write eind entries (element nodes)
Thomas Witkowski's avatar
Thomas Witkowski committed
106
	for (int i = 0; i < dim + 1; i++) {
107
108
109
110
111
	  *ptr_eind = element->getDOF(i, 0);
	  ptr_eind++;
	}

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

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

  ParMetisMesh::~ParMetisMesh()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
128
129
    if (eptr)
      FREE_MEMORY(eptr, int, nElements + 1);
130
    
Thomas Witkowski's avatar
Thomas Witkowski committed
131
132
    if (eind)     
      FREE_MEMORY(eind, int, nElements * (dim + 1));
133
    
134
135
    if (elmdist)
      FREE_MEMORY(elmdist, int, mpiComm->Get_size() + 1);
136
    
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
    if (xyz)
      FREE_MEMORY(xyz, float, nElements * dim);
139
    
Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
    if (elem_p2a) 
      FREE_MEMORY(elem_p2a, int, nElements);
142
143
  }

144
  ParMetisGraph::ParMetisGraph(ParMetisMesh *parMesh,
145
			       MPI::Comm *comm,
146
			       int ncommonnodes)
147
    : parMetisMesh(parMesh)
148
149
150
  {
    int numflag = 0;

151
152
    if (ncommonnodes == -1) 
      ncommonnodes = parMetisMesh->getDim();
153

154
    MPI_Comm tmpComm = MPI_Comm(*comm);
155

156
157
158
    ParMETIS_V3_Mesh2Dual(parMetisMesh->getElementDist(),
			  parMetisMesh->getElementPtr(),
			  parMetisMesh->getElementInd(),
159
160
			  &numflag,
			  &ncommonnodes,
Thomas Witkowski's avatar
Thomas Witkowski committed
161
162
			  &xadj,
			  &adjncy,
163
			  &tmpComm);
164
165
166
167
  }

  ParMetisGraph::~ParMetisGraph()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
169
    free(xadj);
    free(adjncy);
170
171
  }

172

173
174
175
176
177
  void ParMetisPartitioner::deletePartitionData() 
  {
    TraverseStack stack;
    ElInfo *elInfo;
    elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER);
178
    while (elInfo) {
179
180
181
182
183
184
185
      Element *element = elInfo->getElement();
      element->deleteElementData(PARTITION_ED);
      elInfo = stack.traverseNext(elInfo);
    }
  }

  void ParMetisPartitioner::createPartitionData() {
186
187
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
188
189
190
191
192
193
194

    TraverseStack stack;
    ElInfo *elInfo;

    // === create initial partitioning on AMDiS mesh ===
    int totalElementCounter = 0;
    elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_LEAF_EL);
195
    while (elInfo) {
196
197
198
199
200
201
202
203
204
      Element *element = elInfo->getElement();

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

      PartitionElementData *elData = 
	NEW PartitionElementData(element->getElementData());
      element->setElementData(elData);

205
      if (totalElementCounter % mpiSize == mpiRank) {
206
207
208
209
210
211
212
213
214
215
216
217
218
219
	elData->setPartitionStatus(IN);
      } else {
	elData->setPartitionStatus(UNDEFINED);
      }
      totalElementCounter++;

      elInfo = stack.traverseNext(elInfo);
    }
  }

  void ParMetisPartitioner::partition(std::map<int, double> *elemWeights,
				      PartitionMode mode,
				      float itr) 
  {
220
    int mpiSize = mpiComm->Get_size();
221
222
223
224
225

    TraverseStack stack;
    ElInfo *elInfo;

    // === create parmetis mesh ===
226
227
228
    if (parMetisMesh) 
      DELETE parMetisMesh;
    parMetisMesh = NEW ParMetisMesh(mesh_, mpiComm);
229

230
    int nElements = parMetisMesh->getNumElements();
231
232

    // === create weight array ===
233
234
    int *wgts = elemWeights ? GET_MEMORY(int, nElements) : NULL;
    float *floatWgts = elemWeights ? GET_MEMORY(float, nElements) : NULL;
235
236
237
238
    float maxWgt = 0.0;
    float *ptr_floatWgts = floatWgts;

    elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER);
239
    while (elInfo) {
240
241
242
243
244
245
      Element *element = elInfo->getElement();

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

246
247
248
      if (partitionData &&
	  partitionData->getPartitionStatus() == IN &&
	  partitionData->getLevel() == 0) {
249
250
251
252
253
254
255
256
257
258
259
260
261
262
	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;
263
    mpiComm->Allreduce(&maxWgt, &tmp, 1, MPI_FLOAT, MPI_MAX);
264
265
266
    maxWgt = tmp;

    // === create dual graph ===
267
    ParMetisGraph parMetisGraph(parMetisMesh, mpiComm);
268
269
270
271
272
273
274
275
276
277

    // === 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
    float *tpwgts = elemWeights ? GET_MEMORY(float, mpiSize) : NULL;
    float ubvec = 1.05;
    int options[3] = {0, 0, 0}; // default options
    int edgecut = -1;
278
    int *part = GET_MEMORY(int, nElements);
279
    
280
281
    if (elemWeights) {
      for (int i = 0; i < mpiSize; i++) {
282
	// set tpwgts
283
	tpwgts[i] = 1.0 / nparts;
284
285
286
      }      

      float scale = 10000 / maxWgt;
287
      for (int i = 0; i < nElements; i++) {
288
289
290
291
292
	// scale wgts
	wgts[i] = static_cast<int>(floatWgts[i] * scale);
      }
    }

293
294
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

295
296
    switch(mode) {
    case INITIAL:
297
      ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
298
299
300
301
302
303
304
305
306
307
308
309
310
			   parMetisGraph.getXAdj(),
			   parMetisGraph.getAdjncy(),
			   wgts,
			   NULL,
			   &wgtflag,
			   &numflag,
			   &ncon,
			   &nparts,
			   tpwgts,
			   &ubvec,
			   options,
			   &edgecut,
			   part,
311
			   &tmpComm);
312
313
314
      break;
    case ADAPTIVE_REPART:
      {
315
316
	int *vsize = GET_MEMORY(int, nElements);
	for (int i = 0; i < nElements; i++) {
317
318
	  vsize[i] = 1;
	}
319
	ParMETIS_V3_AdaptiveRepart(parMetisMesh->getElementDist(),
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
				   parMetisGraph.getXAdj(),
				   parMetisGraph.getAdjncy(),
				   wgts,
				   NULL,
				   vsize,
				   &wgtflag,
				   &numflag,
				   &ncon,
				   &nparts,
				   tpwgts,
				   &ubvec,
				   &itr,
				   options,
				   &edgecut,
				   part,
335
				   &tmpComm);
336
	FREE_MEMORY(vsize, int, nElements);
337
338
339
      }
      break;
    case REFINE_PART:
340
      ParMETIS_V3_RefineKway(parMetisMesh->getElementDist(),
341
342
343
344
345
346
347
348
349
350
351
352
353
			     parMetisGraph.getXAdj(),
			     parMetisGraph.getAdjncy(),
			     wgts,
			     NULL,
			     &wgtflag,
			     &numflag,
			     &ncon,
			     &nparts,
			     tpwgts,
			     &ubvec,
			     options,
			     &edgecut,
			     part,
354
			     &tmpComm);
355
356
357
358
359
360
361
362
      break;
    default: 
      ERROR_EXIT("unknown partitioning mode\n");
    }

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

363
    if (floatWgts) 
364
      FREE_MEMORY(floatWgts, float, nElements);
365
    if (wgts) 
366
      FREE_MEMORY(wgts, int, nElements);    
367
368
    if (tpwgts) 
      FREE_MEMORY(tpwgts, float, mpiSize);
369
    FREE_MEMORY(part, int, nElements);
370
371
372
373
374
375
376
377
378
  }

  void ParMetisPartitioner::fillCoarsePartitionVec(std::map<int, int> *partitionVec)
  {
    TEST_EXIT(partitionVec)("no partition vector\n");

    partitionVec->clear();

    // update ParMETIS mesh to new partitioning
379
380
    if (!parMetisMesh) 
      parMetisMesh = NEW ParMetisMesh(mesh_, mpiComm);
381

382
383
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
384
385
    int *nPartitionElements = GET_MEMORY(int, mpiSize);
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
386

387
    for (int i = 0;  i < mpiSize; i++) {
388
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
389
390
391
    }

    // === count number of elements ===
392
393
394
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
395

396
    int *partitionElements = GET_MEMORY(int, nElements);
397
398

    // distribute partition elements
399
400
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
401
402
			MPI_INT, 
			partitionElements, 
403
			nPartitionElements, 
404
405
			elmdist, 
			MPI_INT);
406
407

    // fill partitionVec
408
    for (int i = 0; i < mpiSize; i++) {
409
      for (int j = 0; j < nPartitionElements[i]; j++) {
410
411
412
413
	(*partitionVec)[partitionElements[elmdist[i] + j]] = i;
      }
    }

414
415
    FREE_MEMORY(partitionElements, int, nElements);
    FREE_MEMORY(nPartitionElements, int, mpiSize);
416
417
418
419
  }

  void ParMetisPartitioner::distributePartitioning(int *part) 
  {
420
421
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
422
    int nElements = parMetisMesh->getNumElements();
423

424
425
    // nPartitionElements[i] is the number of elements for the i-th partition
    int *nPartitionElements = GET_MEMORY(int, mpiSize);
426
    for (int i = 0; i < mpiSize; i++) 
427
428
429
      nPartitionElements[i] = 0;
    for (int i = 0; i < nElements; i++) {
      nPartitionElements[part[i]]++;
430
431
432
    }

    // collect number of partition elements from all ranks for this rank
433
434
435
    int *nRankElements = GET_MEMORY(int, mpiSize);
    mpiComm->Alltoall(nPartitionElements, 1, MPI_INT,
		      nRankElements, 1, MPI_INT);
436
437
438

    // sum up partition elements over all ranks
    int *sumPartitionElements = GET_MEMORY(int, mpiSize);
439
    mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize,
440
		       MPI_INT, MPI_SUM);
441
442
443
444
445

    
    // prepare distribution (fill partitionElements with AMDiS indices)
    int *bufferOffset = GET_MEMORY(int, mpiSize);
    bufferOffset[0] = 0;
446
    for (int i = 1; i < mpiSize; i++) {
447
      bufferOffset[i] = bufferOffset[i - 1] + nPartitionElements[i - 1];
448
449
    }

450
    int *partitionElements = GET_MEMORY(int, nElements);
451
452
    int **partitionPtr = GET_MEMORY(int*, mpiSize);

453
    for (int i = 0; i < mpiSize; i++) {
454
455
456
      partitionPtr[i] = partitionElements + bufferOffset[i];
    }

457
    for (int i = 0; i < nElements; i++) {
458
      int partition = part[i];
459
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
460
461
462
463
464
465
466
467
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

    // all to all: partition elements to rank elements
    int *rankElements = GET_MEMORY(int, sumPartitionElements[mpiRank]);
    int *recvBufferOffset = GET_MEMORY(int, mpiSize);
    recvBufferOffset[0] = 0;
468
    for (int i = 1; i < mpiSize; i++) {
469
      recvBufferOffset[i] = recvBufferOffset[i - 1] + nRankElements[i - 1];
470
471
    }

472
    mpiComm->Alltoallv(partitionElements, 
473
		       nPartitionElements,
474
475
476
		       bufferOffset,
		       MPI_INT,
		       rankElements,
477
		       nRankElements,
478
479
		       recvBufferOffset,
		       MPI_INT);
480
481
    

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

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER);
495
    while (elInfo) {
496
497
498
499
500
501
      Element *element = elInfo->getElement();

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

502
      if (partitionData && partitionData->getLevel() == 0) {
503
	int amdisIndex = element->getIndex();
504
	if (elementInPartition[amdisIndex]) {
505
506
507
508
509
510
511
512
513
514
	  partitionData->setPartitionStatus(IN);
	} else {
	  partitionData->setPartitionStatus(OUT);
	}
	descendPartitionData(element);
      }

      elInfo = stack.traverseNext(elInfo);
    }

515
516
    DELETE parMetisMesh;
    parMetisMesh = NULL;
517
518

    FREE_MEMORY(rankElements, int, sumPartitionElements[mpiRank]);
519
520
    FREE_MEMORY(nPartitionElements, int, mpiSize);
    FREE_MEMORY(nRankElements, int, mpiSize);
521
    FREE_MEMORY(sumPartitionElements, int, mpiSize);
522
    FREE_MEMORY(partitionElements, int, nElements);
523
524
525
526
527
528
529
    FREE_MEMORY(partitionPtr, int*, mpiSize);
    FREE_MEMORY(bufferOffset, int, mpiSize);
    FREE_MEMORY(recvBufferOffset, int, mpiSize);
  }

  void ParMetisPartitioner::descendPartitionData(Element *element) 
  {
530
    if (!element->isLeaf()) {
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
      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;
    ElInfo *elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER);
559
    while (elInfo) {
560
561
562
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));
563
564
      if (partitionData) {
	if (partitionData->getLevel() == 0) {
565
566
567
568
569
570
571
572
573
574
	  partition = (*(coarseVec))[element->getIndex()];
	}
	if(element->isLeaf()) {
	  (*(fineVec))[element->getIndex()] = partition;
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }
  }
}