MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_RebalanceTransferFactory_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_REBALANCETRANSFERFACTORY_DEF_HPP
47 #define MUELU_REBALANCETRANSFERFACTORY_DEF_HPP
48 
49 #include <sstream>
50 #include <Teuchos_Tuple.hpp>
51 
52 #include "Xpetra_MultiVector.hpp"
53 #include "Xpetra_MultiVectorFactory.hpp"
54 #include "Xpetra_Vector.hpp"
55 #include "Xpetra_VectorFactory.hpp"
56 #include <Xpetra_Matrix.hpp>
57 #include <Xpetra_MapFactory.hpp>
58 #include <Xpetra_MatrixFactory.hpp>
59 #include <Xpetra_Import.hpp>
60 #include <Xpetra_ImportFactory.hpp>
61 #include <Xpetra_IO.hpp>
62 
64 
65 #include "MueLu_Level.hpp"
66 #include "MueLu_MasterList.hpp"
67 #include "MueLu_Monitor.hpp"
68 #include "MueLu_PerfUtils.hpp"
69 
70 namespace MueLu {
71 
72 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74  RCP<ParameterList> validParamList = rcp(new ParameterList());
75 
76 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
77  SET_VALID_ENTRY("repartition: rebalance P and R");
78  SET_VALID_ENTRY("repartition: explicit via new copy rebalance P and R");
79  SET_VALID_ENTRY("repartition: rebalance Nullspace");
80  SET_VALID_ENTRY("transpose: use implicit");
81  SET_VALID_ENTRY("repartition: use subcommunicators");
82 #undef SET_VALID_ENTRY
83 
84  {
86  RCP<validatorType> typeValidator = rcp(new validatorType(Teuchos::tuple<std::string>("Interpolation", "Restriction"), "type"));
87  validParamList->set("type", "Interpolation", "Type of the transfer operator that need to be rebalanced (Interpolation or Restriction)", typeValidator);
88  }
89 
90  validParamList->set<RCP<const FactoryBase> >("P", null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
91  validParamList->set<RCP<const FactoryBase> >("R", null, "Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
92  validParamList->set<RCP<const FactoryBase> >("Nullspace", null, "Factory of the nullspace that need to be rebalanced (only used if type=Interpolation)");
93  validParamList->set<RCP<const FactoryBase> >("Coordinates", null, "Factory of the coordinates that need to be rebalanced (only used if type=Interpolation)");
94  validParamList->set<RCP<const FactoryBase> >("BlockNumber", null, "Factory of the block ids that need to be rebalanced (only used if type=Interpolation)");
95  validParamList->set<RCP<const FactoryBase> >("Importer", null, "Factory of the importer object used for the rebalancing");
96  validParamList->set<int>("write start", -1, "First level at which coordinates should be written to file");
97  validParamList->set<int>("write end", -1, "Last level at which coordinates should be written to file");
98 
99  // TODO validation: "P" parameter valid only for type="Interpolation" and "R" valid only for type="Restriction". Like so:
100  // if (paramList.isEntry("type") && paramList.get("type) == "Interpolation) {
101  // validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
102 
103  return validParamList;
104 }
105 
106 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
108  const ParameterList& pL = GetParameterList();
109 
110  if (pL.get<std::string>("type") == "Interpolation") {
111  Input(coarseLevel, "P");
112  if (pL.get<bool>("repartition: rebalance Nullspace"))
113  Input(coarseLevel, "Nullspace");
114  if (pL.get<RCP<const FactoryBase> >("Coordinates") != Teuchos::null)
115  Input(coarseLevel, "Coordinates");
116  if (pL.get<RCP<const FactoryBase> >("BlockNumber") != Teuchos::null)
117  Input(coarseLevel, "BlockNumber");
118 
119  } else {
120  if (pL.get<bool>("transpose: use implicit") == false)
121  Input(coarseLevel, "R");
122  }
123 
124  Input(coarseLevel, "Importer");
125 }
126 
127 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
129  FactoryMonitor m(*this, "Build", coarseLevel);
131 
132  const ParameterList& pL = GetParameterList();
133 
134  RCP<Matrix> originalP = Get<RCP<Matrix> >(coarseLevel, "P");
135  // If we don't have a valid P (e.g., # global aggregates is 0), skip this rebalancing. This level will
136  // ultimately be removed in MueLu_Hierarchy_defs.h via a resize()
137  if (originalP == Teuchos::null) {
138  Set(coarseLevel, "P", originalP);
139  return;
140  }
141  int implicit = !pL.get<bool>("repartition: rebalance P and R");
142  int reallyExplicit = pL.get<bool>("repartition: explicit via new copy rebalance P and R");
143  int writeStart = pL.get<int>("write start");
144  int writeEnd = pL.get<int>("write end");
145 
146  if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel, "Coordinates")) {
147  std::string fileName = "coordinates_level_0.m";
148  RCP<xdMV> fineCoords = fineLevel.Get<RCP<xdMV> >("Coordinates");
149  if (fineCoords != Teuchos::null)
150  Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO>::Write(fileName, *fineCoords);
151  }
152 
153  if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel, "BlockNumber")) {
154  std::string fileName = "BlockNumber_level_0.m";
155  RCP<LocalOrdinalVector> fineBlockNumber = fineLevel.Get<RCP<LocalOrdinalVector> >("BlockNumber");
156  if (fineBlockNumber != Teuchos::null)
157  Xpetra::IO<LO, LO, GO, NO>::Write(fileName, *fineBlockNumber);
158  }
159 
160  RCP<const Import> importer = Get<RCP<const Import> >(coarseLevel, "Importer");
161  if (implicit) {
162  // Save the importer, we'll need it for solve
163  coarseLevel.Set("Importer", importer, NoFactory::get());
164  }
165 
166  RCP<ParameterList> params = rcp(new ParameterList());
167  if (IsPrint(Statistics2)) {
168  params->set("printLoadBalancingInfo", true);
169  params->set("printCommInfo", true);
170  }
171 
172  std::string transferType = pL.get<std::string>("type");
173  if (transferType == "Interpolation") {
174  originalP = Get<RCP<Matrix> >(coarseLevel, "P");
175 
176  {
177  // This line must be after the Get call
178  SubFactoryMonitor m1(*this, "Rebalancing prolongator", coarseLevel);
179 
180  if (implicit || importer.is_null()) {
181  GetOStream(Runtime0) << "Using original prolongator" << std::endl;
182  Set(coarseLevel, "P", originalP);
183 
184  } else {
185  // There are two version of an explicit rebalanced P and R.
186  // The !reallyExplicit way, is sufficient for all MueLu purposes
187  // with the exception of the CombinePFactory that needs true domain
188  // and column maps.
189  // !reallyExplicit:
190  // Rather than calling fillComplete (which would entail creating a new
191  // column map), it's sufficient to replace the domain map and importer.
192  // Note that this potentially violates the assumption that in the
193  // column map, local IDs appear before any off-rank IDs.
194  //
195  // reallyExplicit:
196  // P transfers from coarse grid to the fine grid. Here, we change
197  // the domain map (coarse) of Paccording to the new partition. The
198  // range map (fine) is kept unchanged.
199  //
200  // The domain map of P must match the range map of R. To change the
201  // domain map of P, P needs to be fillCompleted again with the new
202  // domain map. To achieve this, P is copied into a new matrix that
203  // is not fill-completed. The doImport() operation is just used
204  // here to make a copy of P: the importer is trivial and there is
205  // no data movement involved. The reordering actually happens during
206  // fillComplete() with domainMap == importer->getTargetMap().
207 
208  RCP<Matrix> rebalancedP;
209  if (reallyExplicit) {
210  size_t totalMaxPerRow = 0;
211  ArrayRCP<size_t> nnzPerRow(originalP->getRowMap()->getLocalNumElements(), 0);
212  for (size_t i = 0; i < originalP->getRowMap()->getLocalNumElements(); ++i) {
213  nnzPerRow[i] = originalP->getNumEntriesInLocalRow(i);
214  if (nnzPerRow[i] > totalMaxPerRow) totalMaxPerRow = nnzPerRow[i];
215  }
216 
217  rebalancedP = MatrixFactory::Build(originalP->getRowMap(), totalMaxPerRow);
218 
219  {
220  RCP<Import> trivialImporter = ImportFactory::Build(originalP->getRowMap(), originalP->getRowMap());
221  SubFactoryMonitor m2(*this, "Rebalancing prolongator -- import only", coarseLevel);
222  rebalancedP->doImport(*originalP, *trivialImporter, Xpetra::INSERT);
223  }
224  rebalancedP->fillComplete(importer->getTargetMap(), originalP->getRangeMap());
225 
226  } else {
227  rebalancedP = originalP;
228  RCP<const CrsMatrixWrap> crsOp = rcp_dynamic_cast<const CrsMatrixWrap>(originalP);
229  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
230 
231  RCP<CrsMatrix> rebalancedP2 = crsOp->getCrsMatrix();
232  TEUCHOS_TEST_FOR_EXCEPTION(rebalancedP2 == Teuchos::null, std::runtime_error, "Xpetra::CrsMatrixWrap doesn't have a CrsMatrix");
233 
234  {
235  SubFactoryMonitor subM(*this, "Rebalancing prolongator -- fast map replacement", coarseLevel);
236 
237  RCP<const Import> newImporter;
238  {
239  SubFactoryMonitor subM2(*this, "Import construction", coarseLevel);
240  newImporter = ImportFactory::Build(importer->getTargetMap(), rebalancedP->getColMap());
241  }
242  rebalancedP2->replaceDomainMapAndImporter(importer->getTargetMap(), newImporter);
243  }
244  }
246  // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
247  // That is probably something for an external permutation factory
248  // if (originalP->IsView("stridedMaps"))
249  // rebalancedP->CreateView("stridedMaps", originalP);
251  if (!rebalancedP.is_null()) {
252  std::ostringstream oss;
253  oss << "P_" << coarseLevel.GetLevelID();
254  rebalancedP->setObjectLabel(oss.str());
255  }
256  Set(coarseLevel, "P", rebalancedP);
257 
258  if (IsPrint(Statistics2))
259  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedP, "P (rebalanced)", params);
260  }
261  }
262 
263  if (importer.is_null()) {
264  if (IsAvailable(coarseLevel, "Nullspace"))
265  Set(coarseLevel, "Nullspace", Get<RCP<MultiVector> >(coarseLevel, "Nullspace"));
266 
267  if (pL.isParameter("Coordinates") && pL.get<RCP<const FactoryBase> >("Coordinates") != Teuchos::null)
268  if (IsAvailable(coarseLevel, "Coordinates"))
269  Set(coarseLevel, "Coordinates", Get<RCP<xdMV> >(coarseLevel, "Coordinates"));
270 
271  if (pL.isParameter("BlockNumber") && pL.get<RCP<const FactoryBase> >("BlockNumber") != Teuchos::null)
272  if (IsAvailable(coarseLevel, "BlockNumber"))
273  Set(coarseLevel, "BlockNumber", Get<RCP<LocalOrdinalVector> >(coarseLevel, "BlockNumber"));
274 
275  return;
276  }
277 
278  if (pL.isParameter("Coordinates") &&
279  pL.get<RCP<const FactoryBase> >("Coordinates") != Teuchos::null &&
280  IsAvailable(coarseLevel, "Coordinates")) {
281  RCP<xdMV> coords = Get<RCP<xdMV> >(coarseLevel, "Coordinates");
282 
283  // This line must be after the Get call
284  SubFactoryMonitor subM(*this, "Rebalancing coordinates", coarseLevel);
285 
286  LO nodeNumElts = coords->getMap()->getLocalNumElements();
287 
288  // If a process has no matrix rows, then we can't calculate blocksize using the formula below.
289  LO myBlkSize = 0, blkSize = 0;
290  if (nodeNumElts > 0)
291  myBlkSize = importer->getSourceMap()->getLocalNumElements() / nodeNumElts;
292  MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
293 
294  RCP<const Import> coordImporter;
295  if (blkSize == 1) {
296  coordImporter = importer;
297 
298  } else {
299  // NOTE: there is an implicit assumption here: we assume that dof any node are enumerated consequently
300  // Proper fix would require using decomposition similar to how we construct importer in the
301  // RepartitionFactory
302  RCP<const Map> origMap = coords->getMap();
303  GO indexBase = origMap->getIndexBase();
304 
305  ArrayView<const GO> OEntries = importer->getTargetMap()->getLocalElementList();
306  LO numEntries = OEntries.size() / blkSize;
307  ArrayRCP<GO> Entries(numEntries);
308  for (LO i = 0; i < numEntries; i++)
309  Entries[i] = (OEntries[i * blkSize] - indexBase) / blkSize + indexBase;
310 
311  RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
312  coordImporter = ImportFactory::Build(origMap, targetMap);
313  }
314 
315  RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
316  permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
317 
318  if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
319  permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
320 
321  if (permutedCoords->getMap() == Teuchos::null)
322  permutedCoords = Teuchos::null;
323 
324  Set(coarseLevel, "Coordinates", permutedCoords);
325 
326  std::string fileName = "rebalanced_coordinates_level_" + toString(coarseLevel.GetLevelID()) + ".m";
327  if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedCoords->getMap() != Teuchos::null)
328  Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO>::Write(fileName, *permutedCoords);
329  }
330 
331  if (pL.isParameter("BlockNumber") &&
332  pL.get<RCP<const FactoryBase> >("BlockNumber") != Teuchos::null &&
333  IsAvailable(coarseLevel, "BlockNumber")) {
334  RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(coarseLevel, "BlockNumber");
335 
336  // This line must be after the Get call
337  SubFactoryMonitor subM(*this, "Rebalancing BlockNumber", coarseLevel);
338 
339  RCP<LocalOrdinalVector> permutedBlockNumber = LocalOrdinalVectorFactory::Build(importer->getTargetMap(), false);
340  permutedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
341 
342  if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
343  permutedBlockNumber->replaceMap(permutedBlockNumber->getMap()->removeEmptyProcesses());
344 
345  if (permutedBlockNumber->getMap() == Teuchos::null)
346  permutedBlockNumber = Teuchos::null;
347 
348  Set(coarseLevel, "BlockNumber", permutedBlockNumber);
349 
350  std::string fileName = "rebalanced_BlockNumber_level_" + toString(coarseLevel.GetLevelID()) + ".m";
351  if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedBlockNumber->getMap() != Teuchos::null)
352  Xpetra::IO<LO, LO, GO, NO>::Write(fileName, *permutedBlockNumber);
353  }
354 
355  if (IsAvailable(coarseLevel, "Nullspace")) {
356  RCP<MultiVector> nullspace = Get<RCP<MultiVector> >(coarseLevel, "Nullspace");
357 
358  // This line must be after the Get call
359  SubFactoryMonitor subM(*this, "Rebalancing nullspace", coarseLevel);
360 
361  RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(importer->getTargetMap(), nullspace->getNumVectors());
362  permutedNullspace->doImport(*nullspace, *importer, Xpetra::INSERT);
363 
364  if (pL.get<bool>("repartition: use subcommunicators") == true)
365  permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
366 
367  if (permutedNullspace->getMap() == Teuchos::null)
368  permutedNullspace = Teuchos::null;
369 
370  Set(coarseLevel, "Nullspace", permutedNullspace);
371  }
372 
373  } else {
374  if (pL.get<bool>("transpose: use implicit") == false) {
375  RCP<Matrix> originalR = Get<RCP<Matrix> >(coarseLevel, "R");
376 
377  SubFactoryMonitor m2(*this, "Rebalancing restrictor", coarseLevel);
378 
379  if (implicit || importer.is_null()) {
380  GetOStream(Runtime0) << "Using original restrictor" << std::endl;
381  Set(coarseLevel, "R", originalR);
382 
383  } else {
384  RCP<Matrix> rebalancedR;
385  {
386  SubFactoryMonitor subM(*this, "Rebalancing restriction -- fusedImport", coarseLevel);
387 
388  RCP<Map> dummy; // meaning: use originalR's domain map.
389  Teuchos::ParameterList listLabel;
390  listLabel.set("Timer Label", "MueLu::RebalanceR-" + Teuchos::toString(coarseLevel.GetLevelID()));
391  rebalancedR = MatrixFactory::Build(originalR, *importer, dummy, importer->getTargetMap(), Teuchos::rcp(&listLabel, false));
392  }
393  if (!rebalancedR.is_null()) {
394  std::ostringstream oss;
395  oss << "R_" << coarseLevel.GetLevelID();
396  rebalancedR->setObjectLabel(oss.str());
397  }
398  Set(coarseLevel, "R", rebalancedR);
399 
401  // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
402  // That is probably something for an external permutation factory
403  // if (originalR->IsView("stridedMaps"))
404  // rebalancedR->CreateView("stridedMaps", originalR);
406 
407  if (IsPrint(Statistics2))
408  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedR, "R (rebalanced)", params);
409  }
410  }
411  }
412 }
413 
414 } // namespace MueLu
415 
416 #endif // MUELU_REBALANCETRANSFERFACTORY_DEF_HPP
Exception indicating invalid cast attempted.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
#define MueLu_maxAll(rcpComm, in, out)
GlobalOrdinal GO
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 void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
size_type size() const
LocalOrdinal LO
One-liner description of what is happening.
static const NoFactory * get()
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Print even more statistics.
bool isParameter(const std::string &name) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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 Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
#define SET_VALID_ENTRY(name)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Node NO
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
std::string toString(const T &t)
bool is_null() const