46#ifndef MUELU_PROJECTORSMOOTHER_DEF_HPP
47#define MUELU_PROJECTORSMOOTHER_DEF_HPP
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_MultiVectorFactory.hpp>
56#include "MueLu_Utilities.hpp"
61 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
63 : coarseSolver_(coarseSolver)
66 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 coarseSolver_->DeclareInput(currentLevel);
77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 coarseSolver_->Setup(currentLevel);
84 this->GetOStream(
Warnings0) <<
"MueLu::ProjectorSmoother::Setup(): Setup() has already been called" << std::endl;
86 RCP<Matrix> A = Factory::Get< RCP<Matrix> > (currentLevel,
"A");
87 RCP<MultiVector> B = Factory::Get< RCP<MultiVector> >(currentLevel,
"Nullspace");
89 int m = B->getNumVectors();
92 Array<Scalar> br(m), bb(m);
93 RCP<MultiVector> R = MultiVectorFactory::Build(B->getMap(), m);
98 Array<size_t> selectedIndices;
99 for (
int i = 0; i < m; i++) {
100 Scalar rayleigh = br[i] / bb[i];
102 if (Teuchos::ScalarTraits<Scalar>::magnitude(rayleigh) < 1e-12)
103 selectedIndices.push_back(i);
105 this->GetOStream(
Runtime0) <<
"Coarse level orth indices: " << selectedIndices << std::endl;
107#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_XPETRA_TPETRA)
108#ifdef HAVE_MUELU_TPETRA_INST_INT_INT
113 RCP<Tpetra::MultiVector<SC,LO,GO,NO> > Borth = B_->subCopy(selectedIndices);
114 for (
int i = 0; i < selectedIndices.size(); i++) {
115 RCP<Tpetra::Vector<SC,LO,GO,NO> > Bi = Borth->getVectorNonConst(i);
118 typename Teuchos::ScalarTraits<Scalar>::magnitudeType norm;
119 for (
int k = 0; k < i; k++) {
120 RCP<const Tpetra::Vector<SC,LO,GO,NO> > Bk = Borth->getVector(k);
123 Bi->update(-dot, *Bk, Teuchos::ScalarTraits<Scalar>::one());
127 Bi->scale(Teuchos::ScalarTraits<Scalar>::one()/norm);
130 Borth_ = rcp(
static_cast<MultiVector*
>(
new TpetraMultiVector(Borth)));
132 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"Tpetra with GO=int not available. The code in ProjectorSmoother should be rewritten!");
139 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
141 coarseSolver_->Apply(X, B, InitialGuessIsZero);
143 int m = Borth_->getNumVectors();
144 int n = X.getNumVectors();
146 RCP<Xpetra::MultiVector<SC,LO,GO,NO> > X_ = Teuchos::rcpFromRef(X);
147 for (
int i = 0; i < n; i++) {
148 RCP<Xpetra::Vector<SC,LO,GO,NO> > Xi = X_->getVectorNonConst(i);
150 Array<Scalar> dot(1);
151 for (
int k = 0; k < m; k++) {
152 RCP<const Xpetra::Vector<SC,LO,GO,NO> > Bk = Borth_->getVector(k);
155 Xi->update(-dot[0], *Bk, Teuchos::ScalarTraits<Scalar>::one());
160 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
165 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
167 std::ostringstream out;
172 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
MueLu::DefaultScalar Scalar
virtual std::string description() const
Return a simple one-line description of this object.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Input(Level &level, const std::string &varName) const
Class that holds all level-specific information.
This class enables the elimination of the nullspace component of the solution through the use of proj...
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver. Solves the linear system AX=B using the constructed solver.
RCP< SmootherPrototype > Copy() const
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
virtual ~ProjectorSmoother()
Destructor.
void Setup(Level ¤tLevel)
Set up the direct solver.
std::string description() const
Return a simple one-line description of this object.
ProjectorSmoother(RCP< SmootherPrototype > coarseSolver)
Constructor.
void DeclareInput(Level ¤tLevel) const
Input.
bool IsSetup() const
Get the state of a smoother prototype.
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Runtime0
One-liner description of what is happening.