46 #ifndef MUELU_SAPFACTORY_KOKKOS_DEF_HPP 47 #define MUELU_SAPFACTORY_KOKKOS_DEF_HPP 49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 60 #include "MueLu_PerfUtils.hpp" 62 #include "MueLu_TentativePFactory.hpp" 63 #include "MueLu_Utilities_kokkos.hpp" 69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class DeviceType>
70 RCP<const ParameterList> SaPFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::GetValidParameterList()
const {
71 RCP<ParameterList> validParamList =
rcp(
new ParameterList());
73 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 77 #undef SET_VALID_ENTRY 79 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
80 validParamList->set< RCP<const FactoryBase> >(
"P", Teuchos::null,
"Tentative prolongator factory");
83 ParameterList norecurse;
84 norecurse.disableRecursiveValidation();
85 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
88 return validParamList;
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class DeviceType>
92 void SaPFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::DeclareInput(Level& fineLevel, Level& coarseLevel)
const {
93 Input(fineLevel,
"A");
97 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
98 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory(
"Ptent"); }
99 coarseLevel.DeclareInput(
"P", initialPFact.get(),
this);
102 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class DeviceType>
103 void SaPFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::Build(Level& fineLevel, Level& coarseLevel)
const {
104 return BuildP(fineLevel, coarseLevel);
107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class DeviceType>
108 void SaPFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::BuildP(Level& fineLevel, Level& coarseLevel)
const {
109 FactoryMonitor m(*
this,
"Prolongator smoothing", coarseLevel);
112 DeviceType::execution_space::print_configuration(GetOStream(
Runtime1));
119 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
120 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory(
"Ptent"); }
121 const ParameterList& pL = GetParameterList();
124 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
125 RCP<Matrix> Ptent = coarseLevel.Get< RCP<Matrix> >(
"P", initialPFact.get());
127 if(restrictionMode_) {
128 SubFactoryMonitor m2(*
this,
"Transpose A", coarseLevel);
130 A = Utilities_kokkos::Transpose(*A,
true);
137 RCP<ParameterList> APparams;
138 if(pL.isSublist(
"matrixmatrix: kernel params"))
139 APparams=
rcp(
new ParameterList(pL.sublist(
"matrixmatrix: kernel params")));
141 APparams=
rcp(
new ParameterList);
142 if (coarseLevel.IsAvailable(
"AP reuse data",
this)) {
143 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
145 APparams = coarseLevel.Get< RCP<ParameterList> >(
"AP reuse data",
this);
147 if (APparams->isParameter(
"graph"))
148 finalP = APparams->get< RCP<Matrix> >(
"graph");
151 APparams->set(
"compute global constants: temporaries",APparams->get(
"compute global constants: temporaries",
false));
152 APparams->set(
"compute global constants",APparams->get(
"compute global constants",
false));
154 SC dampingFactor = as<SC>(pL.get<
double>(
"sa: damping factor"));
155 LO maxEigenIterations = as<LO>(pL.get<
int>(
"sa: eigenvalue estimate num iterations"));
156 bool estimateMaxEigen = pL.get<
bool>(
"sa: calculate eigenvalue estimate");
161 SubFactoryMonitor m2(*
this,
"Eigenvalue estimate", coarseLevel);
162 lambdaMax = A->GetMaxEigenvalueEstimate();
164 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = "<< maxEigenIterations <<
")" << std::endl;
165 Magnitude stopTol = 1e-4;
166 lambdaMax = Utilities_kokkos::PowerMethod(*A,
true, maxEigenIterations, stopTol);
167 A->SetMaxEigenvalueEstimate(lambdaMax);
169 GetOStream(
Statistics1) <<
"Using cached max eigenvalue estimate" << std::endl;
171 GetOStream(
Statistics0) <<
"Prolongator damping factor = " << dampingFactor/lambdaMax <<
" (" << dampingFactor <<
" / " << lambdaMax <<
")" << std::endl;
175 SubFactoryMonitor m2(*
this,
"Fused (I-omega*D^{-1} A)*Ptent", coarseLevel);
176 RCP<Vector> invDiag = Utilities_kokkos::GetMatrixDiagonalInverse(*A);
178 SC omega = dampingFactor / lambdaMax;
181 finalP =
Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(
Statistics2), std::string(
"MueLu::SaP-") +
toString(coarseLevel.GetLevelID()), APparams);
189 if (!restrictionMode_) {
191 Set(coarseLevel,
"P", finalP);
194 if (Ptent->IsView(
"stridedMaps"))
195 finalP->CreateView(
"stridedMaps", Ptent);
199 RCP<Matrix> R = Utilities_kokkos::Transpose(*finalP,
true);
200 Set(coarseLevel,
"R", R);
203 if (Ptent->IsView(
"stridedMaps"))
204 R->CreateView(
"stridedMaps", Ptent,
true);
208 RCP<ParameterList> params =
rcp(
new ParameterList());
209 params->set(
"printLoadBalancingInfo",
true);
210 params->set(
"printCommInfo",
true);
218 #endif // HAVE_MUELU_KOKKOS_REFACTOR 219 #endif // MUELU_SAPFACTORY_KOKKOS_DEF_HPP
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
Print skeleton for the run, i.e. factory calls and used parameters.
Print even more statistics.
Print statistics that do not involve significant additional computation.
static RCP< Matrix > Jacobi(SC omega, const Vector &Dinv, const Matrix &A, const Matrix &B, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, const std::string &label, RCP< ParameterList > ¶ms)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Description of what is happening (more verbose)
#define SET_VALID_ENTRY(name)