46#ifndef MUELU_BLOCKEDRAPFACTORY_DEF_HPP
47#define MUELU_BLOCKEDRAPFACTORY_DEF_HPP
49#include <Xpetra_BlockedCrsMatrix.hpp>
50#include <Xpetra_MatrixFactory.hpp>
51#include <Xpetra_Matrix.hpp>
52#include <Xpetra_MatrixMatrix.hpp>
58#include "MueLu_PerfUtils.hpp"
60#include "MueLu_Utilities.hpp"
64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 : checkAc_(false), repairZeroDiagonals_(false)
69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 RCP<ParameterList> validParamList = rcp(
new ParameterList());
73#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
76 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A used during the prolongator smoothing process");
77 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Prolongator factory");
78 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Restrictor factory");
80 return validParamList;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 const Teuchos::ParameterList& pL = GetParameterList();
86 if (pL.get<
bool>(
"transpose: use implicit") ==
false)
87 Input(coarseLevel,
"R");
89 Input(fineLevel,
"A");
90 Input(coarseLevel,
"P");
93 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
94 (*it)->CallDeclareInput(coarseLevel);
97 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
178 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
219 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
222 TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null,
Exceptions::BadCast,
223 "Transfer factory is not derived from TwoLevelFactoryBase. This is very strange. "
224 "(Note: you can remove this exception if there's a good reason for)");
225 transferFacts_.push_back(factory);
230#define MUELU_BLOCKEDRAPFACTORY_SHORT
#define SET_VALID_ENTRY(name)
static void CheckMainDiagonal(RCP< BlockedCrsMatrix > &bAc, bool repairZeroDiagonals=false)
void Build(Level &fineLevel, Level &coarseLevel) const override
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const override
Input.
RCP< const ParameterList > GetValidParameterList() const override
Return a const parameter list of valid parameters that setParameterList() will accept.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
Exception indicating invalid cast attempted.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.
@ Runtime0
One-liner description of what is happening.