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>
73 finalP = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(rcpFromRef(P0));
79 RCP<const Matrix> A = rcpFromRef(Aref);
81 RCP<Matrix> X, P, R, Z, AP;
82 RCP<Matrix> newX, tmpAP;
83 #ifndef TWO_ARG_MATRIX_ADD 84 RCP<Matrix> newR, newP;
87 SC oldRZ, newRZ, alpha, beta, app;
89 bool useTpetra = (A->getRowMap()->lib() == Xpetra::UseTpetra);
94 RCP<CrsMatrix> T_ = Xpetra::CrsMatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(C.
GetPattern());
95 T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
96 RCP<Matrix> T = rcp(
new CrsMatrixWrap(T_));
98 SC one = Teuchos::ScalarTraits<SC>::one();
103 X = rcp_const_cast<Matrix>(rcpFromRef(P0));
105 bool doFillComplete =
true;
107 bool optimizeStorage =
true;
109 tmpAP = MatrixMatrix::Multiply(*A,
false, *X,
false, mmfancy, doFillComplete, optimizeStorage);
113 R = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(T);
114 #ifdef HAVE_MUELU_TPETRA 115 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT 122 this->
GetOStream(
Warnings0) <<
"WARNING: MueLu_CGSolver: calling Xpetra::CrsMatrix::resumeFill instead of Tpetra::CrsMatrix::resumeFill. The results should be verified in this case." << std::endl;
128 R->fillComplete(R->getDomainMap(), R->getRangeMap());
131 Z = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(R);
135 P = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(Z);
139 for (
size_t k = 0; k <
nIts_; k++) {
141 if (k == 0 || useTpetra) {
146 tmpAP = MatrixMatrix::Multiply(*A,
false, *P,
false, mmfancy, doFillComplete, optimizeStorage);
149 tmpAP = MatrixMatrix::Multiply(*A,
false, *P,
false, tmpAP, mmfancy, doFillComplete, optimizeStorage);
155 if (Teuchos::ScalarTraits<SC>::magnitude(app) < Teuchos::ScalarTraits<SC>::sfmin()) {
160 X = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(rcpFromRef(P0));
169 #ifndef TWO_ARG_MATRIX_ADD 170 newX = Teuchos::null;
171 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*P,
false, alpha, *X,
false, Teuchos::ScalarTraits<Scalar>::one(), newX, mmfancy);
172 newX->fillComplete(P0.getDomainMap(), P0.getRangeMap());
175 MatrixMatrix::TwoMatrixAdd(*P,
false, alpha, *X, one);
182 #ifndef TWO_ARG_MATRIX_ADD 183 newR = Teuchos::null;
184 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*AP,
false, -alpha, *R,
false, Teuchos::ScalarTraits<Scalar>::one(), newR, mmfancy);
185 newR->fillComplete(P0.getDomainMap(), P0.getRangeMap());
188 MatrixMatrix::TwoMatrixAdd(*AP,
false, -alpha, *R, one);
192 Z = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(R);
197 beta = newRZ / oldRZ;
200 #ifndef TWO_ARG_MATRIX_ADD 201 newP = Teuchos::null;
202 MatrixMatrix::TwoMatrixAdd(*P,
false, beta, *Z,
false, Teuchos::ScalarTraits<Scalar>::one(), newP, mmfancy);
203 newP->fillComplete(P0.getDomainMap(), P0.getRangeMap());
206 MatrixMatrix::TwoMatrixAdd(*Z,
false, one, *P, beta);
217 #endif //ifndef MUELU_CGSOLVER_DECL_HPP static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Important warning messages (one line)
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.
RCP< const CrsGraph > GetPattern() const
Print even more statistics.
void Iterate(const Matrix &A, const Constraint &C, const Matrix &P0, RCP< Matrix > &P) const
Iterate.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
size_t nIts_
Number of performed iterations.
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)