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_MultiVector.hpp"
52 #include "Xpetra_MultiVectorFactory.hpp"
53 #include "Xpetra_Vector.hpp"
54 #include "Xpetra_VectorFactory.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  {
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=Interpolation)");
90  validParamList->set< RCP<const FactoryBase> >("Coordinates", null, "Factory of the coordinates that need to be rebalanced (only used if type=Interpolation)");
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  if (IsPrint(Statistics2)) {
146  params->set("printLoadBalancingInfo", true);
147  params->set("printCommInfo", true);
148  }
149 
150  std::string transferType = pL.get<std::string>("type");
151  if (transferType == "Interpolation") {
152  RCP<Matrix> originalP = Get< RCP<Matrix> >(coarseLevel, "P");
153 
154  {
155  // This line must be after the Get call
156  SubFactoryMonitor m1(*this, "Rebalancing prolongator", coarseLevel);
157 
158  if (implicit || importer.is_null()) {
159  GetOStream(Runtime0) << "Using original prolongator" << std::endl;
160  Set(coarseLevel, "P", originalP);
161 
162  } else {
163  // P is the transfer operator from the coarse grid to the fine grid.
164  // P must transfer the data from the newly reordered coarse A to the
165  // (unchanged) fine A. This means that the domain map (coarse) of P
166  // must be changed according to the new partition. The range map
167  // (fine) is kept unchanged.
168  //
169  // The domain map of P must match the range map of R. See also note
170  // below about domain/range map of R and its implications for P.
171  //
172  // To change the domain map of P, P needs to be fillCompleted again
173  // with the new domain map. To achieve this, P is copied into a new
174  // matrix that is not fill-completed. The doImport() operation is
175  // just used here to make a copy of P: the importer is trivial and
176  // there is no data movement involved. The reordering actually
177  // happens during the fillComplete() with domainMap == importer->getTargetMap().
178  RCP<Matrix> rebalancedP = originalP;
179  RCP<const CrsMatrixWrap> crsOp = rcp_dynamic_cast<const CrsMatrixWrap>(originalP);
180  TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
181 
182  RCP<CrsMatrix> rebalancedP2 = crsOp->getCrsMatrix();
183  TEUCHOS_TEST_FOR_EXCEPTION(rebalancedP2 == Teuchos::null, std::runtime_error, "Xpetra::CrsMatrixWrap doesn't have a CrsMatrix");
184 
185  {
186  SubFactoryMonitor subM(*this, "Rebalancing prolongator -- fast map replacement", coarseLevel);
187 
188  RCP<const Import> newImporter;
189  {
190  SubFactoryMonitor(*this, "Import construction", coarseLevel);
191  newImporter = ImportFactory::Build(importer->getTargetMap(), rebalancedP->getColMap());
192  }
193  rebalancedP2->replaceDomainMapAndImporter(importer->getTargetMap(), newImporter);
194  }
195 
197  // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
198  // That is probably something for an external permutation factory
199  // if (originalP->IsView("stridedMaps"))
200  // rebalancedP->CreateView("stridedMaps", originalP);
202 
203  Set(coarseLevel, "P", rebalancedP);
204 
205  if (IsPrint(Statistics2))
206  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedP, "P (rebalanced)", params);
207  }
208  }
209 
210  if (importer.is_null()) {
211  if (IsAvailable(coarseLevel, "Nullspace"))
212  Set(coarseLevel, "Nullspace", Get<RCP<MultiVector> >(coarseLevel, "Nullspace"));
213 
214  if (pL.isParameter("Coordinates") && pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null)
215  if (IsAvailable(coarseLevel, "Coordinates"))
216  Set(coarseLevel, "Coordinates", Get< RCP<xdMV> >(coarseLevel, "Coordinates"));
217 
218  return;
219  }
220 
221  if (pL.isParameter("Coordinates") &&
222  pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null &&
223  IsAvailable(coarseLevel, "Coordinates")) {
224  RCP<xdMV> coords = Get<RCP<xdMV> >(coarseLevel, "Coordinates");
225 
226  // This line must be after the Get call
227  SubFactoryMonitor subM(*this, "Rebalancing coordinates", coarseLevel);
228 
229  LO nodeNumElts = coords->getMap()->getNodeNumElements();
230 
231  // If a process has no matrix rows, then we can't calculate blocksize using the formula below.
232  LO myBlkSize = 0, blkSize = 0;
233  if (nodeNumElts > 0)
234  myBlkSize = importer->getSourceMap()->getNodeNumElements() / nodeNumElts;
235  MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
236 
237  RCP<const Import> coordImporter;
238  if (blkSize == 1) {
239  coordImporter = importer;
240 
241  } else {
242  // NOTE: there is an implicit assumption here: we assume that dof any node are enumerated consequently
243  // Proper fix would require using decomposition similar to how we construct importer in the
244  // RepartitionFactory
245  RCP<const Map> origMap = coords->getMap();
246  GO indexBase = origMap->getIndexBase();
247 
248  ArrayView<const GO> OEntries = importer->getTargetMap()->getNodeElementList();
249  LO numEntries = OEntries.size()/blkSize;
250  ArrayRCP<GO> Entries(numEntries);
251  for (LO i = 0; i < numEntries; i++)
252  Entries[i] = (OEntries[i*blkSize]-indexBase)/blkSize + indexBase;
253 
254  RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
255  coordImporter = ImportFactory::Build(origMap, targetMap);
256  }
257 
258  RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
259  permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
260 
261  if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
262  permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
263 
264  Set(coarseLevel, "Coordinates", permutedCoords);
265 
266  std::string fileName = "rebalanced_coordinates_level_" + toString(coarseLevel.GetLevelID()) + ".m";
267  if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedCoords->getMap() != Teuchos::null)
268  Xpetra::IO<double,LO,GO,NO>::Write(fileName, *permutedCoords);
269  }
270 
271  if (IsAvailable(coarseLevel, "Nullspace")) {
272  RCP<MultiVector> nullspace = Get< RCP<MultiVector> >(coarseLevel, "Nullspace");
273 
274  // This line must be after the Get call
275  SubFactoryMonitor subM(*this, "Rebalancing nullspace", coarseLevel);
276 
277  RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(importer->getTargetMap(), nullspace->getNumVectors());
278  permutedNullspace->doImport(*nullspace, *importer, Xpetra::INSERT);
279 
280  if (pL.get<bool>("repartition: use subcommunicators") == true)
281  permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
282 
283  Set(coarseLevel, "Nullspace", permutedNullspace);
284  }
285 
286  } else {
287  if (pL.get<bool>("transpose: use implicit") == false) {
288  RCP<Matrix> originalR = Get< RCP<Matrix> >(coarseLevel, "R");
289 
290  SubFactoryMonitor m2(*this, "Rebalancing restrictor", coarseLevel);
291 
292  if (implicit || importer.is_null()) {
293  GetOStream(Runtime0) << "Using original restrictor" << std::endl;
294  Set(coarseLevel, "R", originalR);
295 
296  } else {
297  RCP<Matrix> rebalancedR;
298  {
299  SubFactoryMonitor subM(*this, "Rebalancing restriction -- fusedImport", coarseLevel);
300 
301  RCP<Map> dummy; // meaning: use originalR's domain map.
302  Teuchos::ParameterList listLabel;
303  listLabel.set("Timer Label","MueLu::RebalanceR-" + Teuchos::toString(coarseLevel.GetLevelID()));
304  rebalancedR = MatrixFactory::Build(originalR, *importer, dummy, importer->getTargetMap(),Teuchos::rcp(&listLabel,false));
305  }
306  Set(coarseLevel, "R", rebalancedR);
307 
309  // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
310  // That is probably something for an external permutation factory
311  // if (originalR->IsView("stridedMaps"))
312  // rebalancedR->CreateView("stridedMaps", originalR);
314 
315  if (IsPrint(Statistics2))
316  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedR, "R (rebalanced)", params);
317  }
318  }
319  }
320  }
321 
322 } // namespace MueLu
323 
324 #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.
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)
One-liner description of what is happening.
size_type size() const
Namespace for MueLu classes and methods.
static const NoFactory * get()
bool is_null() const
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.
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.
#define SET_VALID_ENTRY(name)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
bool isParameter(const std::string &name) const
std::string toString(const T &t)