46 #ifndef MUELU_RAPFACTORY_DEF_HPP 47 #define MUELU_RAPFACTORY_DEF_HPP 52 #include <Xpetra_Matrix.hpp> 53 #include <Xpetra_MatrixFactory.hpp> 54 #include <Xpetra_MatrixMatrix.hpp> 55 #include <Xpetra_Vector.hpp> 56 #include <Xpetra_VectorFactory.hpp> 62 #include "MueLu_PerfUtils.hpp" 68 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
70 : hasDeclaredInput_(false) { }
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 RCP<ParameterList> validParamList = rcp(
new ParameterList());
76 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 78 #undef SET_VALID_ENTRY 79 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A used during the prolongator smoothing process");
80 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Prolongator factory");
81 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Restrictor factory");
82 validParamList->set< RCP<const FactoryBase> >(
"AP Pattern", null,
"AP pattern factory");
83 validParamList->set< RCP<const FactoryBase> >(
"RAP Pattern", null,
"RAP pattern factory");
84 validParamList->set<
bool > (
"Keep AP Pattern",
false,
"Keep the AP pattern (for reuse)");
85 validParamList->set<
bool > (
"Keep RAP Pattern",
false,
"Keep the RAP pattern (for reuse)");
87 validParamList->set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros (default = false).");
88 validParamList->set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal (default = false).");
90 return validParamList;
93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 if (pL.get<
bool>(
"transpose: use implicit") ==
false)
97 Input(coarseLevel,
"R");
99 Input(fineLevel,
"A");
100 Input(coarseLevel,
"P");
104 (*it)->CallDeclareInput(coarseLevel);
109 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 std::ostringstream levelstr;
117 "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
121 if (pL.get<
bool>(
"Keep AP Pattern"))
122 coarseLevel.
Keep(
"AP Pattern",
this);
123 if (pL.get<
bool>(
"Keep RAP Pattern"))
124 coarseLevel.
Keep(
"RAP Pattern",
this);
126 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
127 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P"), AP, Ac;
130 if (coarseLevel.IsAvailable(
"AP Pattern",
this)) {
133 AP = Get< RCP<Matrix> >(coarseLevel,
"AP Pattern");
139 AP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *P,
false, AP,
GetOStream(
Statistics2),
true,
true,std::string(
"MueLu::A*P-")+levelstr.str());
141 if (pL.get<
bool>(
"Keep AP Pattern"))
142 Set(coarseLevel,
"AP Pattern", AP);
145 if (coarseLevel.IsAvailable(
"RAP Pattern",
this)) {
148 Ac = Get< RCP<Matrix> >(coarseLevel,
"RAP Pattern");
151 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
157 bool doOptimizeStorage =
true;
159 const bool doTranspose =
true;
160 const bool doFillComplete =
true;
161 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
163 Ac = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, doTranspose, *AP, !doTranspose, Ac,
GetOStream(
Statistics2), doFillComplete, doOptimizeStorage,std::string(
"MueLu::R*(AP)-implicit-")+levelstr.str());
166 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
169 Ac = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac,
GetOStream(
Statistics2), doFillComplete, doOptimizeStorage,std::string(
"MueLu::R*(AP)-explicit-")+levelstr.str());
175 RCP<ParameterList> params = rcp(
new ParameterList());;
176 params->set(
"printLoadBalancingInfo",
true);
177 params->set(
"printCommInfo",
true);
181 Set(coarseLevel,
"A", Ac);
182 if (pL.get<
bool>(
"Keep RAP Pattern"))
183 Set(coarseLevel,
"RAP Pattern", Ac);
191 RCP<const FactoryBase> fac = *it;
192 GetOStream(
Runtime0) <<
"RAPFactory: call transfer factory: " << fac->description() << std::endl;
193 fac->CallBuild(coarseLevel);
229 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
232 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal");
233 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal");
235 if (!checkAc && !repairZeroDiagonals)
238 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
240 Teuchos::RCP<Teuchos::ParameterList> p = Teuchos::rcp(
new Teuchos::ParameterList());
241 p->set(
"DoOptimizeStorage",
true);
243 RCP<const Map> rowMap = Ac->getRowMap();
244 RCP<Vector> diagVec = VectorFactory::Build(rowMap);
245 Ac->getLocalDiagCopy(*diagVec);
248 Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
250 for (
size_t i = 0; i < rowMap->getNodeNumElements(); i++) {
251 if (diagVal[i] == zero) {
256 MueLu_sumAll(rowMap->getComm(), Teuchos::as<GO>(lZeroDiags), gZeroDiags);
258 if (repairZeroDiagonals && gZeroDiags > 0) {
273 RCP<Matrix> fixDiagMatrix = MatrixFactory::Build(rowMap, 1);
274 for (
size_t r = 0; r < rowMap->getNodeNumElements(); r++) {
275 if (diagVal[r] == zero) {
276 GO grid = rowMap->getGlobalElement(r);
277 Teuchos::ArrayRCP<GO> indout(1,grid);
278 Teuchos::ArrayRCP<SC> valout(1, one);
279 fixDiagMatrix->insertGlobalValues(grid,indout.view(0, 1), valout.view(0, 1));
283 Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer(
"CheckRepairMainDiagonal: fillComplete1"));
284 if(rowMap->lib() == Xpetra::UseTpetra) Ac->resumeFill();
287 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Ac,
false, 1.0, *fixDiagMatrix, 1.0);
288 if (Ac->IsView(
"stridedMaps"))
289 fixDiagMatrix->CreateView(
"stridedMaps", Ac);
297 Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer(
"CheckRepairMainDiagonal: fillComplete2"));
298 if(rowMap->lib() == Xpetra::UseTpetra) Ac->resumeFill();
307 GetOStream(
Warnings0) <<
"RAPFactory (WARNING): " << (repairZeroDiagonals ?
"repaired " :
"found ")
308 << gZeroDiags <<
" zeros on main diagonal of Ac." << std::endl;
310 #ifdef HAVE_MUELU_DEBUG // only for debugging 312 Ac->getLocalDiagCopy(*diagVec);
313 Teuchos::ArrayRCP< Scalar > diagVal2 = diagVec->getDataNonConst(0);
314 for (
size_t r = 0; r < Ac->getRowMap()->getNodeNumElements(); r++) {
315 if (diagVal2[r] == zero) {
316 GetOStream(
Errors,-1) <<
"Error: there are zeros left on diagonal after repair..." << std::endl;
323 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
326 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null,
Exceptions::BadCast,
327 "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. " 328 "This is very strange. (Note: you can remove this exception if there's a good reason for)");
335 #define MUELU_RAPFACTORY_SHORT 336 #endif // MUELU_RAPFACTORY_DEF_HPP Important warning messages (one line)
void Keep(const std::string &ename, const FactoryBase *factory)
Request to keep variable 'ename' generated by 'factory' after the setup phase.
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
std::vector< RCP< const FactoryBase > > transferFacts_
list of user-defined transfer Factories
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
One-liner description of what is happening.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
Namespace for MueLu classes and methods.
virtual const Teuchos::ParameterList & GetParameterList() const
Print even more statistics.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void Set(Level &level, const std::string &varName, const T &data) const
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
bool IsPrint(MsgType type, int thisProcRankOnly=-1) const
Find out whether we need to print out information for a specific message type.
void CheckRepairMainDiagonal(RCP< Matrix > &Ac) const
int GetLevelID() const
Return level number.
Exception throws to report errors in the internal logical of the program.
void Input(Level &level, const std::string &varName) const
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.