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
#include "Serializer.h"
17
18
19
20
21
22
23
24
25
26
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "mpi.h"

namespace AMDiS {

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

49
    nElements = elementCounter;
50

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

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

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

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

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

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

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

    elementCounter = 0;
82
    int nodeCounter = 0;
83

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

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

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

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

107
108
109
110
	  ptr_eind++;
	}

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

	elementCounter++;
      }
121

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

126

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

145

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

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

156
157
    int numflag = 0;

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

161
    MPI_Comm tmpComm = MPI_Comm(*comm);
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
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
224
  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;
  }


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

230
    int mpiSize = mpiComm->Get_size();
231

232
233
234

    // === Create parmetis mesh ===

235
    if (parMetisMesh) 
Thomas Witkowski's avatar
Thomas Witkowski committed
236
      delete parMetisMesh;
237

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

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

242
    int nElements = parMetisMesh->getNumElements();
243

244
245
246

    // === Create weight array ===

247
248
249
    std::vector<int> wgts(nElements);
    std::vector<float> floatWgts(nElements);
    unsigned int floatWgtsPos = 0;
250
251
    float maxWgt = 0.0;

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

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

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

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

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

275
276
277

    // === Create dual graph ===

278
    ParMetisGraph parMetisGraph(parMetisMesh, mpiComm);
279

280
281
282

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

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

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

294
295
296
    // set tpwgts
    for (int i = 0; i < mpiSize; i++)
      tpwgts[i] = 1.0 / nparts;
297
   
298
    float scale = 10000.0 / maxWgt;
299
300
301
302


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

303
304
305
306
307
    int smin = 9999999;
    int smax = 0;
    int ssum = 0;

    for (int i = 0; i < nElements; i++) {
308
      wgts[i] = static_cast<int>(floatWgts[i] * scale);
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
      smin = std::min(smin, wgts[i]);
      smax = std::max(smax, wgts[i]);
      ssum += wgts[i];
    }

    mpi::globalMin(smin);
    mpi::globalMax(smax);
    mpi::globalAdd(ssum);

    MSG("DATA   SMIN = %d   SMAX = %d    SSUM = %d\n", smin, smax, ssum);

    int kpart = ssum / mpiSize;
    int kpartMax = 0;
    for (int i = 0; i < nElements; i++)
      if (wgts[i] < kpart)
324
	kpartMax = std::max(kpartMax, wgts[i]);
325
326
327
328
329
330
331
332
333
334

    mpi::globalMax(kpartMax);

    MSG("KPART MAX: %d\n", kpartMax);

    smin = 9999999;
    smax = 0;
    ssum = 0;

    for (int i = 0; i < nElements; i++) {
335
336
      wgts[i] = std::min(wgts[i], kpartMax);
      wgts[i] = std::max(wgts[i], 1);
337

338
339
340
341
342
343
344
345
346
347
      smin = std::min(smin, wgts[i]);
      smax = std::max(smax, wgts[i]);
      ssum += wgts[i];
    }

    mpi::globalMin(smin);
    mpi::globalMax(smax);
    mpi::globalAdd(ssum);

    MSG("DATA   SMIN = %d   SMAX = %d    SSUM = %d\n", smin, smax, ssum);
348
349
350


    // === Start ParMETIS. ===
351

352
353
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

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

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

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

419
420
421

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

422
    return distributePartitioning(&(part[0]));
423
424
  }

425

426
  void ParMetisPartitioner::getPartitionMap(std::map<int, int> &partitionMap)
427
  {
428
    FUNCNAME("ParMetisPartitioner::getPartitionMap()");
429

430
    partitionMap.clear();
431
432

    // update ParMETIS mesh to new partitioning
433
    if (!parMetisMesh)
434
      parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
435

436
437
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
438
    std::vector<int> nPartitionElements(mpiSize);
439
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
440

441
    for (int i = 0; i < mpiSize; i++)
442
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
443
444

    // === count number of elements ===
445
446
447
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
448

449
    std::vector<int> partitionElements(nElements);
450
451

    // distribute partition elements
452
453
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
454
			MPI_INT, 
455
456
			&(partitionElements[0]), 
			&(nPartitionElements[0]), 
457
458
			elmdist, 
			MPI_INT);
459
460

    // fill partitionVec
Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
    for (int i = 0; i < mpiSize; i++)
      for (int j = 0; j < nPartitionElements[i]; j++)
463
	partitionMap[partitionElements[elmdist[i] + j]] = i;
464
465
  }

466

467
  bool ParMetisPartitioner::distributePartitioning(int *part) 
468
  {
469
470
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

471
472
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
473
    int nElements = parMetisMesh->getNumElements();
474

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

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

486

487
    // sum up partition elements over all ranks
Thomas Witkowski's avatar
Thomas Witkowski committed
488
    int *sumPartitionElements = new int[mpiSize];
489
    mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize,
490
		       MPI_INT, MPI_SUM);
491
492
493
494
495
496
497
498
499

    // 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;    

500
    // prepare distribution (fill partitionElements with AMDiS indices)
Thomas Witkowski's avatar
Thomas Witkowski committed
501
    int *bufferOffset = new int[mpiSize];
502
    bufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
503
    for (int i = 1; i < mpiSize; i++)
504
      bufferOffset[i] = bufferOffset[i - 1] + nPartitionElements[i - 1];
505

Thomas Witkowski's avatar
Thomas Witkowski committed
506
507
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
508

509
    for (int i = 0; i < mpiSize; i++)
510
511
      partitionPtr[i] = partitionElements + bufferOffset[i];

512
    sendElements.clear();
513
    for (int i = 0; i < nElements; i++) {
514
      int partition = part[i];
515
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
516

517
518
519
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

520
521
522
523
524
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

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

531
    mpiComm->Alltoallv(partitionElements, 
532
		       nPartitionElements,
533
534
535
		       bufferOffset,
		       MPI_INT,
		       rankElements,
536
		       nRankElements,
537
538
		       recvBufferOffset,
		       MPI_INT);
539
    
540
541
542
543
    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;
544

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

Thomas Witkowski's avatar
Thomas Witkowski committed
558
    delete parMetisMesh;
559
    parMetisMesh = NULL;
560

Thomas Witkowski's avatar
Thomas Witkowski committed
561
562
563
564
565
566
567
    delete [] rankElements;
    delete [] nPartitionElements;
    delete [] nRankElements;
    delete [] sumPartitionElements;
    delete [] partitionElements;
    delete [] partitionPtr;
    delete [] bufferOffset;
Thomas Witkowski's avatar
huh    
Thomas Witkowski committed
568
    delete [] recvBufferOffset;
569
570

    return true;
571
  }
572
 
573
}