46 #ifndef MUELU_BLOCKEDRAPFACTORY_DEF_HPP 47 #define MUELU_BLOCKEDRAPFACTORY_DEF_HPP 49 #ifdef HAVE_MUELU_EXPERIMENTAL 51 #include <Xpetra_Matrix.hpp> 52 #include <Xpetra_MatrixMatrix.hpp> 53 #include <Xpetra_CrsMatrix.hpp> 54 #include <Xpetra_CrsMatrixWrap.hpp> 55 #include <Xpetra_MatrixFactory.hpp> 56 #include <Xpetra_BlockedCrsMatrix.hpp> 57 #include <Xpetra_Vector.hpp> 58 #include <Xpetra_VectorFactory.hpp> 64 #include "MueLu_PerfUtils.hpp" 66 #include "MueLu_Utilities.hpp" 70 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 : checkAc_(false), repairZeroDiagonals_(false)
75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 RCP<ParameterList> validParamList = rcp(
new ParameterList());
79 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 81 #undef SET_VALID_ENTRY 82 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A used during the prolongator smoothing process");
83 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Prolongator factory");
84 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Restrictor factory");
86 return validParamList;
89 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 if (pL.get<
bool>(
"transpose: use implicit") ==
false)
93 Input(coarseLevel,
"R");
95 Input(fineLevel,
"A");
96 Input(coarseLevel,
"P");
100 (*it)->CallDeclareInput(coarseLevel);
103 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
109 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
110 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P");
112 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);
113 RCP<BlockedCrsMatrix> bP = rcp_dynamic_cast<BlockedCrsMatrix>(P);
114 TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null() || bP.is_null(),
Exceptions::BadCast,
"Matrices R, A and P must be of type BlockedCrsMatrix.");
116 RCP<BlockedCrsMatrix> bAP;
117 RCP<BlockedCrsMatrix> bAc;
123 "Block matrix dimensions do not match: " 124 "A is " << bA->Rows() <<
"x" << bA->Cols() <<
125 "P is " << bP->Rows() <<
"x" << bP->Cols());
127 bAP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixMultiplyBlock(*bA,
false, *bP,
false,
GetOStream(
Statistics2),
true,
true);
135 const bool doTranspose =
true;
136 const bool doFillComplete =
true;
137 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
139 bAc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixMultiplyBlock(*bP, doTranspose, *bAP, !doTranspose,
GetOStream(
Statistics2), doFillComplete, doOptimizeStorage);
142 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
143 RCP<BlockedCrsMatrix> bR = rcp_dynamic_cast<BlockedCrsMatrix>(R);
144 TEUCHOS_TEST_FOR_EXCEPTION(bR.is_null(),
Exceptions::BadCast,
"Matrix R must be of type BlockedCrsMatrix.");
147 "Block matrix dimensions do not match: " 148 "R is " << bR->Rows() <<
"x" << bR->Cols() <<
149 "A is " << bA->Rows() <<
"x" << bA->Cols());
152 bAc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixMultiplyBlock(*bR, !doTranspose, *bAP, !doTranspose,
GetOStream(
Statistics2), doFillComplete, doOptimizeStorage);
172 Set<RCP <Matrix> >(coarseLevel,
"A", bAc);
179 RCP<const FactoryBase> fac = *it;
180 GetOStream(
Runtime0) <<
"BlockRAPFactory: call transfer factory: " << fac->description() << std::endl;
181 fac->CallBuild(coarseLevel);
184 coarseLevel.Release(*fac);
190 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
192 RCP<CrsMatrix> c00 = bAc->getMatrix(0, 0);
193 RCP<CrsMatrix> Aout = CrsMatrixFactory::Build(c00->getRowMap(), c00->getGlobalMaxNumRowEntries(), Xpetra::StaticProfile);
195 RCP<Vector> diagVec = VectorFactory::Build(c00->getRowMap());
196 c00->getLocalDiagCopy(*diagVec);
197 ArrayRCP<SC> diagVal = diagVec->getDataNonConst(0);
200 for(
size_t row=0; row<c00->getNodeNumRows(); row++) {
202 GO grid = c00->getRowMap()->getGlobalElement(row);
204 ArrayView<const LO> indices;
205 ArrayView<const SC> vals;
206 c00->getLocalRowView(row, indices, vals);
209 ArrayRCP<GO> indout(indices.size(), Teuchos::OrdinalTraits<GO>::zero());
210 ArrayRCP<SC> valout(indices.size(), Teuchos::ScalarTraits<SC>::zero());
213 for (
size_t i = 0; i < as<size_t>(indices.size()); i++) {
214 GO gcid = c00->getColMap()->getGlobalElement(indices[i]);
216 valout [i] = vals[i];
219 Aout->insertGlobalValues(grid, indout.view(0, indout.size()), valout.view(0, valout.size()));
220 if (diagVal[row] == 0.0 && repairZeroDiagonals) {
222 Aout->insertGlobalValues(grid, Teuchos::tuple<GO>(grid), Teuchos::tuple<SC>(1.0));
226 Aout->fillComplete(c00->getDomainMap(), c00->getRangeMap());
228 bAc->setMatrix(0, 0, Aout);
231 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
234 TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null,
Exceptions::BadCast,
235 "Transfer factory is not derived from TwoLevelFactoryBase. This is very strange. (Note: you can remove this exception if there's a good reason for)");
241 #define MUELU_BLOCKEDRAPFACTORY_SHORT 243 #endif // MUELU_BLOCKEDRAPFACTORY_DEF_HPP Exception indicating invalid cast attempted.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
virtual const Teuchos::ParameterList & GetParameterList() const
Print even more statistics.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static void CheckMainDiagonal(RCP< BlockedCrsMatrix > &bAc, bool repairZeroDiagonals=false)
#define SET_VALID_ENTRY(name)
std::vector< RCP< const FactoryBase > > transferFacts_
list of user-defined transfer Factories
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
void Input(Level &level, const std::string &varName) const