46 #ifndef MUELU_SAPFACTORY_KOKKOS_DEF_HPP 47 #define MUELU_SAPFACTORY_KOKKOS_DEF_HPP 49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 53 #include <Xpetra_Matrix.hpp> 54 #include <Xpetra_IteratorOps.hpp> 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");
82 return validParamList;
85 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class DeviceType>
86 void SaPFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::DeclareInput(Level& fineLevel, Level& coarseLevel)
const {
87 Input(fineLevel,
"A");
91 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
92 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory(
"Ptent"); }
93 coarseLevel.DeclareInput(
"P", initialPFact.get(),
this);
96 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class DeviceType>
97 void SaPFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::Build(Level& fineLevel, Level& coarseLevel)
const {
98 return BuildP(fineLevel, coarseLevel);
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class DeviceType>
102 void SaPFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::BuildP(Level& fineLevel, Level& coarseLevel)
const {
103 FactoryMonitor m(*
this,
"Prolongator smoothing", coarseLevel);
106 DeviceType::execution_space::print_configuration(GetOStream(
Runtime1));
108 typedef typename Teuchos::ScalarTraits<SC>::magnitudeType Magnitude;
113 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
114 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory(
"Ptent"); }
117 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
118 RCP<Matrix> Ptent = coarseLevel.Get< RCP<Matrix> >(
"P", initialPFact.get());
120 if(restrictionMode_) {
121 SubFactoryMonitor m2(*
this,
"Transpose A", coarseLevel);
123 A = Utilities_kokkos::Transpose(*A,
true);
130 RCP<ParameterList> APparams = rcp(
new ParameterList);
131 if (coarseLevel.IsAvailable(
"AP reuse data",
this)) {
132 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
134 APparams = coarseLevel.Get< RCP<ParameterList> >(
"AP reuse data",
this);
136 if (APparams->isParameter(
"graph"))
137 finalP = APparams->get< RCP<Matrix> >(
"graph");
140 const ParameterList& pL = GetParameterList();
141 SC dampingFactor = as<SC>(pL.get<
double>(
"sa: damping factor"));
142 LO maxEigenIterations = as<LO>(pL.get<
int>(
"sa: eigenvalue estimate num iterations"));
143 bool estimateMaxEigen = pL.get<
bool>(
"sa: calculate eigenvalue estimate");
144 if (dampingFactor != Teuchos::ScalarTraits<SC>::zero()) {
148 SubFactoryMonitor m2(*
this,
"Eigenvalue estimate", coarseLevel);
149 lambdaMax = A->GetMaxEigenvalueEstimate();
150 if (lambdaMax == -Teuchos::ScalarTraits<SC>::one() || estimateMaxEigen) {
151 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = "<< maxEigenIterations <<
")" << std::endl;
152 Magnitude stopTol = 1e-4;
153 lambdaMax = Utilities_kokkos::PowerMethod(*A,
true, maxEigenIterations, stopTol);
154 A->SetMaxEigenvalueEstimate(lambdaMax);
156 GetOStream(
Statistics1) <<
"Using cached max eigenvalue estimate" << std::endl;
158 GetOStream(
Statistics0) <<
"Prolongator damping factor = " << dampingFactor/lambdaMax <<
" (" << dampingFactor <<
" / " << lambdaMax <<
")" << std::endl;
162 SubFactoryMonitor m2(*
this,
"Fused (I-omega*D^{-1} A)*Ptent", coarseLevel);
163 RCP<Vector> invDiag = Utilities_kokkos::GetMatrixDiagonalInverse(*A);
165 SC omega = dampingFactor / lambdaMax;
168 finalP = Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(
Statistics2), std::string(
"MueLu::SaP-") +
toString(coarseLevel.GetLevelID()), APparams);
176 if (!restrictionMode_) {
178 Set(coarseLevel,
"P", finalP);
181 if (Ptent->IsView(
"stridedMaps"))
182 finalP->CreateView(
"stridedMaps", Ptent);
186 RCP<Matrix> R = Utilities_kokkos::Transpose(*finalP,
true);
187 Set(coarseLevel,
"R", R);
190 if (Ptent->IsView(
"stridedMaps"))
191 R->CreateView(
"stridedMaps", Ptent,
true);
195 RCP<ParameterList> params = rcp(
new ParameterList());
196 params->set(
"printLoadBalancingInfo",
true);
197 params->set(
"printCommInfo",
true);
205 #endif // HAVE_MUELU_KOKKOS_REFACTOR 206 #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 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)