ParMetisPartitioner.cc 13.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
#include "parallel/ParMetisPartitioner.h"
17
#include "parallel/ParallelDofMapping.h"
18
#include "parallel/MpiHelper.h"
19
#include "Serializer.h"
20
21
22
23
24
25
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "FixVec.h"
#include "DOFVector.h"
26

27
28
29

namespace AMDiS {

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

    int mpiSize = mpiComm->Get_size();
40
41
42
43
    int elementCounter = 0;
    int dow = Global::getGeo(WORLD);

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

      elInfo = stack.traverseNext(elInfo);
    }

52
    nElements = elementCounter;
53

54
    TEST_EXIT(nElements > 0)("No elements in ParMETIS mesh!\n");
55
56

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
67
    eptr[0] = 0;
68

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

76
    elmdist[0] = 0;
77
    for (int i = 2; i < mpiSize + 1; i++)
78
      elmdist[i] += elmdist[i - 1];
79
80

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

    elementCounter = 0;
85
    int nodeCounter = 0;
86

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

      // if element in partition
93
      if (elementInRank[index]) {
94
95
96
97
98
	// remember index
	setParMetisIndex(index, elementCounter);
	setAMDiSIndex(elementCounter, index);

	// write eptr entry
99
	nodeCounter += mesh->getGeo(VERTEX);
100
101
102
103
	*ptr_eptr = nodeCounter;
	ptr_eptr++;

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

110
111
112
113
	  ptr_eind++;
	}

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

	elementCounter++;
      }
124

125
126
127
128
      elInfo = stack.traverseNext(elInfo);
    }
  }

129

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

148

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

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

159
160
    int numflag = 0;

161
162
    if (ncommonnodes == -1) 
      ncommonnodes = parMetisMesh->getDim();
163

164
    MPI_Comm tmpComm = MPI_Comm(*comm);
165

166
167
168
    ParMETIS_V3_Mesh2Dual(parMetisMesh->getElementDist(),
			  parMetisMesh->getElementPtr(),
			  parMetisMesh->getElementInd(),
169
170
			  &numflag,
			  &ncommonnodes,
Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
			  &xadj,
			  &adjncy,
173
			  &tmpComm);
174
175
  }

176

177
178
  ParMetisGraph::~ParMetisGraph()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
    free(xadj);
    free(adjncy);
181
182
  }

183

184
185
186
187
  void ParMetisGraph::print()
  {
    FUNCNAME("ParMetisGraph::print()");

188
    stringstream oss;
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
227
    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;
  }


228
  bool ParMetisPartitioner::partition(map<int, double> &elemWeights,
229
				      PartitionMode mode) 
