MueLu  Version of the Day
MueLu_RAPFactory_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_RAPFACTORY_DEF_HPP
47 #define MUELU_RAPFACTORY_DEF_HPP
48 
49 
50 #include <sstream>
51 
52 #include <Xpetra_Matrix.hpp>
53 #include <Xpetra_MatrixFactory.hpp>
54 #include <Xpetra_MatrixMatrix.hpp>
55 #include <Xpetra_MatrixUtils.hpp>
56 #include <Xpetra_TripleMatrixMultiply.hpp>
57 #include <Xpetra_Vector.hpp>
58 #include <Xpetra_VectorFactory.hpp>
59 
61 
62 #include "MueLu_MasterList.hpp"
63 #include "MueLu_Monitor.hpp"
64 #include "MueLu_PerfUtils.hpp"
66 //#include "MueLu_Utilities.hpp"
67 
68 namespace MueLu {
69 
70  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
72  : hasDeclaredInput_(false) { }
73 
74  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
76  RCP<ParameterList> validParamList = rcp(new ParameterList());
77 
78 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
79  SET_VALID_ENTRY("transpose: use implicit");
80  SET_VALID_ENTRY("rap: triple product");
81  SET_VALID_ENTRY("rap: fix zero diagonals");
82  SET_VALID_ENTRY("rap: fix zero diagonals threshold");
83  SET_VALID_ENTRY("rap: fix zero diagonals replacement");
84  SET_VALID_ENTRY("rap: relative diagonal floor");
85 #undef SET_VALID_ENTRY
86  validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A used during the prolongator smoothing process");
87  validParamList->set< RCP<const FactoryBase> >("P", null, "Prolongator factory");
88  validParamList->set< RCP<const FactoryBase> >("R", null, "Restrictor factory");
89 
90  validParamList->set< bool > ("CheckMainDiagonal", false, "Check main diagonal for zeros");
91  validParamList->set< bool > ("RepairMainDiagonal", false, "Repair zeros on main diagonal");
92 
93  // Make sure we don't recursively validate options for the matrixmatrix kernels
94  ParameterList norecurse;
95  norecurse.disableRecursiveValidation();
96  validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
97 
98  return validParamList;
99  }
100 
101  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
103  const Teuchos::ParameterList& pL = GetParameterList();
104  if (pL.get<bool>("transpose: use implicit") == false)
105  Input(coarseLevel, "R");
106 
107  Input(fineLevel, "A");
108  Input(coarseLevel, "P");
109 
110  // call DeclareInput of all user-given transfer factories
111  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
112  (*it)->CallDeclareInput(coarseLevel);
113 
114  hasDeclaredInput_ = true;
115  }
116 
117  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
119  const bool doTranspose = true;
120  const bool doFillComplete = true;
121  const bool doOptimizeStorage = true;
122  RCP<Matrix> Ac;
123  {
124  FactoryMonitor m(*this, "Computing Ac", coarseLevel);
125  std::ostringstream levelstr;
126  levelstr << coarseLevel.GetLevelID();
127  std::string labelstr = FormattingHelper::getColonLabel(coarseLevel.getObjectLabel());
128 
129  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_ == false, Exceptions::RuntimeError,
130  "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
131 
132  const Teuchos::ParameterList& pL = GetParameterList();
133  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
134  RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P"), AP;
135 
136  bool isEpetra = A->getRowMap()->lib() == Xpetra::UseEpetra;
137  bool isGPU =
138 #ifdef KOKKOS_ENABLE_CUDA
139  (typeid(Node).name() == typeid(Kokkos::Compat::KokkosCudaWrapperNode).name()) ||
140 #endif
141 #ifdef KOKKOS_ENABLE_HIP
142  (typeid(Node).name() == typeid(Kokkos::Compat::KokkosHIPWrapperNode).name()) ||
143 #endif
144  false;
145 
146  if (pL.get<bool>("rap: triple product") == false || isEpetra || isGPU) {
147  if (pL.get<bool>("rap: triple product") && isEpetra)
148  GetOStream(Warnings1) << "Switching from triple product to R x (A x P) since triple product has not been implemented for Epetra.\n";
149 #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
150  if (pL.get<bool>("rap: triple product") && isGPU)
151  GetOStream(Warnings1) << "Switching from triple product to R x (A x P) since triple product has not been implemented for "
152  << Node::execution_space::name() << std::endl;
153 #endif
154 
155  // Reuse pattern if available (multiple solve)
156  RCP<ParameterList> APparams = rcp(new ParameterList);
157  if(pL.isSublist("matrixmatrix: kernel params"))
158  APparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
159 
160  // By default, we don't need global constants for A*P
161  APparams->set("compute global constants: temporaries",APparams->get("compute global constants: temporaries",false));
162  APparams->set("compute global constants",APparams->get("compute global constants",false));
163 
164  if (coarseLevel.IsAvailable("AP reuse data", this)) {
165  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
166 
167  APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
168 
169  if (APparams->isParameter("graph"))
170  AP = APparams->get< RCP<Matrix> >("graph");
171  }
172 
173  {
174  SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
175 
176  AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(Statistics2),
177  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::A*P-")+levelstr.str(), APparams);
178  }
179 
180  // Reuse coarse matrix memory if available (multiple solve)
181  RCP<ParameterList> RAPparams = rcp(new ParameterList);
182  if(pL.isSublist("matrixmatrix: kernel params"))
183  RAPparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
184 
185  if (coarseLevel.IsAvailable("RAP reuse data", this)) {
186  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous RAP data" << std::endl;
187 
188  RAPparams = coarseLevel.Get< RCP<ParameterList> >("RAP reuse data", this);
189 
190  if (RAPparams->isParameter("graph"))
191  Ac = RAPparams->get< RCP<Matrix> >("graph");
192 
193  // Some eigenvalue may have been cached with the matrix in the previous run.
194  // As the matrix values will be updated, we need to reset the eigenvalue.
195  Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
196  }
197 
198  // We *always* need global constants for the RAP, but not for the temps
199  RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
200  RAPparams->set("compute global constants",true);
201 
202  // Allow optimization of storage.
203  // This is necessary for new faster Epetra MM kernels.
204  // Seems to work with matrix modifications to repair diagonal entries.
205 
206  if (pL.get<bool>("transpose: use implicit") == true) {
207  SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
208 
209  Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
210  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
211 
212  } else {
213  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
214 
215  SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
216 
217  Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
218  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
219  }
220 
221  Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >("rap: relative diagonal floor")();
222  if(relativeFloor.size() > 0) {
223  Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(Statistics2));
224  }
225 
226  bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
227  bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
228  if (checkAc || repairZeroDiagonals) {
229  using magnitudeType = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
230  magnitudeType threshold;
231  if (pL.isType<magnitudeType>("rap: fix zero diagonals threshold"))
232  threshold = pL.get<magnitudeType>("rap: fix zero diagonals threshold");
233  else
234  threshold = Teuchos::as<magnitudeType>(pL.get<double>("rap: fix zero diagonals threshold"));
235  Scalar replacement = Teuchos::as<Scalar>(pL.get<double>("rap: fix zero diagonals replacement"));
236  Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1), threshold, replacement);
237  }
238 
239  if (IsPrint(Statistics2)) {
240  RCP<ParameterList> params = rcp(new ParameterList());;
241  params->set("printLoadBalancingInfo", true);
242  params->set("printCommInfo", true);
243  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
244  }
245 
246  if(!Ac.is_null()) {std::ostringstream oss; oss << "A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
247  Set(coarseLevel, "A", Ac);
248 
249  APparams->set("graph", AP);
250  Set(coarseLevel, "AP reuse data", APparams);
251  RAPparams->set("graph", Ac);
252  Set(coarseLevel, "RAP reuse data", RAPparams);
253  } else {
254  RCP<ParameterList> RAPparams = rcp(new ParameterList);
255  if(pL.isSublist("matrixmatrix: kernel params"))
256  RAPparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
257 
258  if (coarseLevel.IsAvailable("RAP reuse data", this)) {
259  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous RAP data" << std::endl;
260 
261  RAPparams = coarseLevel.Get< RCP<ParameterList> >("RAP reuse data", this);
262 
263  if (RAPparams->isParameter("graph"))
264  Ac = RAPparams->get< RCP<Matrix> >("graph");
265 
266  // Some eigenvalue may have been cached with the matrix in the previous run.
267  // As the matrix values will be updated, we need to reset the eigenvalue.
268  Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
269  }
270 
271  // We *always* need global constants for the RAP, but not for the temps
272  RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
273  RAPparams->set("compute global constants",true);
274 
275  if (pL.get<bool>("transpose: use implicit") == true) {
276 
277  Ac = MatrixFactory::Build(P->getDomainMap(), Teuchos::as<LO>(0));
278 
279  SubFactoryMonitor m2(*this, "MxMxM: R x A x P (implicit)", coarseLevel);
280 
281  Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
282  MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
283  doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-implicit-")+levelstr.str(),
284  RAPparams);
285 
286  } else {
287  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
288  Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));
289 
290  SubFactoryMonitor m2(*this, "MxMxM: R x A x P (explicit)", coarseLevel);
291 
292  Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
293  MultiplyRAP(*R, !doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
294  doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-explicit-")+levelstr.str(),
295  RAPparams);
296  }
297 
298  Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >("rap: relative diagonal floor")();
299  if(relativeFloor.size() > 0) {
300  Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(Statistics2));
301  }
302 
303  bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
304  bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
305  if (checkAc || repairZeroDiagonals) {
306  using magnitudeType = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
307  magnitudeType threshold;
308  if (pL.isType<magnitudeType>("rap: fix zero diagonals threshold"))
309  threshold = pL.get<magnitudeType>("rap: fix zero diagonals threshold");
310  else
311  threshold = Teuchos::as<magnitudeType>(pL.get<double>("rap: fix zero diagonals threshold"));
312  Scalar replacement = Teuchos::as<Scalar>(pL.get<double>("rap: fix zero diagonals replacement"));
313  Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1), threshold, replacement);
314  }
315 
316 
317  if (IsPrint(Statistics2)) {
318  RCP<ParameterList> params = rcp(new ParameterList());;
319  params->set("printLoadBalancingInfo", true);
320  params->set("printCommInfo", true);
321  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
322  }
323 
324  if(!Ac.is_null()) {std::ostringstream oss; oss << "A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
325  Set(coarseLevel, "A", Ac);
326 
327  RAPparams->set("graph", Ac);
328  Set(coarseLevel, "RAP reuse data", RAPparams);
329  }
330 
331 
332  }
333 
334 #ifdef HAVE_MUELU_DEBUG
335  MatrixUtils::checkLocalRowMapMatchesColMap(*Ac);
336 #endif // HAVE_MUELU_DEBUG
337 
338  if (transferFacts_.begin() != transferFacts_.end()) {
339  SubFactoryMonitor m(*this, "Projections", coarseLevel);
340 
341  // call Build of all user-given transfer factories
342  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
343  RCP<const FactoryBase> fac = *it;
344  GetOStream(Runtime0) << "RAPFactory: call transfer factory: " << fac->description() << std::endl;
345  fac->CallBuild(coarseLevel);
346  // Coordinates transfer is marginally different from all other operations
347  // because it is *optional*, and not required. For instance, we may need
348  // coordinates only on level 4 if we start repartitioning from that level,
349  // but we don't need them on level 1,2,3. As our current Hierarchy setup
350  // assumes propagation of dependencies only through three levels, this
351  // means that we need to rely on other methods to propagate optional data.
352  //
353  // The method currently used is through RAP transfer factories, which are
354  // simply factories which are called at the end of RAP with a single goal:
355  // transfer some fine data to coarser level. Because these factories are
356  // kind of outside of the mainline factories, they behave different. In
357  // particular, we call their Build method explicitly, rather than through
358  // Get calls. This difference is significant, as the Get call is smart
359  // enough to know when to release all factory dependencies, and Build is
360  // dumb. This led to the following CoordinatesTransferFactory sequence:
361  // 1. Request level 0
362  // 2. Request level 1
363  // 3. Request level 0
364  // 4. Release level 0
365  // 5. Release level 1
366  //
367  // The problem is missing "6. Release level 0". Because it was missing,
368  // we had outstanding request on "Coordinates", "Aggregates" and
369  // "CoarseMap" on level 0.
370  //
371  // This was fixed by explicitly calling Release on transfer factories in
372  // RAPFactory. I am still unsure how exactly it works, but now we have
373  // clear data requests for all levels.
374  coarseLevel.Release(*fac);
375  }
376  }
377 
378  }
379 
380  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
382  // check if it's a TwoLevelFactoryBase based transfer factory
383  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast,
384  "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. "
385  "This is very strange. (Note: you can remove this exception if there's a good reason for)");
386  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_, Exceptions::RuntimeError, "MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
387  transferFacts_.push_back(factory);
388  }
389 
390 } //namespace MueLu
391 
392 #define MUELU_RAPFACTORY_SHORT
393 #endif // MUELU_RAPFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
Namespace for MueLu classes and methods.
MueLu::DefaultNode Node
static std::string getColonLabel(const std::string &label)
Helper function for object label.
Print even more statistics.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Additional warnings.
MueLu::DefaultScalar Scalar
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.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.