MueLu  Version of the Day
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
56 #include <Teuchos_DefaultMpiComm.hpp>
57 #include <Teuchos_CommHelpers.hpp>
58 
59 #include <Xpetra_Map.hpp>
60 #include <Xpetra_MapFactory.hpp>
61 #include <Xpetra_VectorFactory.hpp>
62 #include <Xpetra_Import.hpp>
63 #include <Xpetra_ImportFactory.hpp>
64 #include <Xpetra_Export.hpp>
65 #include <Xpetra_ExportFactory.hpp>
66 #include <Xpetra_Matrix.hpp>
67 #include <Xpetra_MatrixFactory.hpp>
68 
69 #include "MueLu_Utilities.hpp"
70 
71 #include "MueLu_Level.hpp"
72 #include "MueLu_MasterList.hpp"
73 #include "MueLu_Monitor.hpp"
74 
75 namespace MueLu {
76 
77  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
79  RCP<ParameterList> validParamList = rcp(new ParameterList());
80 
81 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
82  SET_VALID_ENTRY("repartition: start level");
83  SET_VALID_ENTRY("repartition: min rows per proc");
84  SET_VALID_ENTRY("repartition: max imbalance");
85  SET_VALID_ENTRY("repartition: print partition distribution");
86  SET_VALID_ENTRY("repartition: remap parts");
87  SET_VALID_ENTRY("repartition: remap num values");
88 #undef SET_VALID_ENTRY
89 
90  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory of the matrix A");
91  validParamList->set< RCP<const FactoryBase> >("Partition", Teuchos::null, "Factory of the partition");
92 
93  return validParamList;
94  }
95 
96  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98  Input(currentLevel, "A");
99  Input(currentLevel, "Partition");
100  }
101 
102  template<class T> class MpiTypeTraits { public: static MPI_Datatype getType(); };
103  template<> class MpiTypeTraits<long> { public: static MPI_Datatype getType() { return MPI_LONG; } };
104  template<> class MpiTypeTraits<int> { public: static MPI_Datatype getType() { return MPI_INT; } };
105  template<> class MpiTypeTraits<short> { public: static MPI_Datatype getType() { return MPI_SHORT; } };
106  template<> class MpiTypeTraits<unsigned> { public: static MPI_Datatype getType() { return MPI_UNSIGNED; } };
107  template<> class MpiTypeTraits<long long> { public: static MPI_Datatype getType() { return MPI_LONG_LONG; } };
108 
109  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
111  FactoryMonitor m(*this, "Build", currentLevel);
112 
113  const Teuchos::ParameterList & pL = GetParameterList();
114  // Access parameters here to make sure that we set the parameter entry flag to "used" even in case of short-circuit evaluation.
115  // TODO (JG): I don't really know if we want to do this.
116  const int startLevel = pL.get<int> ("repartition: start level");
117  const LO minRowsPerProcessor = pL.get<LO> ("repartition: min rows per proc");
118  const double nonzeroImbalance = pL.get<double>("repartition: max imbalance");
119  const bool remapPartitions = pL.get<bool> ("repartition: remap parts");
120 
121  // TODO: We only need a CrsGraph. This class does not have to be templated on Scalar types.
122  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
123 
124  // ======================================================================================================
125  // Determine whether partitioning is needed
126  // ======================================================================================================
127  // NOTE: most tests include some global communication, which is why we currently only do tests until we make
128  // a decision on whether to repartition. However, there is value in knowing how "close" we are to having to
129  // rebalance an operator. So, it would probably be beneficial to do and report *all* tests.
130 
131  // Test1: skip repartitioning if current level is less than the specified minimum level for repartitioning
132  if (currentLevel.GetLevelID() < startLevel) {
133  GetOStream(Statistics0) << "Repartitioning? NO:" <<
134  "\n current level = " << Teuchos::toString(currentLevel.GetLevelID()) <<
135  ", first level where repartitioning can happen is " + Teuchos::toString(startLevel) << std::endl;
136 
137  Set<RCP<const Import> >(currentLevel, "Importer", Teuchos::null);
138  return;
139  }
140 
141  RCP<const Map> rowMap = A->getRowMap();
142 
143  // NOTE: Teuchos::MPIComm::duplicate() calls MPI_Bcast inside, so this is
144  // a synchronization point. However, as we do MueLu_sumAll afterwards anyway, it
145  // does not matter.
146  RCP<const Teuchos::Comm<int> > origComm = rowMap->getComm();
147  RCP<const Teuchos::Comm<int> > comm = origComm->duplicate();
148 
149  // Test 2: check whether A is actually distributed, i.e. more than one processor owns part of A
150  // TODO: this global communication can be avoided if we store the information with the matrix (it is known when matrix is created)
151  // TODO: further improvements could be achieved when we use subcommunicator for the active set. Then we only need to check its size
152  {
153  int numActiveProcesses = 0;
154  MueLu_sumAll(comm, Teuchos::as<int>((A->getNodeNumRows() > 0) ? 1 : 0), numActiveProcesses);
155 
156  if (numActiveProcesses == 1) {
157  GetOStream(Statistics0) << "Repartitioning? NO:" <<
158  "\n # processes with rows = " << Teuchos::toString(numActiveProcesses) << std::endl;
159 
160  Set<RCP<const Import> >(currentLevel, "Importer", Teuchos::null);
161  return;
162  }
163  }
164 
165  bool test3 = false, test4 = false;
166  std::string msg3, msg4;
167 
168  // Test3: check whether number of rows on any processor satisfies the minimum number of rows requirement
169  // NOTE: Test2 ensures that repartitionning is not done when there is only one processor (it may or may not satisfy Test3)
170  if (minRowsPerProcessor > 0) {
171  LO numMyRows = Teuchos::as<LO>(A->getNodeNumRows()), minNumRows, LOMAX = Teuchos::OrdinalTraits<LO>::max();
172  LO haveFewRows = (numMyRows < minRowsPerProcessor ? 1 : 0), numWithFewRows = 0;
173  MueLu_sumAll(comm, haveFewRows, numWithFewRows);
174  MueLu_minAll(comm, (numMyRows > 0 ? numMyRows : LOMAX), minNumRows);
175 
176  // TODO: we could change it to repartition only if the number of processors with numRows < minNumRows is larger than some
177  // percentage of the total number. This way, we won't repartition if 2 out of 1000 processors don't have enough elements.
178  // I'm thinking maybe 20% threshold. To implement, simply add " && numWithFewRows < .2*numProcs" to the if statement.
179  if (numWithFewRows > 0)
180  test3 = true;
181 
182  msg3 = "\n min # rows per proc = " + Teuchos::toString(minNumRows) + ", min allowable = " + Teuchos::toString(minRowsPerProcessor);
183  }
184 
185  // Test4: check whether the balance in the number of nonzeros per processor is greater than threshold
186  if (!test3) {
187  GO minNnz, maxNnz, numMyNnz = Teuchos::as<GO>(A->getNodeNumEntries());
188  MueLu_maxAll(comm, numMyNnz, maxNnz);
189  MueLu_minAll(comm, (numMyNnz > 0 ? numMyNnz : maxNnz), minNnz); // min nnz over all active processors
190  double imbalance = Teuchos::as<double>(maxNnz)/minNnz;
191 
192  if (imbalance > nonzeroImbalance)
193  test4 = true;
194 
195  msg4 = "\n nonzero imbalance = " + Teuchos::toString(imbalance) + ", max allowable = " + Teuchos::toString(nonzeroImbalance);
196  }
197 
198  if (!test3 && !test4) {
199  GetOStream(Statistics0) << "Repartitioning? NO:" << msg3 + msg4 << std::endl;
200 
201  Set<RCP<const Import> >(currentLevel, "Importer", Teuchos::null);
202  return;
203  }
204 
205  GetOStream(Statistics0) << "Repartitioning? YES:" << msg3 + msg4 << std::endl;
206 
207  GO indexBase = rowMap->getIndexBase();
208  Xpetra::UnderlyingLib lib = rowMap->lib();
209  int myRank = comm->getRank();
210  int numProcs = comm->getSize();
211 
212  RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
213  TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null, Exceptions::RuntimeError, "Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
214  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
215 
216  // ======================================================================================================
217  // Calculate number of partitions
218  // ======================================================================================================
219  // FIXME Quick way to figure out how many partitions there should be (same algorithm as ML)
220  // FIXME Should take into account nnz? Perhaps only when user is using min #nnz per row threshold.
221  GO numPartitions;
222  if (currentLevel.IsAvailable("number of partitions")) {
223  numPartitions = currentLevel.Get<GO>("number of partitions");
224  GetOStream(Warnings0) << "Using user-provided \"number of partitions\", the performance is unknown" << std::endl;
225 
226  } else {
227  if (Teuchos::as<GO>(A->getGlobalNumRows()) < minRowsPerProcessor) {
228  // System is too small, migrate it to a single processor
229  numPartitions = 1;
230 
231  } else {
232  // Make sure that each processor has approximately minRowsPerProcessor
233  numPartitions = A->getGlobalNumRows() / minRowsPerProcessor;
234  }
235  numPartitions = std::min(numPartitions, Teuchos::as<GO>(numProcs));
236 
237  currentLevel.Set("number of partitions", numPartitions, NoFactory::get());
238  }
239  GetOStream(Statistics0) << "Number of partitions to use = " << numPartitions << std::endl;
240 
241  // ======================================================================================================
242  // Construct decomposition vector
243  // ======================================================================================================
244  RCP<GOVector> decomposition;
245  if (numPartitions == 1) {
246  // Trivial case: decomposition is the trivial one, all zeros. We skip the call to Zoltan_Interface
247  // (this is mostly done to avoid extra output messages, as even if we didn't skip there is a shortcut
248  // in Zoltan[12]Interface).
249  // TODO: We can probably skip more work in this case (like building all extra data structures)
250  GetOStream(Warnings0) << "Only one partition: Skip call to the repartitioner." << std::endl;
251  decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(A->getRowMap(), true);
252 
253  } else {
254  decomposition = Get<RCP<GOVector> >(currentLevel, "Partition");
255 
256  if (decomposition.is_null()) {
257  GetOStream(Warnings0) << "No repartitioning necessary: partitions were left unchanged by the repartitioner" << std::endl;
258  Set<RCP<const Import> >(currentLevel, "Importer", Teuchos::null);
259  return;
260  }
261  }
262 
263  // ======================================================================================================
264  // Remap if necessary
265  // ======================================================================================================
266  // From a user perspective, we want user to not care about remapping, thinking of it as only a performance feature.
267  // There are two problems, however.
268  // (1) Next level aggregation depends on the order of GIDs in the vector, if one uses "natural" or "random" orderings.
269  // This also means that remapping affects next level aggregation, despite the fact that the _set_ of GIDs for
270  // each partition is the same.
271  // (2) Even with the fixed order of GIDs, the remapping may influence the aggregation for the next-next level.
272  // Let us consider the following example. Lets assume that when we don't do remapping, processor 0 would have
273  // GIDs {0,1,2}, and processor 1 GIDs {3,4,5}, and if we do remapping processor 0 would contain {3,4,5} and
274  // processor 1 {0,1,2}. Now, when we run repartitioning algorithm on the next level (say Zoltan1 RCB), it may
275  // be dependent on whether whether it is [{0,1,2}, {3,4,5}] or [{3,4,5}, {0,1,2}]. Specifically, the tie-breaking
276  // algorithm can resolve these differently. For instance, running
277  // mpirun -np 5 ./MueLu_ScalingTestParamList.exe --xml=easy_sa.xml --nx=12 --ny=12 --nz=12
278  // with
279  // <ParameterList name="MueLu">
280  // <Parameter name="coarse: max size" type="int" value="1"/>
281  // <Parameter name="repartition: enable" type="bool" value="true"/>
282  // <Parameter name="repartition: min rows per proc" type="int" value="2"/>
283  // <ParameterList name="level 1">
284  // <Parameter name="repartition: remap parts" type="bool" value="false/true"/>
285  // </ParameterList>
286  // </ParameterList>
287  // produces different repartitioning for level 2.
288  // This different repartitioning may then escalate into different aggregation for the next level.
289  //
290  // We fix (1) by fixing the order of GIDs in a vector by sorting the resulting vector.
291  // Fixing (2) is more complicated.
292  // FIXME: Fixing (2) in Zoltan may not be enough, as we may use some arbitration in MueLu,
293  // for instance with CoupledAggregation. What we really need to do is to use the same order of processors containing
294  // the same order of GIDs. To achieve that, the newly created subcommunicator must be conforming with the order. For
295  // instance, if we have [{0,1,2}, {3,4,5}], we create a subcommunicator where processor 0 gets rank 0, and processor 1
296  // 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.
297  // This rank permutation requires help from Epetra/Tpetra, both of which have no such API in place.
298  // One should also be concerned that if we had such API in place, rank 0 in subcommunicator may no longer be rank 0 in
299  // MPI_COMM_WORLD, which may lead to issues for logging.
300  if (remapPartitions) {
301  SubFactoryMonitor m1(*this, "DeterminePartitionPlacement", currentLevel);
302 
303  DeterminePartitionPlacement(*A, *decomposition, numPartitions);
304  }
305 
306  // ======================================================================================================
307  // Construct importer
308  // ======================================================================================================
309  // At this point, the following is true:
310  // * Each processors owns 0 or 1 partitions
311  // * If a processor owns a partition, that partition number is equal to the processor rank
312  // * The decomposition vector contains the partitions ids that the corresponding GID belongs to
313 
314  ArrayRCP<const GO> decompEntries;
315  if (decomposition->getLocalLength() > 0)
316  decompEntries = decomposition->getData(0);
317 
318 #ifdef HAVE_MUELU_DEBUG
319  // Test range of partition ids
320  int incorrectRank = -1;
321  for (int i = 0; i < decompEntries.size(); i++)
322  if (decompEntries[i] >= numProcs || decompEntries[i] < 0) {
323  incorrectRank = myRank;
324  break;
325  }
326 
327  int incorrectGlobalRank = -1;
328  MueLu_maxAll(comm, incorrectRank, incorrectGlobalRank);
329  TEUCHOS_TEST_FOR_EXCEPTION(incorrectGlobalRank >- 1, Exceptions::RuntimeError, "pid " + Teuchos::toString(incorrectGlobalRank) + " encountered a partition number is that out-of-range");
330 #endif
331 
332  Array<GO> myGIDs;
333  myGIDs.reserve(decomposition->getLocalLength());
334 
335  // Step 0: Construct mapping
336  // part number -> GIDs I own which belong to this part
337  // NOTE: my own part GIDs are not part of the map
338  typedef std::map<GO, Array<GO> > map_type;
339  map_type sendMap;
340  for (LO i = 0; i < decompEntries.size(); i++) {
341  GO id = decompEntries[i];
342  GO GID = rowMap->getGlobalElement(i);
343 
344  if (id == myRank)
345  myGIDs .push_back(GID);
346  else
347  sendMap[id].push_back(GID);
348  }
349  decompEntries = Teuchos::null;
350 
351  if (IsPrint(Statistics2)) {
352  GO numLocalKept = myGIDs.size(), numGlobalKept, numGlobalRows = A->getGlobalNumRows();
353  MueLu_sumAll(comm,numLocalKept, numGlobalKept);
354  GetOStream(Statistics2) << "Unmoved rows: " << numGlobalKept << " / " << numGlobalRows << " (" << 100*Teuchos::as<double>(numGlobalKept)/numGlobalRows << "%)" << std::endl;
355  }
356 
357  int numSend = sendMap.size(), numRecv;
358 
359  // Arrayify map keys
360  Array<GO> myParts(numSend), myPart(1);
361  int cnt = 0;
362  myPart[0] = myRank;
363  for (typename map_type::const_iterator it = sendMap.begin(); it != sendMap.end(); it++)
364  myParts[cnt++] = it->first;
365 
366  // Step 1: Find out how many processors send me data
367  // partsIndexBase starts from zero, as the processors ids start from zero
368  GO partsIndexBase = 0;
369  RCP<Map> partsIHave = MapFactory ::Build(lib, Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), myParts(), partsIndexBase, comm);
370  RCP<Map> partsIOwn = MapFactory ::Build(lib, numProcs, myPart(), partsIndexBase, comm);
371  RCP<Export> partsExport = ExportFactory::Build(partsIHave, partsIOwn);
372 
373  RCP<GOVector> partsISend = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIHave);
374  RCP<GOVector> numPartsIRecv = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIOwn);
375  if (numSend) {
376  ArrayRCP<GO> partsISendData = partsISend->getDataNonConst(0);
377  for (int i = 0; i < numSend; i++)
378  partsISendData[i] = 1;
379  }
380  (numPartsIRecv->getDataNonConst(0))[0] = 0;
381 
382  numPartsIRecv->doExport(*partsISend, *partsExport, Xpetra::ADD);
383  numRecv = (numPartsIRecv->getData(0))[0];
384 
385  // Step 2: Get my GIDs from everybody else
386  MPI_Datatype MpiType = MpiTypeTraits<GO>::getType();
387  int msgTag = 12345; // TODO: use Comm::dup for all internal messaging
388 
389  // Post sends
390  Array<MPI_Request> sendReqs(numSend);
391  cnt = 0;
392  for (typename map_type::iterator it = sendMap.begin(); it != sendMap.end(); it++)
393  MPI_Isend(static_cast<void*>(it->second.getRawPtr()), it->second.size(), MpiType, Teuchos::as<GO>(it->first), msgTag, *rawMpiComm, &sendReqs[cnt++]);
394 
395  map_type recvMap;
396  size_t totalGIDs = myGIDs.size();
397  for (int i = 0; i < numRecv; i++) {
398  MPI_Status status;
399  MPI_Probe(MPI_ANY_SOURCE, msgTag, *rawMpiComm, &status);
400 
401  // Get rank and number of elements from status
402  int fromRank = status.MPI_SOURCE, count;
403  MPI_Get_count(&status, MpiType, &count);
404 
405  recvMap[fromRank].resize(count);
406  MPI_Recv(static_cast<void*>(recvMap[fromRank].getRawPtr()), count, MpiType, fromRank, msgTag, *rawMpiComm, &status);
407 
408  totalGIDs += count;
409  }
410 
411  // Do waits on send requests
412  if (numSend) {
413  Array<MPI_Status> sendStatuses(numSend);
414  MPI_Waitall(numSend, sendReqs.getRawPtr(), sendStatuses.getRawPtr());
415  }
416 
417  // Merge GIDs
418  myGIDs.reserve(totalGIDs);
419  for (typename map_type::const_iterator it = recvMap.begin(); it != recvMap.end(); it++) {
420  int offset = myGIDs.size(), len = it->second.size();
421  if (len) {
422  myGIDs.resize(offset + len);
423  memcpy(myGIDs.getRawPtr() + offset, it->second.getRawPtr(), len*sizeof(GO));
424  }
425  }
426  // NOTE 2: The general sorting algorithm could be sped up by using the knowledge that original myGIDs and all received chunks
427  // (i.e. it->second) are sorted. Therefore, a merge sort would work well in this situation.
428  std::sort(myGIDs.begin(), myGIDs.end());
429 
430  // Step 3: Construct importer
431  RCP<Map> newRowMap = MapFactory ::Build(lib, rowMap->getGlobalNumElements(), myGIDs(), indexBase, origComm);
432  RCP<const Import> rowMapImporter;
433  {
434  SubFactoryMonitor m1(*this, "Import construction", currentLevel);
435  rowMapImporter = ImportFactory::Build(rowMap, newRowMap);
436  }
437 
438  Set(currentLevel, "Importer", rowMapImporter);
439 
440  // ======================================================================================================
441  // Print some data
442  // ======================================================================================================
443  if (pL.get<bool>("repartition: print partition distribution") && IsPrint(Statistics2)) {
444  // Print the grid of processors
445  GetOStream(Statistics2) << "Partition distribution over cores (ownership is indicated by '+')" << std::endl;
446 
447  char amActive = (myGIDs.size() ? 1 : 0);
448  std::vector<char> areActive(numProcs, 0);
449  MPI_Gather(&amActive, 1, MPI_CHAR, &areActive[0], 1, MPI_CHAR, 0, *rawMpiComm);
450 
451  int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
452  for (int proc = 0; proc < numProcs; proc += rowWidth) {
453  for (int j = 0; j < rowWidth; j++)
454  if (proc + j < numProcs)
455  GetOStream(Statistics2) << (areActive[proc + j] ? "+" : ".");
456  else
457  GetOStream(Statistics2) << " ";
458 
459  GetOStream(Statistics2) << " " << proc << ":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
460  }
461  }
462 
463  } // Build
464 
465  //----------------------------------------------------------------------
466  template<typename T, typename W>
467  struct Triplet {
468  T i, j;
469  W v;
470  };
471  template<typename T, typename W>
472  static bool compareTriplets(const Triplet<T,W>& a, const Triplet<T,W>& b) {
473  return (a.v > b.v); // descending order
474  }
475 
476  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
478  DeterminePartitionPlacement(const Matrix& A, GOVector& decomposition, GO numPartitions) const {
479  RCP<const Map> rowMap = A.getRowMap();
480 
481  RCP<const Teuchos::Comm<int> > comm = rowMap->getComm()->duplicate();
482  int numProcs = comm->getSize();
483 
484  RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
485  TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null, Exceptions::RuntimeError, "Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
486  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
487 
488  const Teuchos::ParameterList& pL = GetParameterList();
489 
490  // maxLocal is a constant which determins the number of largest edges which are being exchanged
491  // The idea is that we do not want to construct the full bipartite graph, but simply a subset of
492  // it, which requires less communication. By selecting largest local edges we hope to achieve
493  // similar results but at a lower cost.
494  const int maxLocal = pL.get<int>("repartition: remap num values");
495  const int dataSize = 2*maxLocal;
496 
497  ArrayRCP<GO> decompEntries;
498  if (decomposition.getLocalLength() > 0)
499  decompEntries = decomposition.getDataNonConst(0);
500 
501  // Step 1: Sort local edges by weight
502  // Each edge of a bipartite graph corresponds to a triplet (i, j, v) where
503  // i: processor id that has some piece of part with part_id = j
504  // j: part id
505  // v: weight of the edge
506  // We set edge weights to be the total number of nonzeros in rows on this processor which
507  // correspond to this part_id. The idea is that when we redistribute matrix, this weight
508  // is a good approximation of the amount of data to move.
509  // We use two maps, original which maps a partition id of an edge to the corresponding weight,
510  // and a reverse one, which is necessary to sort by edges.
511  std::map<GO,GO> lEdges;
512  for (LO i = 0; i < decompEntries.size(); i++)
513  lEdges[decompEntries[i]] += A.getNumEntriesInLocalRow(i);
514 
515  // Reverse map, so that edges are sorted by weight.
516  // This results in multimap, as we may have edges with the same weight
517  std::multimap<GO,GO> revlEdges;
518  for (typename std::map<GO,GO>::const_iterator it = lEdges.begin(); it != lEdges.end(); it++)
519  revlEdges.insert(std::make_pair(it->second, it->first));
520 
521  // Both lData and gData are arrays of data which we communicate. The data is stored
522  // in pairs, so that data[2*i+0] is the part index, and data[2*i+1] is the corresponding edge weight.
523  // We do not store processor id in data, as we can compute that by looking on the offset in the gData.
524  Array<GO> lData(dataSize, -1), gData(numProcs * dataSize);
525  int numEdges = 0;
526  for (typename std::multimap<GO,GO>::reverse_iterator rit = revlEdges.rbegin(); rit != revlEdges.rend() && numEdges < maxLocal; rit++) {
527  lData[2*numEdges+0] = rit->second; // part id
528  lData[2*numEdges+1] = rit->first; // edge weight
529  numEdges++;
530  }
531 
532  // Step 2: Gather most edges
533  // Each processors contributes maxLocal edges by providing maxLocal pairs <part id, weight>, which is of size dataSize
534  MPI_Datatype MpiType = MpiTypeTraits<GO>::getType();
535  MPI_Allgather(static_cast<void*>(lData.getRawPtr()), dataSize, MpiType, static_cast<void*>(gData.getRawPtr()), dataSize, MpiType, *rawMpiComm);
536 
537  // Step 3: Construct mapping
538 
539  // Construct the set of triplets
540  std::vector<Triplet<int,int> > gEdges(numProcs * maxLocal);
541  size_t k = 0;
542  for (LO i = 0; i < gData.size(); i += 2) {
543  GO part = gData[i+0];
544  GO weight = gData[i+1];
545  if (part != -1) { // skip nonexistent edges
546  gEdges[k].i = i/dataSize; // determine the processor by its offset (since every processor sends the same amount)
547  gEdges[k].j = part;
548  gEdges[k].v = weight;
549  k++;
550  }
551  }
552  gEdges.resize(k);
553 
554  // Sort edges by weight
555  // NOTE: compareTriplets is actually a reverse sort, so the edges weight is in decreasing order
556  std::sort(gEdges.begin(), gEdges.end(), compareTriplets<int,int>);
557 
558  // Do matching
559  std::map<int,int> match;
560  std::vector<char> matchedRanks(numProcs, 0), matchedParts(numProcs, 0);
561  int numMatched = 0;
562  for (typename std::vector<Triplet<int,int> >::const_iterator it = gEdges.begin(); it != gEdges.end(); it++) {
563  GO rank = it->i;
564  GO part = it->j;
565  if (matchedRanks[rank] == 0 && matchedParts[part] == 0) {
566  matchedRanks[rank] = 1;
567  matchedParts[part] = 1;
568  match[part] = rank;
569  numMatched++;
570  }
571  }
572  GetOStream(Statistics0) << "Number of unassigned paritions before cleanup stage: " << (numPartitions - numMatched) << " / " << numPartitions << std::endl;
573 
574  // Step 4: Assign unassigned partitions
575  // We do that through random matching for remaining partitions. Not all part numbers are valid, but valid parts are a subset of [0, numProcs).
576  // The reason it is done this way is that we don't need any extra communication, as we don't need to know which parts are valid.
577  for (int part = 0, matcher = 0; part < numProcs; part++)
578  if (match.count(part) == 0) {
579  // Find first non-matched rank
580  while (matchedRanks[matcher])
581  matcher++;
582 
583  match[part] = matcher++;
584  }
585 
586  // Step 5: Permute entries in the decomposition vector
587  for (LO i = 0; i < decompEntries.size(); i++)
588  decompEntries[i] = match[decompEntries[i]];
589  }
590 
591 } // namespace MueLu
592 
593 #endif //ifdef HAVE_MPI
594 
595 #endif // MUELU_REPARTITIONFACTORY_DEF_HPP
Important warning messages (one line)
void Build(Level &currentLevel) const
Build an object with this factory.
#define MueLu_sumAll(rcpComm, in, out)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
#define MueLu_maxAll(rcpComm, in, out)
Timer to be used in factories. Similar to Monitor but with additional timers.
static MPI_Datatype getType()
Namespace for MueLu classes and methods.
#define SET_VALID_ENTRY(name)
#define MueLu_minAll(rcpComm, in, out)
static const NoFactory * get()
Print even more statistics.
Print statistics that do not involve significant additional computation.
static bool compareTriplets(const Triplet< T, W > &a, const Triplet< T, W > &b)
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 DeclareInput(Level &currentLevel) const
Determines the data that RepartitionFactory needs, and the factories that generate that data...
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.
void DeterminePartitionPlacement(const Matrix &A, GOVector &decomposition, GO numPartitions) const
Determine which process should own each partition.