MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
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 = rcp(new ParameterList);
136  if(pL.isSublist("matrixmatrix: kernel params"))
137  APparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
138 
139  // By default, we don't need global constants for A*P
140  APparams->set("compute global constants: temporaries",APparams->get("compute global constants: temporaries",false));
141  APparams->set("compute global constants",APparams->get("compute global constants",false));
142 
143  if (coarseLevel.IsAvailable("AP reuse data", this)) {
144  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
145 
146  APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
147 
148  if (APparams->isParameter("graph"))
149  AP = APparams->get< RCP<Matrix> >("graph");
150  }
151 
152  {
153  SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
154 
155  AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(Statistics2),
156  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::A*P-")+levelstr.str(), APparams);
157  }
158 
159  // Reuse coarse matrix memory if available (multiple solve)
160  RCP<ParameterList> RAPparams = rcp(new ParameterList);
161  if(pL.isSublist("matrixmatrix: kernel params"))
162  RAPparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
163 
164  if (coarseLevel.IsAvailable("RAP reuse data", this)) {
165  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous RAP data" << std::endl;
166 
167  RAPparams = coarseLevel.Get< RCP<ParameterList> >("RAP reuse data", this);
168 
169  if (RAPparams->isParameter("graph"))
170  Ac = RAPparams->get< RCP<Matrix> >("graph");
171 
172  // Some eigenvalue may have been cached with the matrix in the previous run.
173  // As the matrix values will be updated, we need to reset the eigenvalue.
174  Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
175  }
176 
177  // We *always* need global constants for the RAP, but not for the temps
178  RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
179  RAPparams->set("compute global constants",true);
180 
181  // Allow optimization of storage.
182  // This is necessary for new faster Epetra MM kernels.
183  // Seems to work with matrix modifications to repair diagonal entries.
184 
185  if (pL.get<bool>("transpose: use implicit") == true) {
186  SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
187 
188  Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
189  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
190 
191  } else {
192  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
193 
194  SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
195 
196  Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
197  doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
198  }
199  bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
200  bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
201  if (checkAc || repairZeroDiagonals)
202  Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1));
203 
204  if (IsPrint(Statistics2)) {
205  RCP<ParameterList> params = rcp(new ParameterList());;
206  params->set("printLoadBalancingInfo", true);
207  params->set("printCommInfo", true);
208  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
209  }
210 
211  if(!Ac.is_null()) {std::ostringstream oss; oss << "A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
212  Set(coarseLevel, "A", Ac);
213 
214  APparams->set("graph", AP);
215  Set(coarseLevel, "AP reuse data", APparams);
216  RAPparams->set("graph", Ac);
217  Set(coarseLevel, "RAP reuse data", RAPparams);
218  } else {
219  RCP<ParameterList> RAPparams = rcp(new ParameterList);
220  if(pL.isSublist("matrixmatrix: kernel params"))
221  RAPparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
222 
223  // We *always* need global constants for the RAP, but not for the temps
224  RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
225  RAPparams->set("compute global constants",true);
226 
227  if (pL.get<bool>("transpose: use implicit") == true) {
228 
229  Ac = MatrixFactory::Build(P->getDomainMap(), Teuchos::as<LO>(0));
230 
231  SubFactoryMonitor m2(*this, "MxMxM: R x A x P (implicit)", coarseLevel);
232 
234  MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
235  doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-implicit-")+levelstr.str(),
236  RAPparams);
237 
238  } else {
239  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
240  Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));
241 
242  SubFactoryMonitor m2(*this, "MxMxM: R x A x P (explicit)", coarseLevel);
243 
245  MultiplyRAP(*R, !doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
246  doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-explicit-")+levelstr.str(),
247  RAPparams);
248  }
249  bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
250  bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
251  if (checkAc || repairZeroDiagonals)
252  Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1));
253 
254  if (IsPrint(Statistics2)) {
255  RCP<ParameterList> params = rcp(new ParameterList());;
256  params->set("printLoadBalancingInfo", true);
257  params->set("printCommInfo", true);
258  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
259  }
260 
261  if(!Ac.is_null()) {std::ostringstream oss; oss << "A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
262  Set(coarseLevel, "A", Ac);
263 
264  // RAPparams->set("graph", Ac);
265  // Set(coarseLevel, "RAP reuse data", RAPparams);
266  }
267 
268 
269  }
270 
271  if (transferFacts_.begin() != transferFacts_.end()) {
272  SubFactoryMonitor m(*this, "Projections", coarseLevel);
273 
274  // call Build of all user-given transfer factories
275  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
276  RCP<const FactoryBase> fac = *it;
277  GetOStream(Runtime0) << "RAPFactory: call transfer factory: " << fac->description() << std::endl;
278  fac->CallBuild(coarseLevel);
279  // Coordinates transfer is marginally different from all other operations
280  // because it is *optional*, and not required. For instance, we may need
281  // coordinates only on level 4 if we start repartitioning from that level,
282  // but we don't need them on level 1,2,3. As our current Hierarchy setup
283  // assumes propagation of dependencies only through three levels, this
284  // means that we need to rely on other methods to propagate optional data.
285  //
286  // The method currently used is through RAP transfer factories, which are
287  // simply factories which are called at the end of RAP with a single goal:
288  // transfer some fine data to coarser level. Because these factories are
289  // kind of outside of the mainline factories, they behave different. In
290  // particular, we call their Build method explicitly, rather than through
291  // Get calls. This difference is significant, as the Get call is smart
292  // enough to know when to release all factory dependencies, and Build is
293  // dumb. This led to the following CoordinatesTransferFactory sequence:
294  // 1. Request level 0
295  // 2. Request level 1
296  // 3. Request level 0
297  // 4. Release level 0
298  // 5. Release level 1
299  //
300  // The problem is missing "6. Release level 0". Because it was missing,
301  // we had outstanding request on "Coordinates", "Aggregates" and
302  // "CoarseMap" on level 0.
303  //
304  // This was fixed by explicitly calling Release on transfer factories in
305  // RAPFactory. I am still unsure how exactly it works, but now we have
306  // clear data requests for all levels.
307  coarseLevel.Release(*fac);
308  }
309  }
310 
311  }
312 
313  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
315  // check if it's a TwoLevelFactoryBase based transfer factory
316  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast,
317  "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. "
318  "This is very strange. (Note: you can remove this exception if there's a good reason for)");
319  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_, Exceptions::RuntimeError, "MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
320  transferFacts_.push_back(factory);
321  }
322 
323 } //namespace MueLu
324 
325 #define MUELU_RAPFACTORY_SHORT
326 #endif // MUELU_RAPFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
virtual std::string getObjectLabel() const
virtual void CallBuild(Level &requestedLevel) const =0
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)
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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.
static void CheckRepairMainDiagonal(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &Ac, bool const &repairZeroDiagonals, Teuchos::FancyOStream &fos, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType threshold=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< Scalar >::magnitudeType >::zero())
static std::string getColonLabel(const std::string &label)
Helper function for object label.
Print even more statistics.
Additional warnings.
bool isParameter(const std::string &name) const
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
bool isSublist(const std::string &name) const
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)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
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.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.