ParMetisPartitioner.cc 13.2 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;
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    int options[4] = {0, 0, 15, PARMETIS_PSR_COUPLED}; // default options
294
    int edgecut = -1;
295
    vector<int> part(nElements);
296

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


    // === Start ParMETIS. ===
308

309
310
    MPI_Comm tmpComm = MPI_Comm(*mpiComm);

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

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

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

376
377
378

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

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

382

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

387
    partitionMap.clear();
388
389

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

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

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

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

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

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

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

    pMap = partitionMap;
423
424
  }

425

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

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

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

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

445

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

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

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

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

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

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

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

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

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

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

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

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

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

    return true;
530
  }
531
 
532
}