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>
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;
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,
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 
176 } //namespace
177 
178 #endif //ifndef MUELU_CREATE_TPETRA_PRECONDITIONER_HPP
179 
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
MueLu::DefaultLocalOrdinal LocalOrdinal
Namespace for MueLu classes and methods.
MueLu::DefaultNode Node
RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetHierarchy() const
Direct access to the underlying MueLu::Hierarchy.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
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.
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.Given a Tpetra::Operator, this function returns a constructed MueLu preconditioner.
Adapter for AmgX library from Nvidia.