46 #ifndef MUELU_PROJECTORSMOOTHER_DEF_HPP 47 #define MUELU_PROJECTORSMOOTHER_DEF_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>
180 #endif // MUELU_PROJECTORSMOOTHER_DEF_HPP Important warning messages (one line)
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.
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.
bool IsSetup() const
Get the state of a smoother prototype.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
virtual ~ProjectorSmoother()
Destructor.
std::string description() const
Return a simple one-line description of this object.
void DeclareInput(Level ¤tLevel) const
Input.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
Class that holds all level-specific information.
RCP< SmootherPrototype > Copy() const
void Setup(Level ¤tLevel)
Set up the direct solver.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
This class enables the elimination of the nullspace component of the solution through the use of proj...
void Input(Level &level, const std::string &varName) const
Exception throws to report errors in the internal logical of the program.
ProjectorSmoother(RCP< SmootherPrototype > coarseSolver)
Constructor.
virtual std::string description() const
Return a simple one-line description of this object.