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_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 
45  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
46  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
47  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
48  Teuchos::ParameterList& inParamList)
49  {
50  typedef Scalar SC;
51  typedef LocalOrdinal LO;
52  typedef GlobalOrdinal GO;
53  typedef Node NO;
54 
55  using Teuchos::ParameterList;
56 
57  typedef Xpetra::MultiVector<SC,LO,GO,NO> MultiVector;
58  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
59  typedef Hierarchy<SC,LO,GO,NO> Hierarchy;
60  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> crs_matrix_type;
61  typedef Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> block_crs_matrix_type;
62 
63 #if defined(HAVE_MUELU_AMGX)
64  std::string externalMG = "use external multigrid package";
65  if (inParamList.isParameter(externalMG) && inParamList.get<std::string>(externalMG) == "amgx"){
66  const RCP<crs_matrix_type> constCrsA = rcp_dynamic_cast<crs_matrix_type>(inA);
67  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.");
68  return rcp(new AMGXOperator<SC,LO,GO,NO>(constCrsA,inParamList));
69  }
70 #endif
71 
72  // Wrap A
73  RCP<Matrix> A;
74  RCP<block_crs_matrix_type> bcrsA = rcp_dynamic_cast<block_crs_matrix_type>(inA);
75  RCP<crs_matrix_type> crsA = rcp_dynamic_cast<crs_matrix_type>(inA);
76  if (crsA != Teuchos::null)
77  A = TpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(crsA);
78  else if (bcrsA != Teuchos::null) {
79  RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > temp = rcp(new Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>(bcrsA));
80  TEUCHOS_TEST_FOR_EXCEPTION(temp==Teuchos::null, Exceptions::RuntimeError, "CreateTpetraPreconditioner: cast from Tpetra::BlockCrsMatrix to Xpetra::TpetraBlockCrsMatrix failed.");
81  A = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(temp));
82  }
83  else {
84  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "CreateTpetraPreconditioner: only Tpetra CrsMatrix and BlockCrsMatrix types are supported.");
85  }
86 
87  Teuchos::ParameterList& userList = inParamList.sublist("user data");
88  if (userList.isParameter("Coordinates")) {
89  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> > coordinates = Teuchos::null;
90  try {
91  coordinates = TpetraMultiVector_To_XpetraMultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>(userList.get<RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > >("Coordinates"));
92  } catch(Teuchos::Exceptions::InvalidParameterType&) {
93  coordinates = userList.get<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > >("Coordinates");
94  }
95  userList.set<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> > >("Coordinates", coordinates);
96  }
97 
98  if (userList.isParameter("Nullspace")) {
99  RCP<MultiVector> nullspace = Teuchos::null;
100  try {
101  nullspace = TpetraMultiVector_To_XpetraMultiVector<SC,LO,GO,NO>(userList.get<RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >("Nullspace"));
102  } catch(Teuchos::Exceptions::InvalidParameterType&) {
103  nullspace = userList.get<RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >("Nullspace");
104  }
105  userList.set<RCP<MultiVector> >("Nullspace", nullspace);
106  }
107 
108  RCP<Hierarchy> H = MueLu::CreateXpetraPreconditioner<SC,LO,GO,NO>(A, inParamList);
109  return rcp(new TpetraOperator<SC,LO,GO,NO>(H));
110  }
111 
112 
122  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
123  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
124  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
125  const std::string& xmlFileName)
126  {
127  Teuchos::ParameterList paramList;
128  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getDomainMap()->getComm());
129  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList);
130  }
131 
132 
141  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
143  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA)
144  {
145  Teuchos::ParameterList paramList;
146  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList);
147  }
148 
149 
157  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
158  void ReuseTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
159  MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
160  typedef Scalar SC;
161  typedef LocalOrdinal LO;
162  typedef GlobalOrdinal GO;
163  typedef Node NO;
164 
165  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
166  typedef MueLu ::Hierarchy<SC,LO,GO,NO> Hierarchy;
167 
168  RCP<Hierarchy> H = Op.GetHierarchy();
169  RCP<Matrix> A = TpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(inA);
170 
171  MueLu::ReuseXpetraPreconditioner<SC,LO,GO,NO>(A, H);
172  }
173 
174 
175 #ifdef HAVE_MUELU_DEPRECATED_CODE
176 
186  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
187  MUELU_DEPRECATED
188  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
189  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
190  Teuchos::ParameterList& inParamList,
191  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>>& inCoords,
192  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace)
193  {
194  typedef Scalar SC;
195  typedef LocalOrdinal LO;
196  typedef GlobalOrdinal GO;
197  typedef Node NO;
198 
199  using Teuchos::ParameterList;
200 
201  // Here the assumption is that the nullspace and coordinates passed
202  // as multivectors will overwrite data on the parameter list.
203  // If you want the data on the parameter list to prevail, then call the
204  // specialization that only accept an operator and a parameterList
205  Teuchos::ParameterList& userList = inParamList.sublist("user data");
206  if (Teuchos::nonnull(inCoords)) {
207  userList.set<RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> > >("Coordinates", inCoords);
208  }
209  if (Teuchos::nonnull(inNullspace)) {
210  userList.set<RCP<Tpetra::MultiVector<SC,LO,GO,NO> > >("Nullspace", inNullspace);
211  }
212 
213  return CreateTpetraPreconditioner<SC,LO,GO,NO>(inA, inParamList);
214  }
215 
216 
226  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
227  MUELU_DEPRECATED
228  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
229  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
230  Teuchos::ParameterList& inParamList,
231  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> >& inCoords)
232  {
233  typedef Scalar SC;
234  typedef LocalOrdinal LO;
235  typedef GlobalOrdinal GO;
236  typedef Node NO;
237 
238  using Teuchos::ParameterList;
239 
240  // Here the assumption is that the nullspace and coordinates passed
241  // as multivectors will overwrite data on the parameter list.
242  // If you want the data on the parameter list to prevail, then call the
243  // specialization that only accept an operator and a parameterList
244  Teuchos::ParameterList& userList = inParamList.sublist("user data");
245  if (Teuchos::nonnull(inCoords)) {
246  userList.set<RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> > >("Coordinates", inCoords);
247  }
248 
249  return CreateTpetraPreconditioner<SC,LO,GO,NO>(inA, inParamList);
250  }
251 
252 
264  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
265  MUELU_DEPRECATED
266  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
267  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
268  const std::string& xmlFileName,
269  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>>& inCoords,
270  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace)
271  {
272  // Create the parameter list and updated is using the xml file
273  Teuchos::ParameterList paramList;
274  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getDomainMap()->getComm());
275 
276  // Use the parameter list interface
277  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
278  }
279 
280 
291  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
292  MUELU_DEPRECATED
293  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
294  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
295  const std::string& xmlFileName,
296  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>>& inCoords)
297  {
298  // Create the parameter list and updated is using the xml file
299  Teuchos::ParameterList paramList;
300  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getDomainMap()->getComm());
301 
302  // Use the parameter list interface
303  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, Teuchos::null);
304  }
305 
306 
317  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
318  MUELU_DEPRECATED
319  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
320  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
321  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>>& inCoords,
322  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace)
323  {
324  // Create the parameter list and put user data on it
325  Teuchos::ParameterList paramList;
326  Teuchos::ParameterList& userList = paramList.sublist("user data");
327  if (Teuchos::nonnull(inCoords)) {
328  userList.set<RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> > >("Coordinates", inCoords);
329  }
330  if (Teuchos::nonnull(inNullspace)) {
331  userList.set<RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >("Nullspace", inNullspace);
332  }
333 
334  // Use the parameter list interface
335  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
336  }
337 
338 
351  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
352  MUELU_DEPRECATED
353  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
354  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &inA,
355  Teuchos::ParameterList& inParamList,
356  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
357  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
358  {
359  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
360  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, inParamList, inCoords, inNullspace);
361  }
362 
363 
376  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
377  MUELU_DEPRECATED
378  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
379  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix <Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
380  const std::string& xmlFileName,
381  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
382  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
383  {
384  Teuchos::ParameterList paramList;
385  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(&paramList), *inA->getComm());
386 
387  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
388  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, paramList, inCoords, inNullspace);
389  }
390 
391 
403  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
404  MUELU_DEPRECATED
405  Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
406  CreateTpetraPreconditioner(const Teuchos::RCP<Tpetra::CrsMatrix <Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
407  const Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>>& inCoords = Teuchos::null,
408  const Teuchos::RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& inNullspace = Teuchos::null)
409  {
410  Teuchos::ParameterList paramList;
411  RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
412  return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, paramList, inCoords, inNullspace);
413  }
414 
415 #endif // HAVE_MUELU_DEPRECATED_CODE
416 
417 } //namespace
418 
419 #endif //ifndef MUELU_CREATE_TPETRA_PRECONDITIONER_HPP
420 
MueLu_AMGXOperator.hpp
MueLu_CreateXpetraPreconditioner.hpp
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
MueLu_Exceptions.hpp
MueLu_MasterList.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:50
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:49
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:51
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:52
MueLu::AMGXOperator
Adapter for AmgX library from Nvidia.
Definition: MueLu_AMGXOperator_decl.hpp:81
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:104
MueLu::TpetraOperator
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
Definition: MueLu_TpetraOperator_decl.hpp:65
MueLu::TpetraOperator::GetHierarchy
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetHierarchy() const
Direct access to the underlying MueLu::Hierarchy.
Definition: MueLu_TpetraOperator_def.hpp:137
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::CreateTpetraPreconditioner
Teuchos::RCP< MueLu::TpetraOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateTpetraPreconditioner(const Teuchos::RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &inA, Teuchos::ParameterList &inParamList)
Helper function to create a MueLu or AMGX preconditioner that can be used by Tpetra....
Definition: MueLu_CreateTpetraPreconditioner.hpp:47
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:158