46 #ifndef MUELU_STRATIMIKOSSMOOTHER_DEF_HPP 47 #define MUELU_STRATIMIKOSSMOOTHER_DEF_HPP 51 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_THYRA) 53 #include <Teuchos_ParameterList.hpp> 55 #include <Xpetra_CrsMatrix.hpp> 56 #include <Xpetra_CrsMatrixWrap.hpp> 57 #include <Xpetra_Matrix.hpp> 62 #include "MueLu_Utilities.hpp" 68 #include <Stratimikos_DefaultLinearSolverBuilder.hpp> 69 #include "Teuchos_AbstractFactoryStd.hpp" 70 #if defined(HAVE_MUELU_THYRATPETRAADAPTERS) && defined(HAVE_MUELU_IFPACK2) 71 #include <Thyra_Ifpack2PreconditionerFactory.hpp> 73 #include <Teuchos_ParameterList.hpp> 74 #include <unordered_map> 79 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::StratimikosSmoother(
const std::string type,
const Teuchos::ParameterList& paramList)
83 std::transform(type_.begin(), type_.end(), type_.begin(), ::toupper);
84 ParameterList& pList =
const_cast<ParameterList&
>(paramList);
86 if (pList.isParameter(
"smoother: recurMgOnFilteredA")) {
87 recurMgOnFilteredA_ =
true;
88 pList.remove(
"smoother: recurMgOnFilteredA");
90 bool isSupported = type_ ==
"STRATIMIKOS";
91 this->declareConstructionOutcome(!isSupported,
"Stratimikos does not provide the smoother '" + type_ +
"'.");
93 SetParameterList(paramList);
96 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetParameterList(
const Teuchos::ParameterList& paramList) {
98 Factory::SetParameterList(paramList);
101 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
102 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel)
const {
103 this->Input(currentLevel,
"A");
106 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
107 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Setup(Level& currentLevel) {
108 FactoryMonitor m(*
this,
"Setup Smoother", currentLevel);
110 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
111 SetupStratimikos(currentLevel);
112 SmootherPrototype::IsSetup(
true);
113 this->GetOStream(
Statistics1) << description() << std::endl;
116 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetupStratimikos(Level& currentLevel) {
119 RCP<const Thyra::LinearOpBase<Scalar> > thyraA;
120 if (recurMgOnFilteredA_) {
121 RCP<Matrix> filteredA;
122 ExperimentalDropVertConnections(filteredA, currentLevel);
123 thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(filteredA)->getCrsMatrix());
125 else thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(A_)->getCrsMatrix());
128 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
129 if (recurMgOnFilteredA_) {
131 Stratimikos::enableMueLu<LocalOrdinal,GlobalOrdinal,Node>(linearSolverBuilder);
133 TEUCHOS_TEST_FOR_EXCEPTION(
true , Exceptions::RuntimeError,
"MueLu::StratimikosSmoother:: must compile with MUELU_RECURMG defined. Unfortunately, cmake does not always produce a proper link.txt file (which sometimes requires libmuelu.a before and after libmuelu-interface.a). After configuring, run script muelu/utils/misc/patchLinkForRecurMG to change link.txt files manually. If you want to create test example, add -DMUELU_RECURMG=ON to cmake arguments.");
137 typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
138 #if defined(HAVE_MUELU_THYRATPETRAADAPTERS) && defined(HAVE_MUELU_IFPACK2) 139 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Impl;
140 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(),
"Ifpack2");
142 linearSolverBuilder.setParameterList(rcpFromRef(const_cast<ParameterList&>(this->GetParameterList())));
146 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
147 solver_ = Thyra::linearOpWithSolve(*solverFactory, thyraA);
148 #ifdef dumpOutRecurMGDebug 150 sprintf(mystring,
"for i in A_[0123456789].m P_[0123456789].m; do T=Xecho $i | sed Xs/.m$/%d.m/XX; mv $i $T; done", (
int) currentLevel.GetLevelID()); fflush(stdout);
151 mystring[50]=
'`'; mystring[65]=
'"'; mystring[76]=
'"'; mystring[77]=
'`';
156 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
157 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::ExperimentalDropVertConnections(RCP<Matrix> & filteredA, Level& currentLevel) {
194 using TST =
typename Teuchos::ScalarTraits<SC>;
196 bool sumDropped =
false;
198 LO dofsPerNode = A_->GetFixedBlockSize();
201 RCP<ParameterList> fillCompleteParams(
new ParameterList);
202 fillCompleteParams->set(
"No Nonlocal Changes",
true);
203 filteredA = MatrixFactory::Build(A_->getCrsGraph());
204 filteredA->resumeFill();
206 ArrayView<const LocalOrdinal> inds;
207 ArrayView<const Scalar> valsA;
208 ArrayView<Scalar> vals;
209 Teuchos::ArrayRCP<LocalOrdinal> TVertLineId = Factory::Get< Teuchos::ArrayRCP<LocalOrdinal> > (currentLevel,
"LineDetection_VertLineIds");
210 Teuchos::ArrayRCP<LocalOrdinal> TLayerId = Factory::Get< Teuchos::ArrayRCP<LocalOrdinal> > (currentLevel,
"LineDetection_Layers");
213 TEUCHOS_TEST_FOR_EXCEPTION( (LayerId == NULL) || (VertLineId == NULL), Exceptions::RuntimeError,
"MueLu::StratimikosSmoother:: no line information found on this level. Cannot use recurMgOnFilteredA on this level.");
215 Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
216 for (
size_t i = 0; i < A_->getRowMap()->getNodeNumElements(); i++) {
217 A_->getLocalRowView(i, inds, valsA);
218 size_t nnz = inds.size();
219 ArrayView<const Scalar> vals1;
220 filteredA->getLocalRowView(i, inds, vals1);
221 vals = ArrayView<Scalar>(
const_cast<Scalar*
>(vals1.getRawPtr()), nnz);
222 memcpy(vals.getRawPtr(), valsA.getRawPtr(), nnz*
sizeof(
Scalar));
223 size_t inode, jdof, jnode, jdof_offset;
224 inode = i/dofsPerNode;
226 std::unordered_map<LocalOrdinal, LocalOrdinal> umap;
230 for (
size_t j = 0; j < nnz; j++) {
232 jnode = jdof/dofsPerNode;
233 jdof_offset = jdof - jnode*dofsPerNode;
234 if ( LayerId[jnode] == LayerId[inode] ) umap[dofsPerNode*VertLineId[jnode]+jdof_offset] = j;
239 for (
size_t j = 0; j < nnz; j++) {
241 jnode = jdof/dofsPerNode;
242 jdof_offset = jdof - jnode*dofsPerNode;
243 if ( LayerId[jnode] != LayerId[inode] ) {
245 if (umap.find(dofsPerNode*VertLineId[jnode + jdof_offset]) != umap.end())
246 vals[umap[dofsPerNode*VertLineId[jnode + jdof_offset]]] += vals[j];
252 filteredA->fillComplete(fillCompleteParams);
256 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
257 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X,
const MultiVector& B,
bool InitialGuessIsZero)
const {
258 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() ==
false, Exceptions::RuntimeError,
"MueLu::StratimikosSmoother::Apply(): Setup() has not been called");
261 if (InitialGuessIsZero) {
263 RCP< Thyra::MultiVectorBase<Scalar> > thyraX = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpFromRef(X)));
264 RCP<const Thyra::MultiVectorBase<Scalar> > thyraB = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpFromRef(B));
265 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraB, thyraX.ptr());
266 RCP<MultiVector> thyXpX = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraX, X.getMap()->getComm());
270 typedef Teuchos::ScalarTraits<Scalar> TST;
271 RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
273 RCP<MultiVector> Cor = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X.getMap(),X.getNumVectors(),
true);
274 RCP< Thyra::MultiVectorBase<Scalar> > thyraCor = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(Cor));
275 RCP<const Thyra::MultiVectorBase<Scalar> > thyraRes = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(Residual);
276 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraRes, thyraCor.ptr());
277 RCP<MultiVector> thyXpCor = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraCor, X.getMap()->getComm());
278 X.update(TST::one(), *thyXpCor, TST::one());
284 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
285 RCP<MueLu::SmootherPrototype<double, LocalOrdinal, GlobalOrdinal, Node> > StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Copy()
const {
286 RCP<StratimikosSmoother> smoother = rcp(
new StratimikosSmoother(*
this) );
287 smoother->SetParameterList(this->GetParameterList());
291 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
292 std::string StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::description()
const {
293 std::ostringstream out;
294 if (SmootherPrototype::IsSetup()) {
295 out << solver_->description();
297 out <<
"STRATIMIKOS {type = " << type_ <<
"}";
302 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
303 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out,
const VerbLevel verbLevel)
const {
307 out0 <<
"Parameter list: " << std::endl;
308 Teuchos::OSTab tab2(out);
309 out << this->GetParameterList();
313 if (solver_ != Teuchos::null) {
314 Teuchos::OSTab tab2(out);
315 out << *solver_ << std::endl;
318 if (verbLevel &
Debug) {
319 out0 <<
"IsSetup: " <<
Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
321 <<
"RCP<solver_>: " << solver_ << std::endl;
325 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
326 size_t StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::getNodeSmootherComplexity()
const {
327 return Teuchos::OrdinalTraits<size_t>::invalid();
333 #endif // HAVE_MUELU_STRATIMIKOS 334 #endif // MUELU_STRATIMIKOSSMOOTHER_DEF_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Print external lib objects.
Print additional debugging information.
Namespace for MueLu classes and methods.
MueLu::DefaultScalar Scalar
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters (more parameters, more verbose)