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 #undef SET_VALID_ENTRY
83  validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A used during the prolongator smoothing process");
84  validParamList->set< RCP<const FactoryBase> >("P", null, "Prolongator factory");
85  validParamList->set< RCP<const FactoryBase> >("R", null, "Restrictor factory");
86 
87  validParamList->set< bool > ("CheckMainDiagonal", false, "Check main diagonal for zeros");
88  validParamList->set< bool > ("RepairMainDiagonal", false, "Repair zeros on main diagonal");
89 
90  // Make sure we don't recursively validate options for the matrixmatrix kernels
91  ParameterList norecurse;
92  norecurse.disableRecursiveValidation();
93  validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
94 
95  return validParamList;
96  }
97 
98  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100  const Teuchos::ParameterList& pL = GetParameterList();
101  if (pL.get<bool>("transpose: use implicit") == false)
102  Input(coarseLevel, "R");
103 
104  Input(fineLevel, "A");
105  Input(coarseLevel, "P");
106 
107  // call DeclareInput of all user-given transfer factories
108  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
109  (*it)->CallDeclareInput(coarseLevel);
110 
111  hasDeclaredInput_ = true;
112  }
113 
114  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
116  const bool doTranspose = true;
117  const bool doFillComplete = true;
118  const bool doOptimizeStorage = true;
119  {
120  FactoryMonitor m(*this, "Computing Ac", coarseLevel);
121  std::ostringstream levelstr;
122  levelstr << coarseLevel.GetLevelID();
123  std::string labelstr = FormattingHelper::getColonLabel(coarseLevel.getObjectLabel());
124 
125  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_ == false, Exceptions::RuntimeError,
126  "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
127 
128  const Teuchos::ParameterList& pL = GetParameterList();
129  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
130  RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P"), AP, Ac;
131 
132 
133  if (pL.get<bool>("rap: triple product") == false) {
134  // Reuse pattern if available (multiple solve)
135  RCP<ParameterList> APparams;
136  if(pL.isSublist("matrixmatrix: kernel params"))
137  APparams=rcp(new ParameterList(pL.sublist("matrixmatrix: kernel params")));
138  else
139  APparams= rcp(new ParameterList);
140 
141  // By default, we don't need global constants for A*P
142  APparams->set("compute global constants: temporaries",APparams->get("compute global constants: temporaries",false));
143  APparams->set("compute global constants",APparams->get("compute global constants",false));
144 
145  if (coarseLevel.IsAvailable("AP reuse data", this)) {
146  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
147 
148  APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
149 
150  if (APparams->isParameter("graph"))
151  AP = APparams->get< RCP<Matrix> >("graph");
152  }
153 
154  {
155  SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
156 
157  AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(Statistics2),
158  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::A*P-")+levelstr.str(), APparams);
159  }
160 
161  // Reuse coarse matrix memory if available (multiple solve)
162  RCP<ParameterList> RAPparams;
163  if(pL.isSublist("matrixmatrix: kernel params")) RAPparams=rcp(new ParameterList(pL.sublist("matrixmatrix: kernel params")));
164  else RAPparams= rcp(new ParameterList);
165 
166 
167 
168  if (coarseLevel.IsAvailable("RAP reuse data", this)) {
169  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous RAP data" << std::endl;
170 
171  RAPparams = coarseLevel.Get< RCP<ParameterList> >("RAP reuse data", this);
172 
173  if (RAPparams->isParameter("graph"))
174  Ac = RAPparams->get< RCP<Matrix> >("graph");
175 
176  // Some eigenvalue may have been cached with the matrix in the previous run.
177  // As the matrix values will be updated, we need to reset the eigenvalue.
178  Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
179  }
180 
181  // We *always* need global constants for the RAP, but not for the temps
182  RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
183  RAPparams->set("compute global constants",true);
184 
185  // Allow optimization of storage.
186  // This is necessary for new faster Epetra MM kernels.
187  // Seems to work with matrix modifications to repair diagonal entries.
188 
189  if (pL.get<bool>("transpose: use implicit") == true) {
190  SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
191 
192  Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
193  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
194 
195  } else {
196  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
197 
198  SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
199 
200  Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
201  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
202  }
203  bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
204  bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
205  if (checkAc || repairZeroDiagonals)
206  Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1));
207 
208  if (IsPrint(Statistics2)) {
209  RCP<ParameterList> params = rcp(new ParameterList());;
210  params->set("printLoadBalancingInfo", true);
211  params->set("printCommInfo", true);
212  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
213  }
214 
215  Set(coarseLevel, "A", Ac);
216 
217  APparams->set("graph", AP);
218  Set(coarseLevel, "AP reuse data", APparams);
219  RAPparams->set("graph", Ac);
220  Set(coarseLevel, "RAP reuse data", RAPparams);
221  } else {
222  RCP<ParameterList> RAPparams;
223  RAPparams= rcp(new ParameterList);
224 
225  // We *always* need global constants for the RAP, but not for the temps
226  RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
227  RAPparams->set("compute global constants",true);
228 
229  if (pL.get<bool>("transpose: use implicit") == true) {
230 
231  Ac = MatrixFactory::Build(P->getDomainMap(), Teuchos::as<LO>(0));
232 
233  SubFactoryMonitor m2(*this, "MxMxM: R x A x P (implicit)", coarseLevel);
234 
235  Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
236  MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
237  doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-implicit-")+levelstr.str(),
238  RAPparams);
239 
240  } else {
241  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
242  Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));
243 
244  SubFactoryMonitor m2(*this, "MxMxM: R x A x P (explicit)", coarseLevel);
245 
246  Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
247  MultiplyRAP(*R, !doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
248  doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-explicit-")+levelstr.str(),
249  RAPparams);
250  }
251  bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
252  bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
253  if (checkAc || repairZeroDiagonals)
254  Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1));
255 
256  if (IsPrint(Statistics2)) {
257  RCP<ParameterList> params = rcp(new ParameterList());;
258  params->set("printLoadBalancingInfo", true);
259  params->set("printCommInfo", true);
260  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
261  }
262 
263  Set(coarseLevel, "A", Ac);
264 
265  // RAPparams->set("graph", Ac);
266  // Set(coarseLevel, "RAP reuse data", RAPparams);
267  }
268 
269 
270  }
271 
272  if (transferFacts_.begin() != transferFacts_.end()) {
273  SubFactoryMonitor m(*this, "Projections", coarseLevel);
274 
275  // call Build of all user-given transfer factories
276  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
277  RCP<const FactoryBase> fac = *it;
278  GetOStream(Runtime0) << "RAPFactory: call transfer factory: " << fac->description() << std::endl;
279  fac->CallBuild(coarseLevel);
280  // Coordinates transfer is marginally different from all other operations
281  // because it is *optional*, and not required. For instance, we may need
282  // coordinates only on level 4 if we start repartitioning from that level,
283  // but we don't need them on level 1,2,3. As our current Hierarchy setup
284  // assumes propagation of dependencies only through three levels, this
285  // means that we need to rely on other methods to propagate optional data.
286  //
287  // The method currently used is through RAP transfer factories, which are
288  // simply factories which are called at the end of RAP with a single goal:
289  // transfer some fine data to coarser level. Because these factories are
290  // kind of outside of the mainline factories, they behave different. In
291  // particular, we call their Build method explicitly, rather than through
292  // Get calls. This difference is significant, as the Get call is smart
293  // enough to know when to release all factory dependencies, and Build is
294  // dumb. This led to the following CoordinatesTransferFactory sequence:
295  // 1. Request level 0
296  // 2. Request level 1
297  // 3. Request level 0
298  // 4. Release level 0
299  // 5. Release level 1
300  //
301  // The problem is missing "6. Release level 0". Because it was missing,
302  // we had outstanding request on "Coordinates", "Aggregates" and
303  // "CoarseMap" on level 0.
304  //
305  // This was fixed by explicitly calling Release on transfer factories in
306  // RAPFactory. I am still unsure how exactly it works, but now we have
307  // clear data requests for all levels.
308  coarseLevel.Release(*fac);
309  }
310  }
311 
312  }
313 
314  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
316  // check if it's a TwoLevelFactoryBase based transfer factory
317  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast,
318  "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. "
319  "This is very strange. (Note: you can remove this exception if there's a good reason for)");
320  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_, Exceptions::RuntimeError, "MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
321  transferFacts_.push_back(factory);
322  }
323 
324 } //namespace MueLu
325 
326 #define MUELU_RAPFACTORY_SHORT
327 #endif // MUELU_RAPFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
virtual void CallBuild(Level &requestedLevel) const =0
virtual std::string getObjectLabel() const
ParameterList & disableRecursiveValidation()
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
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)
bool isSublist(const std::string &name) const
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.
Print skeleton for the run, i.e. factory calls and used parameters.
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.
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.
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
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
virtual std::string description() const
Return a simple one-line description of this object.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.