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 <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
302
303
304


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

305
306
307
308
309
    int smin = 9999999;
    int smax = 0;
    int ssum = 0;

    for (int i = 0; i < nElements; i++) {
310
      wgts[i] = static_cast<int>(floatWgts[i] * scale);
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
      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)
326
	kpartMax = std::max(kpartMax, wgts[i]);
327
328
329
330
331
332
333
334
335
336

    mpi::globalMax(kpartMax);

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

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

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

340
341
342
343
344
345
346
347
348
349
      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);
350
351
352


    // === 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
	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
    return distributePartitioning(&(part[0]));
425
426
  }

427

428
  void ParMetisPartitioner::createPartitionMap(map<int, int>& pMap)
429
  {
430
    FUNCNAME("ParMetisPartitioner::createPartitionMap()");
431

432
    partitionMap.clear();
433
434

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

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

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

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

451
    vector<int> partitionElements(nElements);
452
453

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

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

    pMap = partitionMap;
468
469
  }

470

471
  bool ParMetisPartitioner::distributePartitioning(int *part) 
472
  {
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
497
498
499
500
501
502
503

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
510
511
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
512

513
    for (int i = 0; i < mpiSize; i++)
514
515
      partitionPtr[i] = partitionElements + bufferOffset[i];

516
    sendElements.clear();
517
    for (int i = 0; i < nElements; i++) {
518
      int partition = part[i];
519
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
520

521
522
523
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

524
525
526
527
528
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

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

535
    mpiComm->Alltoallv(partitionElements, 
536
		       nPartitionElements,
537
538
539
		       bufferOffset,
		       MPI_INT,
		       rankElements,
540
		       nRankElements,
541
542
		       recvBufferOffset,
		       MPI_INT);
543
    
544
    TEST_EXIT(elementInRank.size() != 0)("Should not happen!\n");
545
    for (map<int, bool>::iterator it = elementInRank.begin();
546
547
	 it != elementInRank.end(); ++it)
      elementInRank[it->first] = false;
548

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

Thomas Witkowski's avatar
Thomas Witkowski committed
562
    delete parMetisMesh;
563
    parMetisMesh = NULL;
564

Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
567
568
569
570
571
    delete [] rankElements;
    delete [] nPartitionElements;
    delete [] nRankElements;
    delete [] sumPartitionElements;
    delete [] partitionElements;
    delete [] partitionPtr;
    delete [] bufferOffset;
Thomas Witkowski's avatar
huh    
Thomas Witkowski committed
572
    delete [] recvBufferOffset;
573
574

    return true;
575
  }
576
 
577
}