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

namespace AMDiS {

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

48
    nElements = elementCounter;
49

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

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

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

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

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

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

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

    elementCounter = 0;
81
    int nodeCounter = 0;
82

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

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

	// write eptr entry
95
	nodeCounter += mesh->getGeo(VERTEX);
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
	  if (mapLocalGlobal)
102
	    *ptr_eind = (*mapLocalGlobal)[element->getDof(i, 0)];
103
104
	  else
	    *ptr_eind = element->getDof(i, 0);	  
105

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

	elementCounter++;
      }
120

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

125

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

144

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

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

155
156
    int numflag = 0;

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

160
    MPI_Comm tmpComm = MPI_Comm(*comm);
161

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

172

173
174
  ParMetisGraph::~ParMetisGraph()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
175
176
    free(xadj);
    free(adjncy);
177
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
  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;
  }


224
225
  void ParMetisPartitioner::createPartitionData() 
  {
226
227
    FUNCNAME("ParMetrisPartitioner::createPartitionData()");

228
229
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
230
231
    int nLeaves = mesh->getNumberOfLeaves();
    int elPerRank = nLeaves / mpiSize;
232

233
    // === Create initial partitioning of the AMDiS mesh. ===
234

235
236
    elementInRank.clear();

237
    TraverseStack stack;
238
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
239
    while (elInfo) {
240
241
      Element *element = elInfo->getElement();

242
243
244
245
      if ((element->getIndex() >= mpiRank * elPerRank &&
	   element->getIndex() < (mpiRank + 1) * elPerRank) ||
	  (element->getIndex() >= mpiSize * elPerRank &&
	   mpiRank == mpiSize - 1))
246
	elementInRank[element->getIndex()] = true;
247
      else
248
	elementInRank[element->getIndex()] = false;
249
      
250
251
252
253
      elInfo = stack.traverseNext(elInfo);
    }
  }

254

255
  bool ParMetisPartitioner::partition(std::map<int, double> &elemWeights,
256
257
258
				      PartitionMode mode,
				      float itr) 
  {
259
    FUNCNAME("ParMetisPartitioner::partition()");
260

261
    int mpiSize = mpiComm->Get_size();
262

263
264
265

    // === Create parmetis mesh ===

266
    if (parMetisMesh) 
Thomas Witkowski's avatar
Thomas Witkowski committed
267
      delete parMetisMesh;
268

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

271
    parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
272

273
    int nElements = parMetisMesh->getNumElements();
274

275
276
277

    // === Create weight array ===

278
279
280
    std::vector<int> wgts(nElements);
    std::vector<float> floatWgts(nElements);
    unsigned int floatWgtsPos = 0;
281
282
    float maxWgt = 0.0;

283
    TraverseStack stack;
284
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
285
    while (elInfo) {
286
      int index = elInfo->getElement()->getIndex();
287

288
      if (elementInRank[index]) {
289
	// get weight 
290
	float wgt = static_cast<float>(elemWeights[index]);
291
292
293
	maxWgt = max(wgt, maxWgt);

	// write float weight
294
295
	TEST_EXIT_DBG(floatWgtsPos < floatWgts.size())("Should not happen!\n");
	floatWgts[floatWgtsPos++] = wgt;
296
297
298
299
      }
      elInfo = stack.traverseNext(elInfo);
    }

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

302
    float tmp;
303
    mpiComm->Allreduce(&maxWgt, &tmp, 1, MPI_FLOAT, MPI_MAX);
304
305
    maxWgt = tmp;

306
307
308

    // === Create dual graph ===

309
    ParMetisGraph parMetisGraph(parMetisMesh, mpiComm);
310

311
312
313

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

314
    int wgtflag = 2; // weights at vertices only!
315
    int numflag = 0; // c numbering style!
316
    int ncon = 1; // one weight at each vertex!
317
    int nparts = mpiSize; // number of partitions
318

319
    std::vector<float> tpwgts(mpiSize);
320
    float ubvec = 1.05;
321
    int options[4] = {0, 0, 15, 1}; // default options
322
    int edgecut = -1;
323
    std::vector<int> part(nElements);
324

325
326
327
    // set tpwgts
    for (int i = 0; i < mpiSize; i++)
      tpwgts[i] = 1.0 / nparts;
328
   
329
    float scale = 10000.0 / maxWgt;
330
331
332
333


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

334
    for (int i = 0; i < nElements; i++)
335
      wgts[i] = static_cast<int>(floatWgts[i] * scale);
336
337
338


    // === Start ParMETIS. ===
339

340
341
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

342
    switch (mode) {
343
    case INITIAL:
344
      ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
345
346
			   parMetisGraph.getXAdj(),
			   parMetisGraph.getAdjncy(),
347
			   &(wgts[0]),
348
			   NULL,
349
			   &wgtflag,
350
351
352
			   &numflag,
			   &ncon,
			   &nparts,
353
			   &(tpwgts[0]),
354
355
356
			   &ubvec,
			   options,
			   &edgecut,
357
			   &(part[0]),
358
			   &tmpComm);
359
360
361
      break;
    case ADAPTIVE_REPART:
      {
362
	std::vector<int> vsize(nElements);
Thomas Witkowski's avatar
Thomas Witkowski committed
363
	for (int i = 0; i < nElements; i++)
364
365
	  vsize[i] = static_cast<int>(floatWgts[i]);

366
	ParMETIS_V3_AdaptiveRepart(parMetisMesh->getElementDist(),
367
368
				   parMetisGraph.getXAdj(),
				   parMetisGraph.getAdjncy(),
369
				   &(wgts[0]),
370
				   NULL,
371
				   &(vsize[0]),
372
373
374
375
				   &wgtflag,
				   &numflag,
				   &ncon,
				   &nparts,
376
				   &(tpwgts[0]),
377
378
379
380
				   &ubvec,
				   &itr,
				   options,
				   &edgecut,
381
				   &(part[0]),
382
				   &tmpComm);
383
384
385
      }
      break;
    case REFINE_PART:
386
      ParMETIS_V3_RefineKway(parMetisMesh->getElementDist(),
387
388
			     parMetisGraph.getXAdj(),
			     parMetisGraph.getAdjncy(),
389
			     &(wgts[0]),
390
391
392
393
394
			     NULL,
			     &wgtflag,
			     &numflag,
			     &ncon,
			     &nparts,
395
			     &(tpwgts[0]),
396
397
398
			     &ubvec,
			     options,
			     &edgecut,
399
			     &(part[0]),
400
			     &tmpComm);
401

402
403
404
405
406
      break;
    default: 
      ERROR_EXIT("unknown partitioning mode\n");
    }