230
  {
231
    FUNCNAME("ParMetisPartitioner::partition()");
232

233
    int mpiSize = mpiComm->Get_size();
234

235
236
237

    // === Create parmetis mesh ===

238
    if (parMetisMesh) 
Thomas Witkowski's avatar
Thomas Witkowski committed
239
      delete parMetisMesh;
240

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

243
    parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
244

245
    int nElements = parMetisMesh->getNumElements();
246

247
248
249

    // === Create weight array ===

250
251
    vector<int> wgts(nElements);
    vector<float> floatWgts(nElements);
252
    unsigned int floatWgtsPos = 0;
253
254
    float maxWgt = 0.0;

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

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

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

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

274
    float tmp;
275
    mpiComm->Allreduce(&maxWgt, &tmp, 1, MPI_FLOAT, MPI_MAX);
276
277
    maxWgt = tmp;

278
279
280

    // === Create dual graph ===

281
    ParMetisGraph parMetisGraph(parMetisMesh, mpiComm);
282

283
284
285

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

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

291
292
    vector<double> tpwgts(mpiSize);
    double ubvec = 1.05;
293
    int options[4] = {0, 0, 15, 1}; // default options
294
    int edgecut = -1;
295
    vector<int> part(nElements);
296

297
298
299
    // set tpwgts
    for (int i = 0; i < mpiSize; i++)
      tpwgts[i] = 1.0 / nparts;
300
   
301
    float scale = 10000.0 / maxWgt;
302
    for (int i = 0; i < nElements; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
303
      wgts[i] = static_cast<int>(floatWgts[i] * scale);
304
305
306


    // === Start ParMETIS. ===
307

308
309
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

310
    switch (mode) {
311
    case INITIAL:
312
      ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
313
314
			   parMetisGraph.getXAdj(),
			   parMetisGraph.getAdjncy(),
315
			   &(wgts[0]),
316
			   NULL,
317
			   &wgtflag,
318
319
320
			   &numflag,
			   &ncon,
			   &nparts,
321
			   &(tpwgts[0]),
322
323
324
			   &ubvec,
			   options,
			   &edgecut,
325
			   &(part[0]),
326
			   &tmpComm);
327
328
329
      break;
    case ADAPTIVE_REPART:
      {
330
	vector<int> vsize(nElements);
Thomas Witkowski's avatar
Thomas Witkowski committed
331
	for (int i = 0; i < nElements; i++)
332
333
	  vsize[i] = static_cast<int>(floatWgts[i]);

334
	ParMETIS_V3_AdaptiveRepart(parMetisMesh->getElementDist(),
335
336
				   parMetisGraph.getXAdj(),
				   parMetisGraph.getAdjncy(),
337
				   &(wgts[0]),
338
				   NULL,
339
				   &(vsize[0]),
340
341
342
343
				   &wgtflag,
				   &numflag,
				   &ncon,
				   &nparts,
344
				   &(tpwgts[0]),
345
346
347
348
				   &ubvec,
				   &itr,
				   options,
				   &edgecut,
349
				   &(part[0]),
350
				   &tmpComm);
351
352
353
      }
      break;
    case REFINE_PART:
354
      ParMETIS_V3_RefineKway(parMetisMesh->getElementDist(),
355
356
			     parMetisGraph.getXAdj(),
			     parMetisGraph.getAdjncy(),
357
			     &(wgts[0]),
358
359
360
361
362
			     NULL,
			     &wgtflag,
			     &numflag,
			     &ncon,
			     &nparts,
363
			     &(tpwgts[0]),
364
365
366
			     &ubvec,
			     options,
			     &edgecut,
367
			     &(part[0]),
368
			     &tmpComm);
369

370
371
372
373
374
      break;
    default: 
      ERROR_EXIT("unknown partitioning mode\n");
    }

375
376
377

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

378
    return distributePartitioning(&(part[0]));
379
380
  }

381

382
  void ParMetisPartitioner::createPartitionMap(map<int, int>& pMap)
383
  {
384
    FUNCNAME("ParMetisPartitioner::createPartitionMap()");
385

386
    partitionMap.clear();
387
388

    // update ParMETIS mesh to new partitioning
389
    if (!parMetisMesh)
390
      parMetisMesh = new ParMetisMesh(mesh, mpiComm, elementInRank, mapLocalGlobal);
391

392
393
    int mpiRank = mpiComm->Get_rank();
    int mpiSize = mpiComm->Get_size();
394
    vector<int> nPartitionElements(mpiSize);
395
    int *elmdist = parMetisMesh->getElementDist();
Thomas Witkowski's avatar
Thomas Witkowski committed
396

397
    for (int i = 0; i < mpiSize; i++)
398
      nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
399
400

    // === count number of elements ===
401
402
403
    int nElements = 0;
    int localElements = parMetisMesh->getNumElements();
    mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
404

405
    vector<int> partitionElements(nElements);
406
407

    // distribute partition elements
408
409
    mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
			nPartitionElements[mpiRank], 
410
			MPI_INT, 
411
412
			&(partitionElements[0]), 
			&(nPartitionElements[0]), 
413
414
			elmdist, 
			MPI_INT);
415

416
    // fill partitionMap
Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
    for (int i = 0; i < mpiSize; i++)
      for (int j = 0; j < nPartitionElements[i]; j++)
419
	partitionMap[partitionElements[elmdist[i] + j]] = i;
420
421

    pMap = partitionMap;
422
423
  }

424

425
  bool ParMetisPartitioner::distributePartitioning(int *part) 
