MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_RepartitionFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_REPARTITIONFACTORY_DEF_HPP
47 #define MUELU_REPARTITIONFACTORY_DEF_HPP
48 
49 #include <algorithm>
50 #include <iostream>
51 #include <sstream>
52 
53 #include "MueLu_RepartitionFactory_decl.hpp" // TMP JG NOTE: before other includes, otherwise I cannot test the fwd declaration in _def
54 
55 #ifdef HAVE_MPI
57 #include <Teuchos_CommHelpers.hpp>
59 
60 #include <Xpetra_Map.hpp>
61 #include <Xpetra_MapFactory.hpp>
62 #include <Xpetra_MultiVectorFactory.hpp>
63 #include <Xpetra_VectorFactory.hpp>
64 #include <Xpetra_Import.hpp>
65 #include <Xpetra_ImportFactory.hpp>
66 #include <Xpetra_Export.hpp>
67 #include <Xpetra_ExportFactory.hpp>
68 #include <Xpetra_Matrix.hpp>
69 #include <Xpetra_MatrixFactory.hpp>
70 
71 #include "MueLu_Utilities.hpp"
72 
73 #include "MueLu_CloneRepartitionInterface.hpp"
74 
75 #include "MueLu_Level.hpp"
76 #include "MueLu_MasterList.hpp"
77 #include "MueLu_Monitor.hpp"
78 #include "MueLu_PerfUtils.hpp"
79 
80 namespace MueLu {
81 
82  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84  RCP<ParameterList> validParamList = rcp(new ParameterList());
85 
86 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
87  SET_VALID_ENTRY("repartition: print partition distribution");
88  SET_VALID_ENTRY("repartition: remap parts");
89  SET_VALID_ENTRY("repartition: remap num values");
90  SET_VALID_ENTRY("repartition: remap accept partition");
91  SET_VALID_ENTRY("repartition: node repartition level");
92 #undef SET_VALID_ENTRY
93 
94  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
95  validParamList->set< RCP<const FactoryBase> >("number of partitions", Teuchos::null, "Instance of RepartitionHeuristicFactory.");
96  validParamList->set< RCP<const FactoryBase> >("Partition", Teuchos::null, "Factory of the partition");
97 
98  return validParamList;
99  }
100 
101  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
103  Input(currentLevel, "A");
104  Input(currentLevel, "number of partitions");
105  Input(currentLevel, "Partition");
106  }
107 
108  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
110  FactoryMonitor m(*this, "Build", currentLevel);
111 
112  const Teuchos::ParameterList & pL = GetParameterList();
113  // Access parameters here to make sure that we set the parameter entry flag to "used" even in case of short-circuit evaluation.
114  // TODO (JG): I don't really know if we want to do this.
115  bool remapPartitions = pL.get<bool> ("repartition: remap parts");
116 
117  // TODO: We only need a CrsGraph. This class does not have to be templated on Scalar types.
118  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
119  RCP<const Map> rowMap = A->getRowMap();
120  GO indexBase = rowMap->getIndexBase();
121  Xpetra::UnderlyingLib lib = rowMap->lib();
122 
123  RCP<const Teuchos::Comm<int> > origComm = rowMap->getComm();
124  RCP<const Teuchos::Comm<int> > comm = origComm;
125 
126  int myRank = comm->getRank();
127  int numProcs = comm->getSize();
128 
129  RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
130  TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null, Exceptions::RuntimeError, "Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
131  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
132 
134  int numPartitions = Get<int>(currentLevel, "number of partitions");
135 
136  // ======================================================================================================
137  // Construct decomposition vector
138  // ======================================================================================================
139  RCP<GOVector> decomposition = Get<RCP<GOVector> >(currentLevel, "Partition");
140 
141  // check which factory provides "Partition"
142  if(remapPartitions == true && Teuchos::rcp_dynamic_cast<const CloneRepartitionInterface>(GetFactory("Partition")) != Teuchos::null) {
143  // if "Partition" is provided by a CloneRepartitionInterface class we have to switch of remapPartitions
144  // as we can assume the processor Ids in Partition to be the expected ones. If we would do remapping we
145  // would get different processors for the different blocks which screws up matrix-matrix multiplication.
146  remapPartitions = false;
147  }
148 
149  // check special cases
150  if (numPartitions == 1) {
151  // Trivial case: decomposition is the trivial one, all zeros. We skip the call to Zoltan_Interface
152  // (this is mostly done to avoid extra output messages, as even if we didn't skip there is a shortcut
153  // in Zoltan[12]Interface).
154  // TODO: We can probably skip more work in this case (like building all extra data structures)
155  GetOStream(Warnings0) << "Only one partition: Skip call to the repartitioner." << std::endl;
156  } else if (numPartitions == -1) {
157  // No repartitioning necessary: decomposition should be Teuchos::null
158  GetOStream(Warnings0) << "No repartitioning necessary: partitions were left unchanged by the repartitioner" << std::endl;
159  Set<RCP<const Import> >(currentLevel, "Importer", Teuchos::null);
160  return;
161  }
162 
163  // If we're doing node away, we need to be sure to get the mapping to the NodeComm's rank 0.
164  const int nodeRepartLevel = pL.get<int> ("repartition: node repartition level");
165  if(currentLevel.GetLevelID() == nodeRepartLevel) {
166  // NodePartitionInterface returns the *ranks* of the guy who gets the info, not the *partition number*
167  // In a sense, we've already done remap here.
168 
169  // FIXME: We need a low-comm import construction
170  remapPartitions = false;
171  }
172 
173  // ======================================================================================================
174  // Remap if necessary
175  // ======================================================================================================
176  // From a user perspective, we want user to not care about remapping, thinking of it as only a performance feature.
177  // There are two problems, however.
178  // (1) Next level aggregation depends on the order of GIDs in the vector, if one uses "natural" or "random" orderings.
179  // This also means that remapping affects next level aggregation, despite the fact that the _set_ of GIDs for
180  // each partition is the same.
181  // (2) Even with the fixed order of GIDs, the remapping may influence the aggregation for the next-next level.
182  // Let us consider the following example. Lets assume that when we don't do remapping, processor 0 would have
183  // GIDs {0,1,2}, and processor 1 GIDs {3,4,5}, and if we do remapping processor 0 would contain {3,4,5} and
184  // processor 1 {0,1,2}. Now, when we run repartitioning algorithm on the next level (say Zoltan1 RCB), it may
185  // be dependent on whether whether it is [{0,1,2}, {3,4,5}] or [{3,4,5}, {0,1,2}]. Specifically, the tie-breaking
186  // algorithm can resolve these differently. For instance, running
187  // mpirun -np 5 ./MueLu_ScalingTestParamList.exe --xml=easy_sa.xml --nx=12 --ny=12 --nz=12
188  // with
189  // <ParameterList name="MueLu">
190  // <Parameter name="coarse: max size" type="int" value="1"/>
191  // <Parameter name="repartition: enable" type="bool" value="true"/>
192  // <Parameter name="repartition: min rows per proc" type="int" value="2"/>
193  // <ParameterList name="level 1">
194  // <Parameter name="repartition: remap parts" type="bool" value="false/true"/>
195  // </ParameterList>
196  // </ParameterList>
197  // produces different repartitioning for level 2.
198  // This different repartitioning may then escalate into different aggregation for the next level.
199  //
200  // We fix (1) by fixing the order of GIDs in a vector by sorting the resulting vector.
201  // Fixing (2) is more complicated.
202  // FIXME: Fixing (2) in Zoltan may not be enough, as we may use some arbitration in MueLu,
203  // for instance with CoupledAggregation. What we really need to do is to use the same order of processors containing
204  // the same order of GIDs. To achieve that, the newly created subcommunicator must be conforming with the order. For
205  // instance, if we have [{0,1,2}, {3,4,5}], we create a subcommunicator where processor 0 gets rank 0, and processor 1
206  // gets rank 1. If, on the other hand, we have [{3,4,5}, {0,1,2}], we assign rank 1 to processor 0, and rank 0 to processor 1.
207  // This rank permutation requires help from Epetra/Tpetra, both of which have no such API in place.
208  // One should also be concerned that if we had such API in place, rank 0 in subcommunicator may no longer be rank 0 in
209  // MPI_COMM_WORLD, which may lead to issues for logging.
210  if (remapPartitions) {
211  SubFactoryMonitor m1(*this, "DeterminePartitionPlacement", currentLevel);
212 
213  bool acceptPartition = pL.get<bool>("repartition: remap accept partition");
214  bool allSubdomainsAcceptPartitions;
215  int localNumAcceptPartition = acceptPartition;
216  int globalNumAcceptPartition;
217  MueLu_sumAll(comm, localNumAcceptPartition, globalNumAcceptPartition);
218  GetOStream(Statistics2) << "Number of ranks that accept partitions: " << globalNumAcceptPartition << std::endl;
219  if (globalNumAcceptPartition < numPartitions) {
220  GetOStream(Warnings0) << "Not enough ranks are willing to accept a partition, allowing partitions on all ranks." << std::endl;
221  acceptPartition = true;
222  allSubdomainsAcceptPartitions = true;
223  } else if (numPartitions > numProcs) {
224  // We are trying to repartition to a larger communicator.
225  allSubdomainsAcceptPartitions = true;
226  } else {
227  allSubdomainsAcceptPartitions = false;
228  }
229 
230  DeterminePartitionPlacement(*A, *decomposition, numPartitions, acceptPartition, allSubdomainsAcceptPartitions);
231  }
232 
233  // ======================================================================================================
234  // Construct importer
235  // ======================================================================================================
236  // At this point, the following is true:
237  // * Each processors owns 0 or 1 partitions
238  // * If a processor owns a partition, that partition number is equal to the processor rank
239  // * The decomposition vector contains the partitions ids that the corresponding GID belongs to
240 
241  ArrayRCP<const GO> decompEntries;
242  if (decomposition->getLocalLength() > 0)
243  decompEntries = decomposition->getData(0);
244 
245 #ifdef HAVE_MUELU_DEBUG
246  // Test range of partition ids
247  int incorrectRank = -1;
248  for (int i = 0; i < decompEntries.size(); i++)
249  if (decompEntries[i] >= numProcs || decompEntries[i] < 0) {
250  incorrectRank = myRank;
251  break;
252  }
253 
254  int incorrectGlobalRank = -1;
255  MueLu_maxAll(comm, incorrectRank, incorrectGlobalRank);
256  TEUCHOS_TEST_FOR_EXCEPTION(incorrectGlobalRank >- 1, Exceptions::RuntimeError, "pid " + Teuchos::toString(incorrectGlobalRank) + " encountered a partition number is that out-of-range");
257 #endif
258 
259  Array<GO> myGIDs;
260  myGIDs.reserve(decomposition->getLocalLength());
261 
262  // Step 0: Construct mapping
263  // part number -> GIDs I own which belong to this part
264  // NOTE: my own part GIDs are not part of the map
265  typedef std::map<GO, Array<GO> > map_type;
266  map_type sendMap;
267  for (LO i = 0; i < decompEntries.size(); i++) {
268  GO id = decompEntries[i];
269  GO GID = rowMap->getGlobalElement(i);
270 
271  if (id == myRank)
272  myGIDs .push_back(GID);
273  else
274  sendMap[id].push_back(GID);
275  }
276  decompEntries = Teuchos::null;
277 
278  if (IsPrint(Statistics2)) {
279  GO numLocalKept = myGIDs.size(), numGlobalKept, numGlobalRows = A->getGlobalNumRows();
280  MueLu_sumAll(comm,numLocalKept, numGlobalKept);
281  GetOStream(Statistics2) << "Unmoved rows: " << numGlobalKept << " / " << numGlobalRows << " (" << 100*Teuchos::as<double>(numGlobalKept)/numGlobalRows << "%)" << std::endl;
282  }
283 
284  int numSend = sendMap.size(), numRecv;
285 
286  // Arrayify map keys
287  Array<GO> myParts(numSend), myPart(1);
288  int cnt = 0;
289  myPart[0] = myRank;
290  for (typename map_type::const_iterator it = sendMap.begin(); it != sendMap.end(); it++)
291  myParts[cnt++] = it->first;
292 
293  // Step 1: Find out how many processors send me data
294  // partsIndexBase starts from zero, as the processors ids start from zero
295  GO partsIndexBase = 0;
296  RCP<Map> partsIHave = MapFactory ::Build(lib, Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), myParts(), partsIndexBase, comm);
297  RCP<Map> partsIOwn = MapFactory ::Build(lib, numProcs, myPart(), partsIndexBase, comm);
298  RCP<Export> partsExport = ExportFactory::Build(partsIHave, partsIOwn);
299 
300  RCP<GOVector> partsISend = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIHave);
301  RCP<GOVector> numPartsIRecv = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIOwn);
302  if (numSend) {
303  ArrayRCP<GO> partsISendData = partsISend->getDataNonConst(0);
304  for (int i = 0; i < numSend; i++)
305  partsISendData[i] = 1;
306  }
307  (numPartsIRecv->getDataNonConst(0))[0] = 0;
308 
309  numPartsIRecv->doExport(*partsISend, *partsExport, Xpetra::ADD);
310  numRecv = (numPartsIRecv->getData(0))[0];
311 
312  // Step 2: Get my GIDs from everybody else
313  MPI_Datatype MpiType = Teuchos::Details::MpiTypeTraits<GO>::getType();
314  int msgTag = 12345; // TODO: use Comm::dup for all internal messaging
315 
316  // Post sends
317  Array<MPI_Request> sendReqs(numSend);
318  cnt = 0;
319  for (typename map_type::iterator it = sendMap.begin(); it != sendMap.end(); it++)
320  MPI_Isend(static_cast<void*>(it->second.getRawPtr()), it->second.size(), MpiType, Teuchos::as<GO>(it->first), msgTag, *rawMpiComm, &sendReqs[cnt++]);
321 
322  map_type recvMap;
323  size_t totalGIDs = myGIDs.size();
324  for (int i = 0; i < numRecv; i++) {
325  MPI_Status status;
326  MPI_Probe(MPI_ANY_SOURCE, msgTag, *rawMpiComm, &status);
327 
328  // Get rank and number of elements from status
329  int fromRank = status.MPI_SOURCE, count;
330  MPI_Get_count(&status, MpiType, &count);
331 
332  recvMap[fromRank].resize(count);
333  MPI_Recv(static_cast<void*>(recvMap[fromRank].getRawPtr()), count, MpiType, fromRank, msgTag, *rawMpiComm, &status);
334 
335  totalGIDs += count;
336  }
337 
338  // Do waits on send requests
339  if (numSend) {
340  Array<MPI_Status> sendStatuses(numSend);
341  MPI_Waitall(numSend, sendReqs.getRawPtr(), sendStatuses.getRawPtr());
342  }
343 
344  // Merge GIDs
345  myGIDs.reserve(totalGIDs);
346  for (typename map_type::const_iterator it = recvMap.begin(); it != recvMap.end(); it++) {
347  int offset = myGIDs.size(), len = it->second.size();
348  if (len) {
349  myGIDs.resize(offset + len);
350  memcpy(myGIDs.getRawPtr() + offset, it->second.getRawPtr(), len*sizeof(GO));
351  }
352  }
353  // NOTE 2: The general sorting algorithm could be sped up by using the knowledge that original myGIDs and all received chunks
354  // (i.e. it->second) are sorted. Therefore, a merge sort would work well in this situation.
355  std::sort(myGIDs.begin(), myGIDs.end());
356 
357  // Step 3: Construct importer
358  RCP<Map> newRowMap = MapFactory ::Build(lib, rowMap->getGlobalNumElements(), myGIDs(), indexBase, origComm);
359  RCP<const Import> rowMapImporter;
360 
361  RCP<const BlockedMap> blockedRowMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
362 
363  {
364  SubFactoryMonitor m1(*this, "Import construction", currentLevel);
365  // Generate a nonblocked rowmap if we need one
366  if(blockedRowMap.is_null())
367  rowMapImporter = ImportFactory::Build(rowMap, newRowMap);
368  else {
369  rowMapImporter = ImportFactory::Build(blockedRowMap->getMap(), newRowMap);
370  }
371  }
372 
373  // If we're running BlockedCrs we should chop up the newRowMap into a newBlockedRowMap here (and do likewise for importers)
374  if(!blockedRowMap.is_null()) {
375  SubFactoryMonitor m1(*this, "Blocking newRowMap and Importer", currentLevel);
377 
378  // NOTE: This code qualifies as "correct but not particularly performant" If this needs to be sped up, we can probably read data from the existing importer to
379  // build sub-importers rather than generating new ones ex nihilo
380  size_t numBlocks = blockedRowMap->getNumMaps();
381  std::vector<RCP<const Import> > subImports(numBlocks);
382 
383  for(size_t i=0; i<numBlocks; i++) {
384  RCP<const Map> source = blockedRowMap->getMap(i);
385  RCP<const Map> target = blockedTargetMap->getMap(i);
386  subImports[i] = ImportFactory::Build(source,target);
387  }
388  Set(currentLevel,"SubImporters",subImports);
389  }
390 
391 
392  Set(currentLevel, "Importer", rowMapImporter);
393 
394  // ======================================================================================================
395  // Print some data
396  // ======================================================================================================
397  if (!rowMapImporter.is_null() && IsPrint(Statistics2)) {
398  // int oldRank = SetProcRankVerbose(rebalancedAc->getRowMap()->getComm()->getRank());
399  GetOStream(Statistics2) << PerfUtils::PrintImporterInfo(rowMapImporter, "Importer for rebalancing");
400  // SetProcRankVerbose(oldRank);
401  }
402  if (pL.get<bool>("repartition: print partition distribution") && IsPrint(Statistics2)) {
403  // Print the grid of processors
404  GetOStream(Statistics2) << "Partition distribution over cores (ownership is indicated by '+')" << std::endl;
405 
406  char amActive = (myGIDs.size() ? 1 : 0);
407  std::vector<char> areActive(numProcs, 0);
408  MPI_Gather(&amActive, 1, MPI_CHAR, &areActive[0], 1, MPI_CHAR, 0, *rawMpiComm);
409 
410  int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
411  for (int proc = 0; proc < numProcs; proc += rowWidth) {
412  for (int j = 0; j < rowWidth; j++)
413  if (proc + j < numProcs)
414  GetOStream(Statistics2) << (areActive[proc + j] ? "+" : ".");
415  else
416  GetOStream(Statistics2) << " ";
417 
418  GetOStream(Statistics2) << " " << proc << ":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
419  }
420  }
421 
422  } // Build
423 
424  //----------------------------------------------------------------------
425  template<typename T, typename W>
426  struct Triplet {
427  T i, j;
428  W v;
429  };
430  template<typename T, typename W>
431  static bool compareTriplets(const Triplet<T,W>& a, const Triplet<T,W>& b) {
432  return (a.v > b.v); // descending order
433  }
434 
435  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
437  DeterminePartitionPlacement(const Matrix& A, GOVector& decomposition, GO numPartitions, bool willAcceptPartition, bool allSubdomainsAcceptPartitions) const {
438  RCP<const Map> rowMap = A.getRowMap();
439 
440  RCP<const Teuchos::Comm<int> > comm = rowMap->getComm()->duplicate();
441  int numProcs = comm->getSize();
442 
443  RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
444  TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null, Exceptions::RuntimeError, "Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
445  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
446 
447  const Teuchos::ParameterList& pL = GetParameterList();
448 
449  // maxLocal is a constant which determins the number of largest edges which are being exchanged
450  // The idea is that we do not want to construct the full bipartite graph, but simply a subset of
451  // it, which requires less communication. By selecting largest local edges we hope to achieve
452  // similar results but at a lower cost.
453  const int maxLocal = pL.get<int>("repartition: remap num values");
454  const int dataSize = 2*maxLocal;
455 
456  ArrayRCP<GO> decompEntries;
457  if (decomposition.getLocalLength() > 0)
458  decompEntries = decomposition.getDataNonConst(0);
459 
460  // Step 1: Sort local edges by weight
461  // Each edge of a bipartite graph corresponds to a triplet (i, j, v) where
462  // i: processor id that has some piece of part with part_id = j
463  // j: part id
464  // v: weight of the edge
465  // We set edge weights to be the total number of nonzeros in rows on this processor which
466  // correspond to this part_id. The idea is that when we redistribute matrix, this weight
467  // is a good approximation of the amount of data to move.
468  // We use two maps, original which maps a partition id of an edge to the corresponding weight,
469  // and a reverse one, which is necessary to sort by edges.
470  std::map<GO,GO> lEdges;
471  if (willAcceptPartition)
472  for (LO i = 0; i < decompEntries.size(); i++)
473  lEdges[decompEntries[i]] += A.getNumEntriesInLocalRow(i);
474 
475  // Reverse map, so that edges are sorted by weight.
476  // This results in multimap, as we may have edges with the same weight
477  std::multimap<GO,GO> revlEdges;
478  for (typename std::map<GO,GO>::const_iterator it = lEdges.begin(); it != lEdges.end(); it++)
479  revlEdges.insert(std::make_pair(it->second, it->first));
480 
481  // Both lData and gData are arrays of data which we communicate. The data is stored
482  // in pairs, so that data[2*i+0] is the part index, and data[2*i+1] is the corresponding edge weight.
483  // We do not store processor id in data, as we can compute that by looking on the offset in the gData.
484  Array<GO> lData(dataSize, -1), gData(numProcs * dataSize);
485  int numEdges = 0;
486  for (typename std::multimap<GO,GO>::reverse_iterator rit = revlEdges.rbegin(); rit != revlEdges.rend() && numEdges < maxLocal; rit++) {
487  lData[2*numEdges+0] = rit->second; // part id
488  lData[2*numEdges+1] = rit->first; // edge weight
489  numEdges++;
490  }
491 
492  // Step 2: Gather most edges
493  // Each processors contributes maxLocal edges by providing maxLocal pairs <part id, weight>, which is of size dataSize
494  MPI_Datatype MpiType = Teuchos::Details::MpiTypeTraits<GO>::getType();
495  MPI_Allgather(static_cast<void*>(lData.getRawPtr()), dataSize, MpiType, static_cast<void*>(gData.getRawPtr()), dataSize, MpiType, *rawMpiComm);
496 
497  // Step 3: Construct mapping
498 
499  // Construct the set of triplets
500  Teuchos::Array<Triplet<int,int> > gEdges(numProcs * maxLocal);
501  Teuchos::Array<bool> procWillAcceptPartition(numProcs, allSubdomainsAcceptPartitions);
502  size_t k = 0;
503  for (LO i = 0; i < gData.size(); i += 2) {
504  int procNo = i/dataSize; // determine the processor by its offset (since every processor sends the same amount)
505  GO part = gData[i+0];
506  GO weight = gData[i+1];
507  if (part != -1) { // skip nonexistent edges
508  gEdges[k].i = procNo;
509  gEdges[k].j = part;
510  gEdges[k].v = weight;
511  procWillAcceptPartition[procNo] = true;
512  k++;
513  }
514  }
515  gEdges.resize(k);
516 
517  // Sort edges by weight
518  // NOTE: compareTriplets is actually a reverse sort, so the edges weight is in decreasing order
519  std::sort(gEdges.begin(), gEdges.end(), compareTriplets<int,int>);
520 
521  // Do matching
522  std::map<int,int> match;
523  Teuchos::Array<char> matchedRanks(numProcs, 0), matchedParts(numPartitions, 0);
524  int numMatched = 0;
525  for (typename Teuchos::Array<Triplet<int,int> >::const_iterator it = gEdges.begin(); it != gEdges.end(); it++) {
526  GO rank = it->i;
527  GO part = it->j;
528  if (matchedRanks[rank] == 0 && matchedParts[part] == 0) {
529  matchedRanks[rank] = 1;
530  matchedParts[part] = 1;
531  match[part] = rank;
532  numMatched++;
533  }
534  }
535  GetOStream(Statistics1) << "Number of unassigned partitions before cleanup stage: " << (numPartitions - numMatched) << " / " << numPartitions << std::endl;
536 
537  // Step 4: Assign unassigned partitions if necessary.
538  // We do that through desperate matching for remaining partitions:
539  // We select the lowest rank that can still take a partition.
540  // The reason it is done this way is that we don't need any extra communication, as we don't
541  // need to know which parts are valid.
542  if (numPartitions - numMatched > 0) {
543  Teuchos::Array<char> partitionCounts(numPartitions, 0);
544  for (typename std::map<int,int>::const_iterator it = match.begin(); it != match.end(); it++)
545  partitionCounts[it->first] += 1;
546  for (int part = 0, matcher = 0; part < numPartitions; part++) {
547  if (partitionCounts[part] == 0) {
548  // Find first non-matched rank that accepts partitions
549  while (matchedRanks[matcher] || !procWillAcceptPartition[matcher])
550  matcher++;
551 
552  match[part] = matcher++;
553  numMatched++;
554  }
555  }
556  }
557 
558  TEUCHOS_TEST_FOR_EXCEPTION(numMatched != numPartitions, Exceptions::RuntimeError, "MueLu::RepartitionFactory::DeterminePartitionPlacement: Only " << numMatched << " partitions out of " << numPartitions << " got assigned to ranks.");
559 
560  // Step 5: Permute entries in the decomposition vector
561  for (LO i = 0; i < decompEntries.size(); i++)
562  decompEntries[i] = match[decompEntries[i]];
563  }
564 
565 } // namespace MueLu
566 
567 #endif //ifdef HAVE_MPI
568 
569 #endif // MUELU_REPARTITIONFACTORY_DEF_HPP
Important warning messages (one line)
void Build(Level &currentLevel) const
Build an object with this factory.
void reserve(size_type n)
#define MueLu_sumAll(rcpComm, in, out)
#define MueLu_maxAll(rcpComm, in, out)
void DeterminePartitionPlacement(const Matrix &A, GOVector &decomposition, GO numPartitions, bool willAcceptPartition=true, bool allSubdomainsAcceptPartitions=true) const
Determine which process should own each partition.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
Print more statistics.
size_type size() const
#define SET_VALID_ENTRY(name)
static std::string PrintImporterInfo(RCP< const Import > importer, const std::string &msgTag)
Print even more statistics.
static bool compareTriplets(const Triplet< T, W > &a, const Triplet< T, W > &b)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void resize(size_type new_size, const value_type &x=value_type())
iterator end()
void DeclareInput(Level &currentLevel) const
Determines the data that RepartitionFactory needs, and the factories that generate that data...
void push_back(const value_type &x)
size_type size() const
Exception throws to report errors in the internal logical of the program.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
iterator begin()
std::string toString(const T &t)
bool is_null() const