ParMetisPartitioner.cc 14 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::getPartitionMap(map<int, int> &partitionMap)
429
  {
430
    FUNCNAME("ParMetisPartitioner::getPartitionMap()");
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
  }

468

469
  bool ParMetisPartitioner::distributePartitioning(int *part) 
470
  {
471
472
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

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

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

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

488

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

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

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

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

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

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

519
520
521
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

522
523
524
525
526
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
560
    delete parMetisMesh;
561
    parMetisMesh = NULL;
562

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

    return true;
573
  }
574
 
575
}