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);
129 const ParameterList& pL = GetParameterList();
130 SC dampingFactor = as<SC>(pL.get<
double>(
"sa: damping factor"));
131 LO maxEigenIterations = as<LO>(pL.get<
int>(
"sa: eigenvalue estimate num iterations"));
132 bool estimateMaxEigen = pL.get<
bool>(
"sa: calculate eigenvalue estimate");
133 if (dampingFactor != Teuchos::ScalarTraits<SC>::zero()) {
137 SubFactoryMonitor m2(*
this,
"Eigenvalue estimate", coarseLevel);
138 lambdaMax = A->GetMaxEigenvalueEstimate();
139 if (lambdaMax == -Teuchos::ScalarTraits<SC>::one() || estimateMaxEigen) {
140 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = "<< maxEigenIterations <<
")" << std::endl;
141 Magnitude stopTol = 1e-4;
142 lambdaMax = Utilities_kokkos::PowerMethod(*A,
true, maxEigenIterations, stopTol);
143 A->SetMaxEigenvalueEstimate(lambdaMax);
145 GetOStream(
Statistics1) <<
"Using cached max eigenvalue estimate" << std::endl;
147 GetOStream(
Statistics0) <<
"Prolongator damping factor = " << dampingFactor/lambdaMax <<
" (" << dampingFactor <<
" / " << lambdaMax <<
")" << std::endl;
151 SubFactoryMonitor m2(*
this,
"Fused (I-omega*D^{-1} A)*Ptent", coarseLevel);
152 RCP<Vector> invDiag = Utilities_kokkos::GetMatrixDiagonalInverse(*A);
154 SC omega = dampingFactor / lambdaMax;
157 finalP = Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(
Statistics2), std::string(
"MueLu::SaP-") +
toString(coarseLevel.GetLevelID()));
165 if (!restrictionMode_) {
167 Set(coarseLevel,
"P", finalP);
170 if (Ptent->IsView(
"stridedMaps"))
171 finalP->CreateView(
"stridedMaps", Ptent);
175 RCP<Matrix> R = Utilities_kokkos::Transpose(*finalP,
true);
176 Set(coarseLevel,
"R", R);
179 if (Ptent->IsView(
"stridedMaps"))
180 R->CreateView(
"stridedMaps", Ptent,
true);
184 RCP<ParameterList> params = rcp(
new ParameterList());
185 params->set(
"printLoadBalancingInfo",
true);
186 params->set(
"printCommInfo",
true);
194 #endif // HAVE_MUELU_KOKKOS_REFACTOR 195 #endif // MUELU_SAPFACTORY_KOKKOS_DEF_HPP std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Namespace for MueLu classes and methods.
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)