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