MueLu  Version of the Day
MueLu_CreateEpetraPreconditioner.hpp
Go to the documentation of this file.
1 #ifndef MUELU_CREATE_EPETRA_PRECONDITIONER_HPP
2 #define MUELU_CREATE_EPETRA_PRECONDITIONER_HPP
3 
4 #include <Teuchos_XMLParameterListHelpers.hpp>
5 #include <Xpetra_CrsMatrix.hpp>
6 #include <Xpetra_MultiVector.hpp>
7 #include <Xpetra_MultiVectorFactory.hpp>
8 
9 #include <MueLu.hpp>
10 
11 #include <MueLu_EpetraOperator.hpp>
12 #include <MueLu_Exceptions.hpp>
13 #include <MueLu_Hierarchy.hpp>
14 #include <MueLu_MasterList.hpp>
15 #include <MueLu_MLParameterListInterpreter.hpp>
16 #include <MueLu_ParameterListInterpreter.hpp>
17 #include <MueLu_Utilities.hpp>
18 #include <MueLu_HierarchyHelpers.hpp>
19 
22 #if defined(HAVE_MUELU_EPETRA)
23 namespace MueLu {
24 
34  Teuchos::RCP<MueLu::EpetraOperator>
35  CreateEpetraPreconditioner(const Teuchos::RCP<Epetra_CrsMatrix>& inA,
36  // FIXME: why is it non-const
37  Teuchos::ParameterList& paramListIn,
38  const Teuchos::RCP<Epetra_MultiVector>& inCoords = Teuchos::null,
39  const Teuchos::RCP<Epetra_MultiVector>& inNullspace = Teuchos::null)
40  {
41  typedef double SC;
42  typedef int LO;
43  typedef int GO;
44  typedef Xpetra::EpetraNode NO;
45 
46  using Teuchos::ParameterList;
47 
48  typedef Xpetra::MultiVector<SC, LO, GO, NO> MultiVector;
49  typedef Xpetra::Matrix<SC, LO, GO, NO> Matrix;
52 
53  bool hasParamList = paramListIn.numParams();
54 
55  RCP<HierarchyManager> mueLuFactory;
56  ParameterList paramList = paramListIn;
57 
58  std::string syntaxStr = "parameterlist: syntax";
59  if (hasParamList && paramList.isParameter(syntaxStr) && paramList.get<std::string>(syntaxStr) == "ml") {
60  paramList.remove(syntaxStr);
61  mueLuFactory = rcp(new MLParameterListInterpreter<SC,LO,GO,NO>(paramList));
62 
63  } else {
64  mueLuFactory = rcp(new ParameterListInterpreter <SC,LO,GO,NO>(paramList,Xpetra::toXpetra(inA->Comm())));
65  }
66 
67  RCP<Hierarchy> H = mueLuFactory->CreateHierarchy();
68  H->setlib(Xpetra::UseEpetra);
69 
70  // Wrap A
71  RCP<Matrix> A = EpetraCrs_To_XpetraMatrix<SC, LO, GO, NO>(inA);
72  H->GetLevel(0)->Set("A", A);
73 
74  // Wrap coordinates if available
75  if (inCoords != Teuchos::null) {
76  RCP<MultiVector> coordinates = EpetraMultiVector_To_XpetraMultiVector<SC,LO,GO,NO>(inCoords);
77  H->GetLevel(0)->Set("Coordinates", coordinates);
78  }
79 
80  // Wrap nullspace if available, otherwise use constants
81  RCP<MultiVector> nullspace;
82  if (inNullspace != Teuchos::null) {
83  nullspace = EpetraMultiVector_To_XpetraMultiVector<SC, LO, GO, NO>(inNullspace);
84 
85  } else {
86  int nPDE = MasterList::getDefault<int>("number of equations");
87  if (paramList.isSublist("Matrix")) {
88  // Factory style parameter list
89  const Teuchos::ParameterList& operatorList = paramList.sublist("Matrix");
90  if (operatorList.isParameter("PDE equations"))
91  nPDE = operatorList.get<int>("PDE equations");
92 
93  } else if (paramList.isParameter("number of equations")) {
94  // Easy style parameter list
95  nPDE = paramList.get<int>("number of equations");
96  }
97 
98  nullspace = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(A->getDomainMap(), nPDE);
99  if (nPDE == 1) {
100  nullspace->putScalar(Teuchos::ScalarTraits<SC>::one());
101 
102  } else {
103  for (int i = 0; i < nPDE; i++) {
104  Teuchos::ArrayRCP<SC> nsData = nullspace->getDataNonConst(i);
105  for (int j = 0; j < nsData.size(); j++) {
106  GO GID = A->getDomainMap()->getGlobalElement(j) - A->getDomainMap()->getIndexBase();
107 
108  if ((GID-i) % nPDE == 0)
109  nsData[j] = Teuchos::ScalarTraits<SC>::one();
110  }
111  }
112  }
113  }
114  H->GetLevel(0)->Set("Nullspace", nullspace);
115 
116  Teuchos::ParameterList nonSerialList,dummyList;
117  ExtractNonSerializableData(paramList, dummyList, nonSerialList);
119 
120  mueLuFactory->SetupHierarchy(*H);
121 
122  return rcp(new EpetraOperator(H));
123  }
124 
133  Teuchos::RCP<MueLu::EpetraOperator>
134  CreateEpetraPreconditioner(const Teuchos::RCP<Epetra_CrsMatrix> & inA,
135  const Teuchos::RCP<Epetra_MultiVector>& inCoords = Teuchos::null,
136  const Teuchos::RCP<Epetra_MultiVector>& inNullspace = Teuchos::null) {
137  Teuchos::ParameterList paramList;
138  return CreateEpetraPreconditioner(inA, paramList, inCoords, inNullspace);
139  }
140 
150  Teuchos::RCP<MueLu::EpetraOperator>
151  CreateEpetraPreconditioner(const Teuchos::RCP<Epetra_CrsMatrix> & A,
152  const std::string& xmlFileName,
153  const Teuchos::RCP<Epetra_MultiVector>& inCoords = Teuchos::null,
154  const Teuchos::RCP<Epetra_MultiVector>& inNullspace = Teuchos::null)
155  {
156  Teuchos::ParameterList paramList;
157  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *Xpetra::toXpetra(A->Comm()));
158 
159  return CreateEpetraPreconditioner(A, paramList, inCoords, inNullspace);
160  }
161 
162  void ReuseEpetraPreconditioner(const Teuchos::RCP<Epetra_CrsMatrix>& inA, MueLu::EpetraOperator& Op) {
163  typedef double SC;
164  typedef int LO;
165  typedef int GO;
166  typedef Xpetra::EpetraNode NO;
167 
168  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
169  typedef Xpetra::Operator<SC,LO,GO,NO> Operator;
171 
172  RCP<Hierarchy> H = Op.GetHierarchy();
173 
174  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), Exceptions::RuntimeError,
175  "ReuseTpetraPreconditioner: Hierarchy has no levels in it");
176  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), Exceptions::RuntimeError,
177  "ReuseTpetraPreconditioner: Hierarchy has no fine level operator");
178  RCP<Level> level0 = H->GetLevel(0);
179 
180  RCP<Operator> O0 = level0->Get<RCP<Operator> >("A");
181  RCP<Matrix> A0 = Teuchos::rcp_dynamic_cast<Matrix>(O0);
182 
183  RCP<Matrix> A = EpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(inA);
184  if (!A0.is_null()) {
185  // If a user provided a "number of equations" argument in a parameter list
186  // during the initial setup, we must honor that settings and reuse it for
187  // all consequent setups.
188  A->SetFixedBlockSize(A0->GetFixedBlockSize());
189  }
190  level0->Set("A", A);
191  H->SetupRe();
192  }
193 
194 } //namespace
195 #endif // HAVE_MUELU_SERIAL and HAVE_MUELU_EPETRA
196 
197 #endif //ifndef MUELU_CREATE_EPETRA_PRECONDITIONER_HPP
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
void ReuseEpetraPreconditioner(const Teuchos::RCP< Epetra_CrsMatrix > &inA, MueLu::EpetraOperator &Op)
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
Namespace for MueLu classes and methods.
Teuchos::RCP< MueLu::EpetraOperator > CreateEpetraPreconditioner(const Teuchos::RCP< Epetra_CrsMatrix > &inA, Teuchos::ParameterList &paramListIn, const Teuchos::RCP< Epetra_MultiVector > &inCoords=Teuchos::null, const Teuchos::RCP< Epetra_MultiVector > &inNullspace=Teuchos::null)
Helper function to create a MueLu preconditioner that can be used by Epetra.Given a EpetraCrs_Matrix...
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.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
long ExtractNonSerializableData(const Teuchos::ParameterList &inList, Teuchos::ParameterList &serialList, Teuchos::ParameterList &nonSerialList)