MueLu  Version of the Day
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: rebalance Nullspace");
79  SET_VALID_ENTRY("transpose: use implicit");
80  SET_VALID_ENTRY("repartition: use subcommunicators");
81 #undef SET_VALID_ENTRY
82 
83  {
84  typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
85  RCP<validatorType> typeValidator = rcp (new validatorType(Teuchos::tuple<std::string>("Interpolation", "Restriction"), "type"));
86  validParamList->set("type", "Interpolation", "Type of the transfer operator that need to be rebalanced (Interpolation or Restriction)", typeValidator);
87  }
88 
89  validParamList->set< RCP<const FactoryBase> >("P", null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
90  validParamList->set< RCP<const FactoryBase> >("R", null, "Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
91  validParamList->set< RCP<const FactoryBase> >("Nullspace", null, "Factory of the nullspace that need to be rebalanced (only used if type=Interpolation)");
92  validParamList->set< RCP<const FactoryBase> >("Coordinates", null, "Factory of the coordinates that need to be rebalanced (only used if type=Interpolation)");
93  validParamList->set< RCP<const FactoryBase> >("BlockNumber", null, "Factory of the block ids that need to be rebalanced (only used if type=Interpolation)");
94  validParamList->set< RCP<const FactoryBase> >("Importer", null, "Factory of the importer object used for the rebalancing");
95  validParamList->set< int > ("write start", -1, "First level at which coordinates should be written to file");
96  validParamList->set< int > ("write end", -1, "Last level at which coordinates should be written to file");
97 
98  // TODO validation: "P" parameter valid only for type="Interpolation" and "R" valid only for type="Restriction". Like so:
99  // if (paramList.isEntry("type") && paramList.get("type) == "Interpolation) {
100  // validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
101 
102  return validParamList;
103  }
104 
105  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
107  const ParameterList& pL = GetParameterList();
108 
109  if (pL.get<std::string>("type") == "Interpolation") {
110  Input(coarseLevel, "P");
111  if (pL.get<bool>("repartition: rebalance Nullspace"))
112  Input(coarseLevel, "Nullspace");
113  if (pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null)
114  Input(coarseLevel, "Coordinates");
115  if (pL.get< RCP<const FactoryBase> >("BlockNumber") != Teuchos::null)
116  Input(coarseLevel, "BlockNumber");
117 
118  } else {
119  if (pL.get<bool>("transpose: use implicit") == false)
120  Input(coarseLevel, "R");
121  }
122 
123  Input(coarseLevel, "Importer");
124  }
125 
126  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
128  FactoryMonitor m(*this, "Build", coarseLevel);
129  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO> xdMV;
130 
131  const ParameterList& pL = GetParameterList();
132 
133  int implicit = !pL.get<bool>("repartition: rebalance P and R");
134  int writeStart = pL.get<int> ("write start");
135  int writeEnd = pL.get<int> ("write end");
136 
137  if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel, "Coordinates")) {
138  std::string fileName = "coordinates_level_0.m";
139  RCP<xdMV> fineCoords = fineLevel.Get< RCP<xdMV> >("Coordinates");
140  if (fineCoords != Teuchos::null)
141  Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *fineCoords);
142  }
143 
144  if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel, "BlockNumber")) {
145  std::string fileName = "BlockNumber_level_0.m";
146  RCP<LocalOrdinalVector> fineBlockNumber = fineLevel.Get< RCP<LocalOrdinalVector> >("BlockNumber");
147  if (fineBlockNumber != Teuchos::null)
148  Xpetra::IO<LO,LO,GO,NO>::Write(fileName, *fineBlockNumber);
149  }
150 
151  RCP<const Import> importer = Get<RCP<const Import> >(coarseLevel, "Importer");
152  if (implicit) {
153  // Save the importer, we'll need it for solve
154  coarseLevel.Set("Importer", importer, NoFactory::get());
155  }
156 
157  RCP<ParameterList> params = rcp(new ParameterList());
158  if (IsPrint(Statistics2)) {
159  params->set("printLoadBalancingInfo", true);
160  params->set("printCommInfo", true);
161  }
162 
163  std::string transferType = pL.get<std::string>("type");
164  if (transferType == "Interpolation") {
165  RCP<Matrix> originalP = Get< RCP<Matrix> >(coarseLevel, "P");
166 
167  {
168  // This line must be after the Get call
169  SubFactoryMonitor m1(*this, "Rebalancing prolongator", coarseLevel);
170 
171  if (implicit || importer.is_null()) {
172  GetOStream(Runtime0) << "Using original prolongator" << std::endl;
173  Set(coarseLevel, "P", originalP);
174 
175  } else {
176  // P is the transfer operator from the coarse grid to the fine grid.
177  // P must transfer the data from the newly reordered coarse A to the
178  // (unchanged) fine A. This means that the domain map (coarse) of P
179  // must be changed according to the new partition. The range map
180  // (fine) is kept unchanged.
181  //
182  // The domain map of P must match the range map of R. See also note
183  // below about domain/range map of R and its implications for P.
184  //
185  // To change the domain map of P, P needs to be fillCompleted again
186  // with the new domain map. To achieve this, P is copied into a new
187  // matrix that is not fill-completed. The doImport() operation is
188  // just used here to make a copy of P: the importer is trivial and
189  // there is no data movement involved. The reordering actually
190  // happens during the fillComplete() with domainMap == importer->getTargetMap().
191  RCP<Matrix> rebalancedP = originalP;
192  RCP<const CrsMatrixWrap> crsOp = rcp_dynamic_cast<const CrsMatrixWrap>(originalP);
193  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
194 
195  RCP<CrsMatrix> rebalancedP2 = crsOp->getCrsMatrix();
196  TEUCHOS_TEST_FOR_EXCEPTION(rebalancedP2 == Teuchos::null, std::runtime_error, "Xpetra::CrsMatrixWrap doesn't have a CrsMatrix");
197 
198  {
199  SubFactoryMonitor subM(*this, "Rebalancing prolongator -- fast map replacement", coarseLevel);
200 
201  RCP<const Import> newImporter;
202  {
203  SubFactoryMonitor(*this, "Import construction", coarseLevel);
204  newImporter = ImportFactory::Build(importer->getTargetMap(), rebalancedP->getColMap());
205  }
206  rebalancedP2->replaceDomainMapAndImporter(importer->getTargetMap(), newImporter);
207  }
208 
210  // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
211  // That is probably something for an external permutation factory
212  // if (originalP->IsView("stridedMaps"))
213  // rebalancedP->CreateView("stridedMaps", originalP);
215  if(!rebalancedP.is_null()) {std::ostringstream oss; oss << "P_" << coarseLevel.GetLevelID(); rebalancedP->setObjectLabel(oss.str());}
216  Set(coarseLevel, "P", rebalancedP);
217 
218  if (IsPrint(Statistics2))
219  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedP, "P (rebalanced)", params);
220  }
221  }
222 
223  if (importer.is_null()) {
224  if (IsAvailable(coarseLevel, "Nullspace"))
225  Set(coarseLevel, "Nullspace", Get<RCP<MultiVector> >(coarseLevel, "Nullspace"));
226 
227  if (pL.isParameter("Coordinates") && pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null)
228  if (IsAvailable(coarseLevel, "Coordinates"))
229  Set(coarseLevel, "Coordinates", Get< RCP<xdMV> >(coarseLevel, "Coordinates"));
230 
231  if (pL.isParameter("BlockNumber") && pL.get< RCP<const FactoryBase> >("BlockNumber") != Teuchos::null)
232  if (IsAvailable(coarseLevel, "BlockNumber"))
233  Set(coarseLevel, "BlockNumber", Get< RCP<LocalOrdinalVector> >(coarseLevel, "BlockNumber"));
234 
235  return;
236  }
237 
238  if (pL.isParameter("Coordinates") &&
239  pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null &&
240  IsAvailable(coarseLevel, "Coordinates")) {
241  RCP<xdMV> coords = Get<RCP<xdMV> >(coarseLevel, "Coordinates");
242 
243  // This line must be after the Get call
244  SubFactoryMonitor subM(*this, "Rebalancing coordinates", coarseLevel);
245 
246  LO nodeNumElts = coords->getMap()->getNodeNumElements();
247 
248  // If a process has no matrix rows, then we can't calculate blocksize using the formula below.
249  LO myBlkSize = 0, blkSize = 0;
250  if (nodeNumElts > 0)
251  myBlkSize = importer->getSourceMap()->getNodeNumElements() / nodeNumElts;
252  MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
253 
254  RCP<const Import> coordImporter;
255  if (blkSize == 1) {
256  coordImporter = importer;
257 
258  } else {
259  // NOTE: there is an implicit assumption here: we assume that dof any node are enumerated consequently
260  // Proper fix would require using decomposition similar to how we construct importer in the
261  // RepartitionFactory
262  RCP<const Map> origMap = coords->getMap();
263  GO indexBase = origMap->getIndexBase();
264 
265  ArrayView<const GO> OEntries = importer->getTargetMap()->getNodeElementList();
266  LO numEntries = OEntries.size()/blkSize;
267  ArrayRCP<GO> Entries(numEntries);
268  for (LO i = 0; i < numEntries; i++)
269  Entries[i] = (OEntries[i*blkSize]-indexBase)/blkSize + indexBase;
270 
271  RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
272  coordImporter = ImportFactory::Build(origMap, targetMap);
273  }
274 
275  RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
276  permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
277 
278  if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
279  permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
280 
281  if (permutedCoords->getMap() == Teuchos::null)
282  permutedCoords = Teuchos::null;
283 
284  Set(coarseLevel, "Coordinates", permutedCoords);
285 
286  std::string fileName = "rebalanced_coordinates_level_" + toString(coarseLevel.GetLevelID()) + ".m";
287  if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedCoords->getMap() != Teuchos::null)
288  Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *permutedCoords);
289  }
290 
291  if (pL.isParameter("BlockNumber") &&
292  pL.get< RCP<const FactoryBase> >("BlockNumber") != Teuchos::null &&
293  IsAvailable(coarseLevel, "BlockNumber")) {
294  RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(coarseLevel, "BlockNumber");
295 
296  // This line must be after the Get call
297  SubFactoryMonitor subM(*this, "Rebalancing BlockNumber", coarseLevel);
298 
299  RCP<LocalOrdinalVector> permutedBlockNumber = LocalOrdinalVectorFactory::Build(importer->getTargetMap(), false);
300  permutedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
301 
302  if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
303  permutedBlockNumber->replaceMap(permutedBlockNumber->getMap()->removeEmptyProcesses());
304 
305  if (permutedBlockNumber->getMap() == Teuchos::null)
306  permutedBlockNumber = Teuchos::null;
307 
308  Set(coarseLevel, "BlockNumber", permutedBlockNumber);
309 
310  std::string fileName = "rebalanced_BlockNumber_level_" + toString(coarseLevel.GetLevelID()) + ".m";
311  if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedBlockNumber->getMap() != Teuchos::null)
312  Xpetra::IO<LO,LO,GO,NO>::Write(fileName, *permutedBlockNumber);
313  }
314 
315  if (IsAvailable(coarseLevel, "Nullspace")) {
316  RCP<MultiVector> nullspace = Get< RCP<MultiVector> >(coarseLevel, "Nullspace");
317 
318  // This line must be after the Get call
319  SubFactoryMonitor subM(*this, "Rebalancing nullspace", coarseLevel);
320 
321  RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(importer->getTargetMap(), nullspace->getNumVectors());
322  permutedNullspace->doImport(*nullspace, *importer, Xpetra::INSERT);
323 
324  if (pL.get<bool>("repartition: use subcommunicators") == true)
325  permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
326 
327  if (permutedNullspace->getMap() == Teuchos::null)
328  permutedNullspace = Teuchos::null;
329 
330  Set(coarseLevel, "Nullspace", permutedNullspace);
331  }
332 
333  } else {
334  if (pL.get<bool>("transpose: use implicit") == false) {
335  RCP<Matrix> originalR = Get< RCP<Matrix> >(coarseLevel, "R");
336 
337  SubFactoryMonitor m2(*this, "Rebalancing restrictor", coarseLevel);
338 
339  if (implicit || importer.is_null()) {
340  GetOStream(Runtime0) << "Using original restrictor" << std::endl;
341  Set(coarseLevel, "R", originalR);
342 
343  } else {
344  RCP<Matrix> rebalancedR;
345  {
346  SubFactoryMonitor subM(*this, "Rebalancing restriction -- fusedImport", coarseLevel);
347 
348  RCP<Map> dummy; // meaning: use originalR's domain map.
349  Teuchos::ParameterList listLabel;
350  listLabel.set("Timer Label","MueLu::RebalanceR-" + Teuchos::toString(coarseLevel.GetLevelID()));
351  rebalancedR = MatrixFactory::Build(originalR, *importer, dummy, importer->getTargetMap(),Teuchos::rcp(&listLabel,false));
352  }
353  if(!rebalancedR.is_null()) {std::ostringstream oss; oss << "R_" << coarseLevel.GetLevelID(); rebalancedR->setObjectLabel(oss.str());}
354  Set(coarseLevel, "R", rebalancedR);
355 
357  // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
358  // That is probably something for an external permutation factory
359  // if (originalR->IsView("stridedMaps"))
360  // rebalancedR->CreateView("stridedMaps", originalR);
362 
363  if (IsPrint(Statistics2))
364  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedR, "R (rebalanced)", params);
365  }
366  }
367  }
368  }
369 
370 } // namespace MueLu
371 
372 #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->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
#define MueLu_maxAll(rcpComm, in, out)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
static const NoFactory * get()
Print even more statistics.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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.
#define SET_VALID_ENTRY(name)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)