MueLu  Version of the Day
MueLu_CreateTpetraPreconditioner.hpp
Go to the documentation of this file.
1 #ifndef MUELU_CREATE_TPETRA_PRECONDITIONER_HPP
2 #define MUELU_CREATE_TPETRA_PRECONDITIONER_HPP
3 
6 
7 #include <Teuchos_XMLParameterListHelpers.hpp>
8 #include <Tpetra_Operator.hpp>
9 #include <Tpetra_RowMatrix.hpp>
10 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
11 #include <Tpetra_Experimental_BlockCrsMatrix.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_TpetraOperator.hpp>
24 #include <MueLu_Utilities.hpp>
25 #include <MueLu_HierarchyHelpers.hpp>
26 
27 
28 #if defined(HAVE_MUELU_EXPERIMENTAL) and defined(HAVE_MUELU_AMGX)
29 #include <MueLu_AMGXOperator.hpp>
30 #include <amgx_c.h>
31 #include "cuda_runtime.h"
32 #endif
33 
34 namespace MueLu {
35 
36 
46  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
47  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
48  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
49  Teuchos::ParameterList& inParamList,
50  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
51  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
52  {
53  typedef Scalar SC;
54  typedef LocalOrdinal LO;
55  typedef GlobalOrdinal GO;
56  typedef Node NO;
57 
58  using Teuchos::ParameterList;
59 
60  typedef Xpetra::MultiVector<SC,LO,GO,NO> MultiVector;
61  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
64  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> crs_matrix_type;
65  typedef Tpetra::Experimental::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> block_crs_matrix_type;
66 
67  bool hasParamList = inParamList.numParams();
68 
69  RCP<HierarchyManager> mueLuFactory;
70  ParameterList paramList = inParamList;
71  RCP<const crs_matrix_type> constCrsA;
72  RCP<crs_matrix_type> crsA;
73 
74 #if defined(HAVE_MUELU_EXPERIMENTAL) and defined(HAVE_MUELU_AMGX)
75  std::string externalMG = "use external multigrid package";
76  if (hasParamList && paramList.isParameter(externalMG) && paramList.get<std::string>(externalMG) == "amgx"){
77  constCrsA = rcp_dynamic_cast<const crs_matrix_type>(inA);
78  TEUCHOS_TEST_FOR_EXCEPTION(constCrsA == Teuchos::null, Exceptions::RuntimeError, "CreateTpetraPreconditioner: failed to dynamic cast to Tpetra::CrsMatrix, which is required to be able to use AmgX.");
79  return rcp(new AMGXOperator<SC,LO,GO,NO>(inA,inParamList));
80  }
81 #endif
82  std::string syntaxStr = "parameterlist: syntax";
83  if (hasParamList && paramList.isParameter(syntaxStr) && paramList.get<std::string>(syntaxStr) == "ml") {
84  paramList.remove(syntaxStr);
85  mueLuFactory = rcp(new MLParameterListInterpreter<SC,LO,GO,NO>(paramList));
86 
87  } else {
88  mueLuFactory = rcp(new ParameterListInterpreter <SC,LO,GO,NO>(paramList,inA->getDomainMap()->getComm()));
89  }
90 
91  RCP<Hierarchy> H = mueLuFactory->CreateHierarchy();
92  H->setlib(Xpetra::UseTpetra);
93 
94  // Wrap A
95  RCP<Matrix> A;
96  RCP<block_crs_matrix_type> bcrsA = rcp_dynamic_cast<block_crs_matrix_type>(inA);
97  crsA = rcp_dynamic_cast<crs_matrix_type>(inA);
98  if (crsA != Teuchos::null)
99  A = TpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(crsA);
100  else if (bcrsA != Teuchos::null) {
101  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > temp = rcp(new Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>(bcrsA));
102  TEUCHOS_TEST_FOR_EXCEPTION(temp==Teuchos::null, Exceptions::RuntimeError, "CreateTpetraPreconditioner: cast from Tpetra::Experimental::BlockCrsMatrix to Xpetra::TpetraBlockCrsMatrix failed.");
103  A = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(temp));
104  }
105  else {
106  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "CreateTpetraPreconditioner: only Tpetra CrsMatrix and BlockCrsMatrix types are supported.");
107  }
108  H->GetLevel(0)->Set("A", A);
109 
110  // Wrap coordinates if available
111  if (inCoords != Teuchos::null) {
112  RCP<Xpetra::MultiVector<double,LO,GO,NO> > coordinates = TpetraMultiVector_To_XpetraMultiVector<double,LO,GO,NO>(inCoords);
113  H->GetLevel(0)->Set("Coordinates", coordinates);
114  }
115 
116  // Wrap nullspace if available, otherwise use constants
117  RCP<MultiVector> nullspace;
118  if (inNullspace != Teuchos::null) {
119  nullspace = TpetraMultiVector_To_XpetraMultiVector<SC,LO,GO,NO>(inNullspace);
120 
121  } else {
122  int nPDE = MasterList::getDefault<int>("number of equations");
123  if (paramList.isSublist("Matrix")) {
124  // Factory style parameter list
125  const Teuchos::ParameterList& operatorList = paramList.sublist("Matrix");
126  if (operatorList.isParameter("PDE equations"))
127  nPDE = operatorList.get<int>("PDE equations");
128 
129  } else if (paramList.isParameter("number of equations")) {
130  // Easy style parameter list
131  nPDE = paramList.get<int>("number of equations");
132  }
133 
134  nullspace = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(A->getDomainMap(), nPDE);
135  if (nPDE == 1) {
136  nullspace->putScalar(Teuchos::ScalarTraits<SC>::one());
137 
138  } else {
139  for (int i = 0; i < nPDE; i++) {
140  Teuchos::ArrayRCP<SC> nsData = nullspace->getDataNonConst(i);
141  for (int j = 0; j < nsData.size(); j++) {
142  GO GID = A->getDomainMap()->getGlobalElement(j) - A->getDomainMap()->getIndexBase();
143 
144  if ((GID-i) % nPDE == 0)
145  nsData[j] = Teuchos::ScalarTraits<SC>::one();
146  }
147  }
148  }
149  }
150  H->GetLevel(0)->Set("Nullspace", nullspace);
151 
152 
153  Teuchos::ParameterList nonSerialList,dummyList;
154  ExtractNonSerializableData(paramList, dummyList, nonSerialList);
156 
157  mueLuFactory->SetupHierarchy(*H);
158  return rcp(new TpetraOperator<SC,LO,GO,NO>(H));
159  }
160 
161 
174  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175  MUELU_DEPRECATED
176  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
177  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
178  Teuchos::ParameterList& inParamList,
179  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
180  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
181  {
182  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
183  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, inParamList, inCoords, inNullspace);
184  }
185 
186 
198  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
199  MUELU_DEPRECATED
200  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
201  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix <Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
202  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
203  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
204  {
205  Teuchos::ParameterList paramList;
206  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
207  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, paramList, inCoords, inNullspace);
208  }
209 
210 
223  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
224  MUELU_DEPRECATED
225  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
226  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix <Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
227  const std::string& xmlFileName,
228  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
229  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
230  {
231  Teuchos::ParameterList paramList;
232  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getComm());
233 
234  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
235  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, paramList, inCoords, inNullspace);
236  }
237 
238 
249  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
250  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
251  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
252  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
253  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
254  {
255  Teuchos::ParameterList paramList;
256  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
257  }
258 
259 
271  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
272  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
273  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
274  const std::string& xmlFileName,
275  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
276  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
277  {
278  Teuchos::ParameterList paramList;
279  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getDomainMap()->getComm());
280  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
281  }
282 
283 
291  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
292  void ReuseTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
294  typedef Scalar SC;
295  typedef LocalOrdinal LO;
296  typedef GlobalOrdinal GO;
297  typedef Node NO;
298 
299  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
300  typedef Xpetra::Operator<SC,LO,GO,NO> Operator;
302 
303  RCP<Hierarchy> H = Op.GetHierarchy();
304 
305  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), Exceptions::RuntimeError,
306  "MueLu::ReuseTpetraPreconditioner: Hierarchy has no levels in it");
307  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), Exceptions::RuntimeError,
308  "MueLu::ReuseTpetraPreconditioner: Hierarchy has no fine level operator");
309  RCP<Level> level0 = H->GetLevel(0);
310 
311  RCP<Operator> O0 = level0->Get<RCP<Operator> >("A");
312  RCP<Matrix> A0 = Teuchos::rcp_dynamic_cast<Matrix>(O0);
313 
314  RCP<Matrix> A = TpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(inA);
315  if (!A0.is_null()) {
316  // If a user provided a "number of equations" argument in a parameter list
317  // during the initial setup, we must honor that settings and reuse it for
318  // all consequent setups.
319  A->SetFixedBlockSize(A0->GetFixedBlockSize());
320  }
321  level0->Set("A", A);
322 
323  H->SetupRe();
324  }
325 
326 } //namespace
327 
328 #endif //ifndef MUELU_CREATE_TPETRA_PRECONDITIONER_HPP
329 
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
Namespace for MueLu classes and methods.
Teuchos::RCP< MueLu::TpetraOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateTpetraPreconditioner(const Teuchos::RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &inA, Teuchos::ParameterList &inParamList, const Teuchos::RCP< Tpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node >> &inCoords=Teuchos::null, const Teuchos::RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &inNullspace=Teuchos::null)
Helper function to create a MueLu or AMGX preconditioner that can be used by Tpetra.Given a Tpetra::Operator, this function returns a constructed MueLu preconditioner.
Class that accepts ML-style parameters and builds a MueLu preconditioner. This interpreter uses the s...
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetHierarchy() const
Direct access to the underlying MueLu::Hierarchy.
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
void ReuseTpetraPreconditioner(const Teuchos::RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &inA, MueLu::TpetraOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Helper function to reuse an existing MueLu preconditioner.
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)
Adapter for AmgX library from Nvidia.