46 #ifndef MUELU_CGSOLVER_DEF_HPP 47 #define MUELU_CGSOLVER_DEF_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);
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);
111 R->fillComplete(R->getDomainMap(), R->getRangeMap());
122 for (
size_t i = 0; i < nIts_; i++) {
124 if (i == 0 || useTpetra) {
129 tmpAP = MatrixMatrix::Multiply(*A,
false, *P,
false, mmfancy,
true,
true);
132 tmpAP = MatrixMatrix::Multiply(*A,
false, *P,
false, tmpAP, mmfancy,
true,
true);
138 if (Teuchos::ScalarTraits<SC>::magnitude(app) < Teuchos::ScalarTraits<SC>::sfmin()) {
143 X = MatrixFactory2::BuildCopy(rcpFromRef(P0));
149 this->GetOStream(
Runtime1,1) <<
"alpha = " << alpha << std::endl;
152 #ifndef TWO_ARG_MATRIX_ADD 153 newX = Teuchos::null;
154 MatrixMatrix::TwoMatrixAdd(*P,
false, alpha, *X,
false, one, newX, mmfancy);
155 newX->fillComplete(P0.getDomainMap(), P0.getRangeMap());
158 MatrixMatrix::TwoMatrixAdd(*P,
false, alpha, *X, one);
165 #ifndef TWO_ARG_MATRIX_ADD 166 newR = Teuchos::null;
167 MatrixMatrix::TwoMatrixAdd(*AP,
false, -alpha, *R,
false, one, newR, mmfancy);
168 newR->fillComplete(P0.getDomainMap(), P0.getRangeMap());
171 MatrixMatrix::TwoMatrixAdd(*AP,
false, -alpha, *R, one);
175 Z = MatrixFactory2::BuildCopy(R);
180 beta = newRZ / oldRZ;
183 #ifndef TWO_ARG_MATRIX_ADD 184 newP = Teuchos::null;
185 MatrixMatrix::TwoMatrixAdd(*P,
false, beta, *Z,
false, one, newP, mmfancy);
186 newP->fillComplete(P0.getDomainMap(), P0.getRangeMap());
189 MatrixMatrix::TwoMatrixAdd(*Z,
false, one, *P, beta);
200 #endif //ifndef MUELU_CGSOLVER_DECL_HPP static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
RCP< const CrsGraph > GetPattern() const
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)
Constraint space information for the potential prolongator.
Namespace for MueLu classes and methods.
Print even more statistics.
void Iterate(const Matrix &A, const Constraint &C, const Matrix &P0, RCP< Matrix > &P) const
Iterate.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > BuildCopy(const RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > A)
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
Description of what is happening (more verbose)
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)