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_CreateXpetraPreconditioner.hpp>
25 #include <MueLu_Utilities.hpp>
26 #include <MueLu_HierarchyUtils.hpp>
27 
28 
29 #if defined(HAVE_MUELU_AMGX)
30 #include <MueLu_AMGXOperator.hpp>
31 #include <amgx_c.h>
32 #include "cuda_runtime.h"
33 #endif
34 
35 namespace MueLu {
36 
37 
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  Teuchos::ParameterList& dummyList)
51  {
52  typedef Scalar SC;
53  typedef LocalOrdinal LO;
54  typedef GlobalOrdinal GO;
55  typedef Node NO;
56 
57  using Teuchos::ParameterList;
58 
59  typedef Xpetra::MultiVector<SC,LO,GO,NO> MultiVector;
60  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
61  typedef Hierarchy<SC,LO,GO,NO> Hierarchy;
62  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> crs_matrix_type;
63  typedef Tpetra::Experimental::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> block_crs_matrix_type;
64 
65 #if defined(HAVE_MUELU_AMGX)
66  std::string externalMG = "use external multigrid package";
67  if (inParamList.isParameter(externalMG) && inParamList.get<std::string>(externalMG) == "amgx"){
68  const RCP<crs_matrix_type> constCrsA = rcp_dynamic_cast<crs_matrix_type>(inA);
69  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.");
70  return rcp(new AMGXOperator<SC,LO,GO,NO>(constCrsA,inParamList));
71  }
72 #endif
73 
74  // Wrap A
75  RCP<Matrix> A;
76  RCP<block_crs_matrix_type> bcrsA = rcp_dynamic_cast<block_crs_matrix_type>(inA);
77  RCP<crs_matrix_type> crsA = rcp_dynamic_cast<crs_matrix_type>(inA);
78  if (crsA != Teuchos::null)
79  A = TpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(crsA);
80  else if (bcrsA != Teuchos::null) {
81  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > temp = rcp(new Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>(bcrsA));
82  TEUCHOS_TEST_FOR_EXCEPTION(temp==Teuchos::null, Exceptions::RuntimeError, "CreateTpetraPreconditioner: cast from Tpetra::Experimental::BlockCrsMatrix to Xpetra::TpetraBlockCrsMatrix failed.");
83  A = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(temp));
84  }
85  else {
86  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "CreateTpetraPreconditioner: only Tpetra CrsMatrix and BlockCrsMatrix types are supported.");
87  }
88 
89  Teuchos::ParameterList& userList = inParamList.sublist("user data");
90  if (userList.isParameter("Coordinates")) {
91  RCP<Xpetra::MultiVector<double,LO,GO,NO> > coordinates = Teuchos::null;
92  try {
93  coordinates = TpetraMultiVector_To_XpetraMultiVector<double,LO,GO,NO>(userList.get<RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> > >("Coordinates"));
94  } catch(Teuchos::Exceptions::InvalidParameterType) {
95  coordinates = userList.get<RCP<Xpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> > >("Coordinates");
96  }
97  userList.set<RCP<Xpetra::MultiVector<double,LO,GO,NO> > >("Coordinates", coordinates);
98  }
99 
100  if (userList.isParameter("Nullspace")) {
101  RCP<MultiVector> nullspace = Teuchos::null;
102  try {
103  nullspace = TpetraMultiVector_To_XpetraMultiVector<SC,LO,GO,NO>(userList.get<RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >("Nullspace"));
104  } catch(Teuchos::Exceptions::InvalidParameterType) {
105  nullspace = userList.get<RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >("Nullspace");
106  }
107  userList.set<RCP<MultiVector> >("Nullspace", nullspace);
108  }
109 
110  RCP<Hierarchy> H = MueLu::CreateXpetraPreconditioner<SC,LO,GO,NO>(A,inParamList,inParamList);
111  return rcp(new TpetraOperator<SC,LO,GO,NO>(H));
112  }
113 
114 
124  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
125  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
126  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
127  Teuchos::ParameterList& inParamList,
128  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
129  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
130  {
131  typedef Scalar SC;
132  typedef LocalOrdinal LO;
133  typedef GlobalOrdinal GO;
134  typedef Node NO;
135 
136  using Teuchos::ParameterList;
137 
138  // Here the assumption is that the nullspace and coordinates passed
139  // as multivectors will overwrite data on the parameter list.
140  // If you want the data on the parameter list to prevail, then call the
141  // specialization that only accept an operator and a parameterList
142  Teuchos::ParameterList& userList = inParamList.sublist("user data");
143  if (inCoords != Teuchos::null) {
144  userList.set<RCP<Tpetra::MultiVector<double,LO,GO,NO> > >("Coordinates", inCoords);
145  }
146  if (inNullspace != Teuchos::null) {
147  userList.set<RCP<Tpetra::MultiVector<SC,LO,GO,NO> > >("Nullspace", inNullspace);
148  }
149 
150  return CreateTpetraPreconditioner<SC,LO,GO,NO>(inA,inParamList,inParamList);
151  }
152 
153 
166  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
167  MUELU_DEPRECATED
168  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
169  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
170  Teuchos::ParameterList& inParamList,
171  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
172  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
173  {
174  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
175  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, inParamList, inCoords, inNullspace);
176  }
177 
178 
190  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
191  MUELU_DEPRECATED
192  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
193  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix <Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
194  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
195  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
196  {
197  Teuchos::ParameterList paramList;
198  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
199  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, paramList, inCoords, inNullspace);
200  }
201 
202 
215  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
216  MUELU_DEPRECATED
217  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
218  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix <Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
219  const std::string& xmlFileName,
220  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
221  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
222  {
223  Teuchos::ParameterList paramList;
224  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getComm());
225 
226  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
227  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, paramList, inCoords, inNullspace);
228  }
229 
230 
241  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
242  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
243  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
244  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
245  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
246  {
247  Teuchos::ParameterList paramList;
248  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
249  }
250 
251 
263  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
264  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
265  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
266  const std::string& xmlFileName,
267  const Teuchos::RCP<Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
268  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
269  {
270  Teuchos::ParameterList paramList;
271  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getDomainMap()->getComm());
272  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
273  }
274 
275 
283  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
284  void ReuseTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
285  MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
286  typedef Scalar SC;
287  typedef LocalOrdinal LO;
288  typedef GlobalOrdinal GO;
289  typedef Node NO;
290 
291  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
292  typedef MueLu ::Hierarchy<SC,LO,GO,NO> Hierarchy;
293 
294  RCP<Hierarchy> H = Op.GetHierarchy();
295  RCP<Matrix> A = TpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(inA);
296 
297  MueLu::ReuseXpetraPreconditioner<SC,LO,GO,NO>(A, H);
298  }
299 
300 } //namespace
301 
302 #endif //ifndef MUELU_CREATE_TPETRA_PRECONDITIONER_HPP
303 
Xpetra_MultiVector.hpp
Xpetra::TpetraBlockCrsMatrix
Xpetra_MultiVectorFactory.hpp
MueLu::CreateTpetraPreconditioner
Teuchos::RCP< MueLu::TpetraOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateTpetraPreconditioner(const Teuchos::RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &inA, Teuchos::ParameterList &inParamList, Teuchos::ParameterList &dummyList)
Helper function to create a MueLu or AMGX preconditioner that can be used by Tpetra.
Definition: MueLu_CreateTpetraPreconditioner.hpp:48
MueLu_CreateXpetraPreconditioner.hpp
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
Xpetra::MultiVector
NO
Node NO
MueLu_AMGXOperator.hpp
MueLu::AMGXOperator
Adapter for AmgX library from Nvidia.
Definition: MueLu_AMGXOperator_fwd.hpp:54
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Tpetra::Experimental::Classes::BlockCrsMatrix
Teuchos::Exceptions::InvalidParameterType
MultiVector
Classes::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
Xpetra_TpetraBlockCrsMatrix.hpp
SC
Scalar SC
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
Teuchos::RCP
MueLu::ReuseTpetraPreconditioner
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.
Definition: MueLu_CreateTpetraPreconditioner.hpp:284
Teuchos::Ptr
LO
LocalOrdinal LO
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Xpetra::Matrix
Tpetra::Classes::CrsMatrix
MueLu.hpp
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
MueLu_Exceptions.hpp
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::TpetraOperator::GetHierarchy
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetHierarchy() const
Direct access to the underlying MueLu::Hierarchy.
Definition: MueLu_TpetraOperator_def.hpp:127
Tpetra::Classes::MultiVector
Xpetra::CrsMatrixWrap
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
MueLu::TpetraOperator
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
Definition: MueLu_TpetraOperator_fwd.hpp:54
GO
GlobalOrdinal GO
Teuchos::ParameterList
MueLu_MasterList.hpp
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:104
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Tpetra::Classes::Operator
Xpetra_CrsMatrix.hpp
Teuchos_XMLParameterListHelpers.hpp