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);
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)