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