407
408
409

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

410
    return distributePartitioning(&(part[0]));
411
412
  }

413

414
415
  void ParMetisPartitioner::fillCoarsePartitionVec(std::map<int, int> *partitionVec)
  {
416
417
    FUNCNAME("ParMetisPartitioner::fillCoarsePartitionVec()");

418
    TEST_EXIT_DBG(partitionVec)("No partition vector!\n");
419
420
421
422

    partitionVec->clear();

    // update ParMETIS mesh to new partitioning
423
    if (!parMetisMesh)
424
      parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
425

426
427
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
428
    std::vector<int> nPartitionElements(mpiSize);
429
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
430

431
    for (int i = 0; i < mpiSize; i++)
432
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
433
434

    // === count number of elements ===
435
436
437
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
438

439
    std::vector<int> partitionElements(nElements);
440
441

    // distribute partition elements
442
443
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
444
			MPI_INT, 
445
446
			&(partitionElements[0]), 
			&(nPartitionElements[0]), 
447
448
			elmdist, 
			MPI_INT);
449
450

    // fill partitionVec
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
    for (int i = 0; i < mpiSize; i++)
      for (int j = 0; j < nPartitionElements[i]; j++)
453
454
455
	(*partitionVec)[partitionElements[elmdist[i] + j]] = i;
  }

456

457
  bool ParMetisPartitioner::distributePartitioning(int *part) 
458
  {
459
460
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

461
462
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
463
    int nElements = parMetisMesh->getNumElements();
464

465
    // nPartitionElements[i] is the number of elements for the i-th partition
Thomas Witkowski's avatar
Thomas Witkowski committed
466
    int *nPartitionElements = new int[mpiSize];
467
    for (int i = 0; i < mpiSize; i++) 
468
      nPartitionElements[i] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
469
    for (int i = 0; i < nElements; i++)
470
      nPartitionElements[part[i]]++;    
471
472

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

476

477
    // sum up partition elements over all ranks
Thomas Witkowski's avatar
Thomas Witkowski committed
478
    int *sumPartitionElements = new int[mpiSize];
479
    mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize,
480
		       MPI_INT, MPI_SUM);
481
482
483
484
485
486
487
488
489

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

490
    // prepare distribution (fill partitionElements with AMDiS indices)
Thomas Witkowski's avatar
Thomas Witkowski committed
491
    int *bufferOffset = new int[mpiSize];
492
    bufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
493
    for (int i = 1; i < mpiSize; i++)
494
      bufferOffset[i] = bufferOffset[i - 1] + nPartitionElements[i - 1];
495

Thomas Witkowski's avatar
Thomas Witkowski committed
496
497
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
498

499
    for (int i = 0; i < mpiSize; i++)
500
501
      partitionPtr[i] = partitionElements + bufferOffset[i];

502
    sendElements.clear();
503
    for (int i = 0; i < nElements; i++) {
504
      int partition = part[i];
505
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
506

507
508
509
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

510
511
512
513
514
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

    // all to all: partition elements to rank elements
Thomas Witkowski's avatar
Thomas Witkowski committed
515
516
    int *rankElements = new int[sumPartitionElements[mpiRank]];
    int *recvBufferOffset = new int[mpiSize];
517
    recvBufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
518
    for (int i = 1; i < mpiSize; i++)
519
      recvBufferOffset[i] = recvBufferOffset[i - 1] + nRankElements[i - 1];
520

521
    mpiComm->Alltoallv(partitionElements, 
522
		       nPartitionElements,
523
524
525
		       bufferOffset,
		       MPI_INT,
		       rankElements,
526
		       nRankElements,
527
528
		       recvBufferOffset,
		       MPI_INT);
529
    
530
531
532
533
    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;
534

535
536
    // Create map which stores for each element index on ther partitioning level
    // if the element is in the partition of this rank.
537
    recvElements.clear();
538
    for (int i = 0; i < mpiSize; i++) {
539
      int *rankStart = rankElements + recvBufferOffset[i];
540
      int *rankEnd = rankStart + nRankElements[i];
541
542
      for (int *rankPtr = rankStart; rankPtr < rankEnd; ++rankPtr) {
	elementInRank[*rankPtr] = true;
543
544
	if (i != mpiRank)
	  recvElements[i].push_back(*rankPtr);
545
546
547
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
548
    delete parMetisMesh;
549
    parMetisMesh = NULL;
550

Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
553
554
555
556
557
    delete [] rankElements;
    delete [] nPartitionElements;
    delete [] nRankElements;
    delete [] sumPartitionElements;
    delete [] partitionElements;
    delete [] partitionPtr;
    delete [] bufferOffset;
Thomas Witkowski's avatar
huh    
Thomas Witkowski committed
558
    delete [] recvBufferOffset;
559
560

    return true;
561
562
563
  }

}