46 #ifndef MUELU_SAPFACTORY_DEF_HPP 47 #define MUELU_SAPFACTORY_DEF_HPP 49 #include <Xpetra_Matrix.hpp> 50 #include <Xpetra_IteratorOps.hpp> 59 #include "MueLu_PerfUtils.hpp" 61 #include "MueLu_TentativePFactory.hpp" 62 #include "MueLu_Utilities.hpp" 66 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
68 RCP<ParameterList> validParamList = rcp(
new ParameterList());
70 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 74 #undef SET_VALID_ENTRY 76 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
77 validParamList->set< RCP<const FactoryBase> >(
"P", Teuchos::null,
"Tentative prolongator factory");
80 ParameterList norecurse;
81 norecurse.disableRecursiveValidation();
82 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
85 return validParamList;
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 Input(fineLevel,
"A");
94 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
95 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
99 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
101 return BuildP(fineLevel, coarseLevel);
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 const std::string prefix =
"MueLu::SaPFactory(" + levelIDs +
"): ";
112 typedef typename Teuchos::ScalarTraits<SC>::magnitudeType Magnitude;
117 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
118 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
119 const ParameterList& pL = GetParameterList();
122 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
123 RCP<Matrix> Ptent = coarseLevel.
Get< RCP<Matrix> >(
"P", initialPFact.get());
125 if (restrictionMode_) {
134 RCP<ParameterList> APparams;
135 if(pL.isSublist(
"matrixmatrix: kernel params"))
136 APparams=rcp(
new ParameterList(pL.sublist(
"matrixmatrix: kernel params")));
138 APparams= rcp(
new ParameterList);
140 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
141 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
143 APparams = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data",
this);
145 if (APparams->isParameter(
"graph"))
146 finalP = APparams->get< RCP<Matrix> >(
"graph");
149 APparams->set(
"compute global constants: temporaries",APparams->get(
"compute global constants: temporaries",
false));
150 APparams->set(
"compute global constants",APparams->get(
"compute global constants",
false));
152 SC
dampingFactor = as<SC>(pL.get<
double>(
"sa: damping factor"));
153 LO maxEigenIterations = as<LO>(pL.get<
int> (
"sa: eigenvalue estimate num iterations"));
154 bool estimateMaxEigen = pL.get<
bool> (
"sa: calculate eigenvalue estimate");
160 lambdaMax = A->GetMaxEigenvalueEstimate();
161 if (lambdaMax == -Teuchos::ScalarTraits<SC>::one() || estimateMaxEigen) {
162 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = "<< maxEigenIterations <<
")" << std::endl;
163 Magnitude stopTol = 1e-4;
165 A->SetMaxEigenvalueEstimate(lambdaMax);
167 GetOStream(
Statistics1) <<
"Using cached max eigenvalue estimate" << std::endl;
179 finalP = Xpetra::IteratorOps<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP,
180 GetOStream(
Statistics2), std::string(
"MueLu::SaP-")+levelIDs, APparams);
188 if (!restrictionMode_) {
190 Set(coarseLevel,
"P", finalP);
192 APparams->set(
"graph", finalP);
193 Set(coarseLevel,
"AP reuse data", APparams);
196 if (Ptent->IsView(
"stridedMaps"))
197 finalP->CreateView(
"stridedMaps", Ptent);
200 RCP<ParameterList> params = rcp(
new ParameterList());
201 params->set(
"printLoadBalancingInfo",
true);
202 params->set(
"printCommInfo",
true);
214 Set(coarseLevel,
"R", R);
217 if (Ptent->IsView(
"stridedMaps"))
218 R->CreateView(
"stridedMaps", Ptent,
true);
221 RCP<ParameterList> params = rcp(
new ParameterList());
222 params->set(
"printLoadBalancingInfo",
true);
223 params->set(
"printCommInfo",
true);
232 #endif // MUELU_SAPFACTORY_DEF_HPP static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixDiagonalInverse(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100)
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
scalar_type dampingFactor
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
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.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
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.
#define SET_VALID_ENTRY(name)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
static Scalar PowerMethod(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)