101 const ParameterList& pL = GetParameterList();
103 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
104 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P");
107 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);
108 RCP<BlockedCrsMatrix> bP = rcp_dynamic_cast<BlockedCrsMatrix>(P);
109 TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null() || bP.is_null(),
Exceptions::BadCast,
"Matrices A and P must be of type BlockedCrsMatrix.");
112 RCP<BlockedCrsMatrix> bAP;
113 RCP<BlockedCrsMatrix> bAc;
119 "Block matrix dimensions do not match: "
120 "A is " << bA->Rows() <<
"x" << bA->Cols() <<
121 "P is " << bP->Rows() <<
"x" << bP->Cols());
123 bAP = MatrixMatrix::TwoMatrixMultiplyBlock(*bA,
false, *bP,
false, GetOStream(
Statistics2),
true,
true);
129 bool doOptimizeStorage = !checkAc_;
131 const bool doTranspose =
true;
132 const bool doFillComplete =
true;
133 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
135 bAc = MatrixMatrix::TwoMatrixMultiplyBlock(*bP, doTranspose, *bAP, !doTranspose, GetOStream(
Statistics2), doFillComplete, doOptimizeStorage);
138 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
139 RCP<BlockedCrsMatrix> bR = rcp_dynamic_cast<BlockedCrsMatrix>(R);
140 TEUCHOS_TEST_FOR_EXCEPTION(bR.is_null(),
Exceptions::BadCast,
"Matrix R must be of type BlockedCrsMatrix.");
143 "Block matrix dimensions do not match: "
144 "R is " << bR->Rows() <<
"x" << bR->Cols() <<
145 "A is " << bA->Rows() <<
"x" << bA->Cols());
148 bAc = MatrixMatrix::TwoMatrixMultiplyBlock(*bR, !doTranspose, *bAP, !doTranspose, GetOStream(
Statistics2), doFillComplete, doOptimizeStorage);
153 CheckMainDiagonal(bAc);
157 Set<RCP <Matrix> >(coarseLevel,
"A", bAc);
159 if (transferFacts_.begin() != transferFacts_.end()) {
163 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
164 RCP<const FactoryBase> fac = *it;
166 GetOStream(
Runtime0) <<
"BlockRAPFactory: call transfer factory: " << fac->description() << std::endl;
168 fac->CallBuild(coarseLevel);
180 RCP<Matrix> c00 = bAc->getMatrix(0, 0);
181 RCP<Matrix> Aout = MatrixFactory::Build(c00->getRowMap(), c00->getGlobalMaxNumRowEntries());
183 RCP<Vector> diagVec = VectorFactory::Build(c00->getRowMap());
184 c00->getLocalDiagCopy(*diagVec);
185 ArrayRCP<SC> diagVal = diagVec->getDataNonConst(0);
188 for (
size_t row = 0; row < c00->getLocalNumRows(); row++) {
190 GO grid = c00->getRowMap()->getGlobalElement(row);
192 ArrayView<const LO> indices;
193 ArrayView<const SC> vals;
194 c00->getLocalRowView(row, indices, vals);
197 ArrayRCP<GO> indout(indices.size(), Teuchos::OrdinalTraits<GO>::zero());
198 ArrayRCP<SC> valout(indices.size(), Teuchos::ScalarTraits<SC>::zero());
201 for (
size_t i = 0; i < as<size_t>(indices.size()); i++) {
202 GO gcid = c00->getColMap()->getGlobalElement(indices[i]);
204 valout [i] = vals[i];
207 Aout->insertGlobalValues(grid, indout.view(0, indout.size()), valout.view(0, valout.size()));
208 if (diagVal[row] == Teuchos::ScalarTraits<SC>::zero() && repairZeroDiagonals) {
210 Aout->insertGlobalValues(grid, Teuchos::tuple<GO>(grid), Teuchos::tuple<SC>(1.0));
214 Aout->fillComplete(c00->getDomainMap(), c00->getRangeMap());
216 bAc->setMatrix(0, 0, Aout);