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
425
426
427
428
429
    bool b = distributePartitioning(&(part[0]));

    if (b)
      createPartitionMap();

    return b;
430
431
  }

432

433
  void ParMetisPartitioner::createPartitionMap()
434
  {
435
    FUNCNAME("ParMetisPartitioner::createPartitionMap()");
436

437
    partitionMap.clear();
438
439

    // update ParMETIS mesh to new partitioning
440
    if (!parMetisMesh)
441
      parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
442

443
444
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
445
    vector<int> nPartitionElements(mpiSize);
446
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
447

448
    for (int i = 0; i < mpiSize; i++)
449
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
450
451

    // === count number of elements ===
452
453
454
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
455

456
    vector<int> partitionElements(nElements);
457
458

    // distribute partition elements
459
460
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
461
			MPI_INT, 
462
463
			&(partitionElements[0]), 
			&(nPartitionElements[0]), 
464
465
			elmdist, 
			MPI_INT);
466

467
    // fill partitionMap
Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
    for (int i = 0; i < mpiSize; i++)
      for (int j = 0; j < nPartitionElements[i]; j++)
470
	partitionMap[partitionElements[elmdist[i] + j]] = i;
471
472
  }

473

474
  bool ParMetisPartitioner::distributePartitioning(int *part) 
475
  {
476
477
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

478
479
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
480
    int nElements = parMetisMesh->getNumElements();
481

482
    // nPartitionElements[i] is the number of elements for the i-th partition
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    int *nPartitionElements = new int[mpiSize];
484
    for (int i = 0; i < mpiSize; i++) 
485
      nPartitionElements[i] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
486
    for (int i = 0; i < nElements; i++)
487
      nPartitionElements[part[i]]++;    
488
489

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

493

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
513
514
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
515

516
    for (int i = 0; i < mpiSize; i++)
517
518
      partitionPtr[i] = partitionElements + bufferOffset[i];

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

524
525
526
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

527
528
529
530
531
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
565
    delete parMetisMesh;
566
    parMetisMesh = NULL;
567

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

    return true;
578
  }
579
 
580
}