426
  {
427
428
    FUNCNAME("ParMetisPartitioner::distributePartitioning()");

429
430
    int mpiSize = mpiComm->Get_size();
    int mpiRank = mpiComm->Get_rank();
431
    int nElements = parMetisMesh->getNumElements();
432

433
    // nPartitionElements[i] is the number of elements for the i-th partition
Thomas Witkowski's avatar
Thomas Witkowski committed
434
    int *nPartitionElements = new int[mpiSize];
435
    for (int i = 0; i < mpiSize; i++) 
436
      nPartitionElements[i] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
437
    for (int i = 0; i < nElements; i++)
438
      nPartitionElements[part[i]]++;    
439
440

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

444

445
    // sum up partition elements over all ranks
Thomas Witkowski's avatar
Thomas Witkowski committed
446
    int *sumPartitionElements = new int[mpiSize];
447
    mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize,
448
		       MPI_INT, MPI_SUM);
449
450
451
452
453
454
455
456
457

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

458
    // prepare distribution (fill partitionElements with AMDiS indices)
Thomas Witkowski's avatar
Thomas Witkowski committed
459
    int *bufferOffset = new int[mpiSize];
460
    bufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
461
    for (int i = 1; i < mpiSize; i++)
462
      bufferOffset[i] = bufferOffset[i - 1] + nPartitionElements[i - 1];
463

Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
    int *partitionElements = new int[nElements];
    int **partitionPtr = new int*[mpiSize];
466

467
    for (int i = 0; i < mpiSize; i++)
468
469
      partitionPtr[i] = partitionElements + bufferOffset[i];

470
    sendElements.clear();
471
    for (int i = 0; i < nElements; i++) {
472
      int partition = part[i];
473
      int amdisIndex = parMetisMesh->getAMDiSIndex(i);
474

475
476
477
      if (partition != mpiRank)
	sendElements[partition].push_back(amdisIndex);

478
479
480
481
482
      *(partitionPtr[partition]) = amdisIndex;
      ++(partitionPtr[partition]);
    }

    // all to all: partition elements to rank elements
Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
    int *rankElements = new int[sumPartitionElements[mpiRank]];
    int *recvBufferOffset = new int[mpiSize];
485
    recvBufferOffset[0] = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
486
    for (int i = 1; i < mpiSize; i++)
487
      recvBufferOffset[i] = recvBufferOffset[i - 1] + nRankElements[i - 1];
488

489
    mpiComm->Alltoallv(partitionElements, 
490
		       nPartitionElements,
491
492
493
		       bufferOffset,
		       MPI_INT,
		       rankElements,
494
		       nRankElements,
495
496
		       recvBufferOffset,
		       MPI_INT);
497
    
498
    TEST_EXIT(elementInRank.size() != 0)("Should not happen!\n");
499
    for (map<int, bool>::iterator it = elementInRank.begin();
500
501
	 it != elementInRank.end(); ++it)
      elementInRank[it->first] = false;
502

503
    // Create map which stores for each element index on macro level
504
    // if the element is in the partition of this rank.
505
    recvElements.clear();
506
    for (int i = 0; i < mpiSize; i++) {
507
      int *rankStart = rankElements + recvBufferOffset[i];
508
      int *rankEnd = rankStart + nRankElements[i];
509
510
      for (int *rankPtr = rankStart; rankPtr < rankEnd; ++rankPtr) {
	elementInRank[*rankPtr] = true;
511
512
	if (i != mpiRank)
	  recvElements[i].push_back(*rankPtr);
513
514
515
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
516
    delete parMetisMesh;
517
    parMetisMesh = NULL;
518

Thomas Witkowski's avatar
Thomas Witkowski committed
519
520
521
522
523
524
525
    delete [] rankElements;
    delete [] nPartitionElements;
    delete [] nRankElements;
    delete [] sumPartitionElements;
    delete [] partitionElements;
    delete [] partitionPtr;
    delete [] bufferOffset;
Thomas Witkowski's avatar
huh  
Thomas Witkowski committed
526
    delete [] recvBufferOffset;
527
528

    return true;
529
  }
530
 
531
}