46#ifndef MUELU_CGSOLVER_DEF_HPP
47#define MUELU_CGSOLVER_DEF_HPP
49#include <Xpetra_MatrixFactory.hpp>
50#include <Xpetra_MatrixMatrix.hpp>
52#include "MueLu_Utilities.hpp"
53#include "MueLu_Constraint.hpp"
57#include "MueLu_CGSolver.hpp"
63 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
68 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 finalP = MatrixFactory2::BuildCopy(rcpFromRef(P0));
79 RCP<const Matrix> A = rcpFromRef(Aref);
81 bool useTpetra = (A->getRowMap()->lib() == Xpetra::UseTpetra);
83 Teuchos::FancyOStream& mmfancy = this->GetOStream(
Statistics2);
85 SC one = Teuchos::ScalarTraits<SC>::one();
87 RCP<Matrix> X, P, R, Z, AP;
88 RCP<Matrix> newX, tmpAP;
89#ifndef TWO_ARG_MATRIX_ADD
90 RCP<Matrix> newR, newP;
93 SC oldRZ, newRZ, alpha, beta, app;
96 RCP<CrsMatrix> T_ = CrsMatrixFactory::Build(C.
GetPattern());
97 T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
98 RCP<Matrix> T = rcp(
new CrsMatrixWrap(T_));
101 X = rcp_const_cast<Matrix>(rcpFromRef(P0));
103 tmpAP = MatrixMatrix::Multiply(*A,
false, *X,
false, mmfancy,
true,
true);
107 R = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(T);
112 Z = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(R);
116 P = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(Z);
120 for (
size_t i = 0; i < nIts_; i++) {
122 if (i == 0 || useTpetra) {
127 tmpAP = MatrixMatrix::Multiply(*A,
false, *P,
false, mmfancy,
true,
true);
130 tmpAP = MatrixMatrix::Multiply(*A,
false, *P,
false, tmpAP, mmfancy,
true,
true);
136 if (Teuchos::ScalarTraits<SC>::magnitude(app) < Teuchos::ScalarTraits<SC>::sfmin()) {
141 X = MatrixFactory2::BuildCopy(rcpFromRef(P0));
147 this->GetOStream(
Runtime1,1) <<
"alpha = " << alpha << std::endl;
150#ifndef TWO_ARG_MATRIX_ADD
151 newX = Teuchos::null;
152 MatrixMatrix::TwoMatrixAdd(*P,
false, alpha, *X,
false, one, newX, mmfancy);
153 newX->fillComplete(P0.getDomainMap(), P0.getRangeMap());
156 MatrixMatrix::TwoMatrixAdd(*P,
false, alpha, *X, one);
163#ifndef TWO_ARG_MATRIX_ADD
164 newR = Teuchos::null;
165 MatrixMatrix::TwoMatrixAdd(*AP,
false, -alpha, *R,
false, one, newR, mmfancy);
166 newR->fillComplete(P0.getDomainMap(), P0.getRangeMap());
169 MatrixMatrix::TwoMatrixAdd(*AP,
false, -alpha, *R, one);
173 Z = MatrixFactory2::BuildCopy(R);
178 beta = newRZ / oldRZ;
181#ifndef TWO_ARG_MATRIX_ADD
182 newP = Teuchos::null;
183 MatrixMatrix::TwoMatrixAdd(*P,
false, beta, *Z,
false, one, newP, mmfancy);
184 newP->fillComplete(P0.getDomainMap(), P0.getRangeMap());
187 MatrixMatrix::TwoMatrixAdd(*Z,
false, one, *P, beta);
void Iterate(const Matrix &A, const Constraint &C, const Matrix &P0, RCP< Matrix > &P) const
Iterate.
Constraint space information for the potential prolongator.
RCP< const CrsGraph > GetPattern() const
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Runtime1
Description of what is happening (more verbose)