1 #ifndef MUELU_CREATE_TPETRA_PRECONDITIONER_HPP 2 #define MUELU_CREATE_TPETRA_PRECONDITIONER_HPP 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> 19 #include <MueLu_Hierarchy.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> 28 #if defined(HAVE_MUELU_EXPERIMENTAL) and defined(HAVE_MUELU_AMGX) 29 #include <MueLu_AMGXOperator.hpp> 31 #include "cuda_runtime.h" 46 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
47 Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
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)
54 typedef LocalOrdinal LO;
55 typedef GlobalOrdinal GO;
58 using Teuchos::ParameterList;
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;
67 bool hasParamList = inParamList.numParams();
69 RCP<HierarchyManager> mueLuFactory;
70 ParameterList paramList = inParamList;
71 RCP<const crs_matrix_type> constCrsA;
72 RCP<crs_matrix_type> crsA;
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.");
82 std::string syntaxStr =
"parameterlist: syntax";
83 if (hasParamList && paramList.isParameter(syntaxStr) && paramList.get<std::string>(syntaxStr) ==
"ml") {
84 paramList.remove(syntaxStr);
91 RCP<Hierarchy> H = mueLuFactory->CreateHierarchy();
92 H->setlib(Xpetra::UseTpetra);
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));
106 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"CreateTpetraPreconditioner: only Tpetra CrsMatrix and BlockCrsMatrix types are supported.");
108 H->GetLevel(0)->Set(
"A", A);
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);
117 RCP<MultiVector> nullspace;
118 if (inNullspace != Teuchos::null) {
119 nullspace = TpetraMultiVector_To_XpetraMultiVector<SC,LO,GO,NO>(inNullspace);
122 int nPDE = MasterList::getDefault<int>(
"number of equations");
123 if (paramList.isSublist(
"Matrix")) {
125 const Teuchos::ParameterList& operatorList = paramList.sublist(
"Matrix");
126 if (operatorList.isParameter(
"PDE equations"))
127 nPDE = operatorList.get<
int>(
"PDE equations");
129 }
else if (paramList.isParameter(
"number of equations")) {
131 nPDE = paramList.get<
int>(
"number of equations");
134 nullspace = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(A->getDomainMap(), nPDE);
136 nullspace->putScalar(Teuchos::ScalarTraits<SC>::one());
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();
144 if ((GID-i) % nPDE == 0)
145 nsData[j] = Teuchos::ScalarTraits<SC>::one();
150 H->GetLevel(0)->Set(
"Nullspace", nullspace);
153 Teuchos::ParameterList nonSerialList,dummyList;
157 mueLuFactory->SetupHierarchy(*H);
174 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
176 Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
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)
182 RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
183 return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, inParamList, inCoords, inNullspace);
198 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
200 Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
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)
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);
223 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
225 Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
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)
231 Teuchos::ParameterList paramList;
232 Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(¶mList), *inA->getComm());
234 RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> opMat(inA);
235 return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(opMat, paramList, inCoords, inNullspace);
249 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
250 Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
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)
255 Teuchos::ParameterList paramList;
256 return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
271 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
272 Teuchos::RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
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)
278 Teuchos::ParameterList paramList;
279 Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<Teuchos::ParameterList>(¶mList), *inA->getDomainMap()->getComm());
280 return CreateTpetraPreconditioner<Scalar, LocalOrdinal, GlobalOrdinal, Node>(inA, paramList, inCoords, inNullspace);
291 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
295 typedef LocalOrdinal LO;
296 typedef GlobalOrdinal GO;
299 typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
300 typedef Xpetra::Operator<SC,LO,GO,NO> Operator;
306 "MueLu::ReuseTpetraPreconditioner: Hierarchy has no levels in it");
308 "MueLu::ReuseTpetraPreconditioner: Hierarchy has no fine level operator");
309 RCP<Level> level0 = H->GetLevel(0);
311 RCP<Operator> O0 = level0->Get<RCP<Operator> >(
"A");
312 RCP<Matrix> A0 = Teuchos::rcp_dynamic_cast<Matrix>(O0);
314 RCP<Matrix> A = TpetraCrs_To_XpetraMatrix<SC,LO,GO,NO>(inA);
319 A->SetFixedBlockSize(A0->GetFixedBlockSize());
328 #endif //ifndef MUELU_CREATE_TPETRA_PRECONDITIONER_HPP 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.