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_MatrixUtils.hpp> 56 #include <Xpetra_TripleMatrixMultiply.hpp> 57 #include <Xpetra_Vector.hpp> 58 #include <Xpetra_VectorFactory.hpp> 64 #include "MueLu_PerfUtils.hpp" 70 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 : hasDeclaredInput_(false) { }
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 RCP<ParameterList> validParamList = rcp(
new ParameterList());
78 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 82 #undef SET_VALID_ENTRY 83 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A used during the prolongator smoothing process");
84 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Prolongator factory");
85 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Restrictor factory");
87 validParamList->set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
88 validParamList->set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
91 ParameterList norecurse;
92 norecurse.disableRecursiveValidation();
93 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
95 return validParamList;
98 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 const Teuchos::ParameterList& pL = GetParameterList();
101 if (pL.get<
bool>(
"transpose: use implicit") ==
false)
102 Input(coarseLevel,
"R");
104 Input(fineLevel,
"A");
105 Input(coarseLevel,
"P");
108 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
109 (*it)->CallDeclareInput(coarseLevel);
111 hasDeclaredInput_ =
true;
114 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 const bool doTranspose =
true;
117 const bool doFillComplete =
true;
118 const bool doOptimizeStorage =
true;
121 std::ostringstream levelstr;
126 "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
128 const Teuchos::ParameterList& pL = GetParameterList();
129 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
130 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P"), AP, Ac;
133 if (pL.get<
bool>(
"rap: triple product") ==
false) {
135 RCP<ParameterList> APparams;
136 if(pL.isSublist(
"matrixmatrix: kernel params"))
137 APparams=rcp(
new ParameterList(pL.sublist(
"matrixmatrix: kernel params")));
139 APparams= rcp(
new ParameterList);
142 APparams->set(
"compute global constants: temporaries",APparams->get(
"compute global constants: temporaries",
false));
143 APparams->set(
"compute global constants",APparams->get(
"compute global constants",
false));
145 if (coarseLevel.IsAvailable(
"AP reuse data",
this)) {
146 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
148 APparams = coarseLevel.Get< RCP<ParameterList> >(
"AP reuse data",
this);
150 if (APparams->isParameter(
"graph"))
151 AP = APparams->get< RCP<Matrix> >(
"graph");
157 AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(
Statistics2),
158 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::A*P-")+levelstr.str(), APparams);
162 RCP<ParameterList> RAPparams;
163 if(pL.isSublist(
"matrixmatrix: kernel params")) RAPparams=rcp(
new ParameterList(pL.sublist(
"matrixmatrix: kernel params")));
164 else RAPparams= rcp(
new ParameterList);
168 if (coarseLevel.IsAvailable(
"RAP reuse data",
this)) {
169 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous RAP data" << std::endl;
171 RAPparams = coarseLevel.Get< RCP<ParameterList> >(
"RAP reuse data",
this);
173 if (RAPparams->isParameter(
"graph"))
174 Ac = RAPparams->get< RCP<Matrix> >(
"graph");
178 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
182 RAPparams->set(
"compute global constants: temporaries",RAPparams->get(
"compute global constants: temporaries",
false));
183 RAPparams->set(
"compute global constants",
true);
189 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
192 Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
193 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
196 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
200 Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
201 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
203 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal") || pL.get<
bool>(
"rap: fix zero diagonals");
204 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal")|| pL.get<
bool>(
"rap: fix zero diagonals"); ;
205 if (checkAc || repairZeroDiagonals)
206 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(
Warnings1));
209 RCP<ParameterList> params = rcp(
new ParameterList());;
210 params->set(
"printLoadBalancingInfo",
true);
211 params->set(
"printCommInfo",
true);
215 Set(coarseLevel,
"A", Ac);
217 APparams->set(
"graph", AP);
218 Set(coarseLevel,
"AP reuse data", APparams);
219 RAPparams->set(
"graph", Ac);
220 Set(coarseLevel,
"RAP reuse data", RAPparams);
222 RCP<ParameterList> RAPparams;
223 RAPparams= rcp(
new ParameterList);
226 RAPparams->set(
"compute global constants: temporaries",RAPparams->get(
"compute global constants: temporaries",
false));
227 RAPparams->set(
"compute global constants",
true);
229 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
231 Ac = MatrixFactory::Build(P->getDomainMap(), Teuchos::as<LO>(0));
235 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
236 MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
237 doOptimizeStorage, labelstr+std::string(
"MueLu::R*A*P-implicit-")+levelstr.str(),
241 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
242 Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));
246 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
247 MultiplyRAP(*R, !doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
248 doOptimizeStorage, labelstr+std::string(
"MueLu::R*A*P-explicit-")+levelstr.str(),
251 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal") || pL.get<
bool>(
"rap: fix zero diagonals");
252 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal")|| pL.get<
bool>(
"rap: fix zero diagonals"); ;
253 if (checkAc || repairZeroDiagonals)
254 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(
Warnings1));
257 RCP<ParameterList> params = rcp(
new ParameterList());;
258 params->set(
"printLoadBalancingInfo",
true);
259 params->set(
"printCommInfo",
true);
263 Set(coarseLevel,
"A", Ac);
272 if (transferFacts_.begin() != transferFacts_.end()) {
276 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
277 RCP<const FactoryBase> fac = *it;
278 GetOStream(
Runtime0) <<
"RAPFactory: call transfer factory: " << fac->description() << std::endl;
279 fac->CallBuild(coarseLevel);
314 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
317 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null,
Exceptions::BadCast,
318 "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. " 319 "This is very strange. (Note: you can remove this exception if there's a good reason for)");
320 TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_,
Exceptions::RuntimeError,
"MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
321 transferFacts_.push_back(factory);
326 #define MUELU_RAPFACTORY_SHORT 327 #endif // MUELU_RAPFACTORY_DEF_HPP #define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
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.
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.
Print skeleton for the run, i.e. factory calls and used parameters.
Print even more statistics.
int GetLevelID() const
Return level number.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.