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!");
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());
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.