MueLu  Version of the Day
MueLu_CreateXpetraPreconditioner.hpp
Go to the documentation of this file.
1 /*
2  * MueLu_CreateXpetraPreconditioner.hpp
3  *
4  * Created on: Feb 5, 2016
5  * Author: tawiesn
6  */
7 
8 #ifndef PACKAGES_MUELU_ADAPTERS_XPETRA_MUELU_CREATEXPETRAPRECONDITIONER_HPP_
9 #define PACKAGES_MUELU_ADAPTERS_XPETRA_MUELU_CREATEXPETRAPRECONDITIONER_HPP_
10 
11 #include <Teuchos_XMLParameterListHelpers.hpp>
12 #include <Xpetra_CrsMatrix.hpp>
13 #include <Xpetra_MultiVector.hpp>
14 #include <Xpetra_MultiVectorFactory.hpp>
15 
16 #include <MueLu.hpp>
17 
18 #include <MueLu_Exceptions.hpp>
19 #include <MueLu_Hierarchy.hpp>
20 #include <MueLu_MasterList.hpp>
21 #include <MueLu_MLParameterListInterpreter.hpp>
22 #include <MueLu_ParameterListInterpreter.hpp>
23 #include <MueLu_Utilities.hpp>
24 #include <MueLu_HierarchyUtils.hpp>
25 
26 namespace MueLu {
27 
28  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
29  Teuchos::RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
30  CreateXpetraPreconditioner(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > op,
31  const Teuchos::ParameterList& inParamList,
32  Teuchos::RCP<Xpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> > coords = Teuchos::null,
33  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > nullspace = Teuchos::null) {
39  typedef Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVectorFactory;
40 
41  std::string timerName = "MueLu setup time";
42  RCP<Teuchos::Time> tm = Teuchos::TimeMonitor::getNewTimer(timerName);
43  tm->start();
44 
45  ParameterList paramList = inParamList;
46 
47  bool hasParamList = paramList.numParams();
48 
49  RCP<HierarchyManager> mueLuFactory;
50 
51  std::string syntaxStr = "parameterlist: syntax";
52  if (hasParamList && paramList.isParameter(syntaxStr) && paramList.get<std::string>(syntaxStr) == "ml") {
53  paramList.remove(syntaxStr);
54  mueLuFactory = rcp(new MLParameterListInterpreter(paramList));
55  } else {
56  mueLuFactory = rcp(new ParameterListInterpreter(paramList,op->getDomainMap()->getComm()));
57  }
58 
59  RCP<Hierarchy> H = mueLuFactory->CreateHierarchy();
60  H->setlib(op->getDomainMap()->lib());
61 
62 
63  // Set fine level operator
64  H->GetLevel(0)->Set("A", op);
65 
66  // Set coordinates if available
67  if (coords != Teuchos::null) {
68  H->GetLevel(0)->Set("Coordinates", coords);
69  }
70 
71  // Wrap nullspace if available, otherwise use constants
72  if (nullspace == Teuchos::null) {
73  int nPDE = MueLu::MasterList::getDefault<int>("number of equations");
74  if (paramList.isSublist("Matrix")) {
75  // Factory style parameter list
76  const Teuchos::ParameterList& operatorList = paramList.sublist("Matrix");
77  if (operatorList.isParameter("PDE equations"))
78  nPDE = operatorList.get<int>("PDE equations");
79 
80  } else if (paramList.isParameter("number of equations")) {
81  // Easy style parameter list
82  nPDE = paramList.get<int>("number of equations");
83  }
84 
85  nullspace = MultiVectorFactory::Build(op->getDomainMap(), nPDE);
86  if (nPDE == 1) {
87  nullspace->putScalar(Teuchos::ScalarTraits<Scalar>::one());
88 
89  } else {
90  for (int i = 0; i < nPDE; i++) {
91  Teuchos::ArrayRCP<Scalar> nsData = nullspace->getDataNonConst(i);
92  for (int j = 0; j < nsData.size(); j++) {
93  GlobalOrdinal GID = op->getDomainMap()->getGlobalElement(j) - op->getDomainMap()->getIndexBase();
94 
95  if ((GID-i) % nPDE == 0)
96  nsData[j] = Teuchos::ScalarTraits<Scalar>::one();
97  }
98  }
99  }
100  }
101  H->GetLevel(0)->Set("Nullspace", nullspace);
102 
103 
104  Teuchos::ParameterList nonSerialList,dummyList;
105  MueLu::ExtractNonSerializableData(paramList, dummyList, nonSerialList);
106  HierarchyUtils::AddNonSerializableDataToHierarchy(*mueLuFactory,*H, nonSerialList);
107 
108  mueLuFactory->SetupHierarchy(*H);
109 
110  tm->stop();
111  tm->incrementNumCalls();
112 
113  if (H->GetVerbLevel() & Statistics0) {
114  const bool alwaysWriteLocal = true;
115  const bool writeGlobalStats = true;
116  const bool writeZeroTimers = false;
117  const bool ignoreZeroTimers = true;
118  const std::string filter = timerName;
119  Teuchos::TimeMonitor::summarize(op->getRowMap()->getComm().ptr(), std::cout, alwaysWriteLocal, writeGlobalStats,
120  writeZeroTimers, Teuchos::Union, filter, ignoreZeroTimers);
121  }
122 
123  tm->reset();
124 
125  return H;
126  }
127 
135  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
136  void ReuseXpetraPreconditioner(const Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A,
138  std::string timerName = "MueLu setup time";
139  RCP<Teuchos::Time> tm = Teuchos::TimeMonitor::getNewTimer(timerName);
140  tm->start();
141 
142  typedef Scalar SC;
143  typedef LocalOrdinal LO;
144  typedef GlobalOrdinal GO;
145  typedef Node NO;
146 
147  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
148  typedef Xpetra::Operator<SC,LO,GO,NO> Operator;
149 
150  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), Exceptions::RuntimeError,
151  "MueLu::ReuseXpetraPreconditioner: Hierarchy has no levels in it");
152  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), Exceptions::RuntimeError,
153  "MueLu::ReuseXpetraPreconditioner: Hierarchy has no fine level operator");
154  RCP<Level> level0 = H->GetLevel(0);
155 
156  RCP<Operator> O0 = level0->Get<RCP<Operator> >("A");
157  RCP<Matrix> A0 = Teuchos::rcp_dynamic_cast<Matrix>(O0);
158 
159  if (!A0.is_null()) {
160  // If a user provided a "number of equations" argument in a parameter list
161  // during the initial setup, we must honor that settings and reuse it for
162  // all consequent setups.
163  A->SetFixedBlockSize(A0->GetFixedBlockSize());
164  }
165  level0->Set("A", A);
166 
167  H->SetupRe();
168 
169  tm->stop();
170  tm->incrementNumCalls();
171 
172  if (H->GetVerbLevel() & Statistics0) {
173  const bool alwaysWriteLocal = true;
174  const bool writeGlobalStats = true;
175  const bool writeZeroTimers = false;
176  const bool ignoreZeroTimers = true;
177  const std::string filter = timerName;
178  Teuchos::TimeMonitor::summarize(A->getRowMap()->getComm().ptr(), std::cout, alwaysWriteLocal, writeGlobalStats,
179  writeZeroTimers, Teuchos::Union, filter, ignoreZeroTimers);
180  }
181 
182  tm->reset();
183  }
184 
185 } //namespace
186 
187 #endif /* PACKAGES_MUELU_ADAPTERS_XPETRA_MUELU_CREATEXPETRAPRECONDITIONER_HPP_ */
void ReuseXpetraPreconditioner(const Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &H)
Helper function to reuse an existing MueLu preconditioner.
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList, Teuchos::RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > coords=Teuchos::null, Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > nullspace=Teuchos::null)
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
Class that accepts ML-style parameters and builds a MueLu preconditioner. This interpreter uses the s...
Exception throws to report errors in the internal logical of the program.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
long ExtractNonSerializableData(const Teuchos::ParameterList &inList, Teuchos::ParameterList &serialList, Teuchos::ParameterList &nonSerialList)