MueLu  Version of the Day
MueLu_SaPFactory_kokkos_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_SAPFACTORY_KOKKOS_DEF_HPP
47 #define MUELU_SAPFACTORY_KOKKOS_DEF_HPP
48 
49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
50 
52 
53 #include <Xpetra_Matrix.hpp>
54 #include <Xpetra_IteratorOps.hpp>
55 
57 #include "MueLu_Level.hpp"
58 #include "MueLu_MasterList.hpp"
59 #include "MueLu_Monitor.hpp"
60 #include "MueLu_PerfUtils.hpp"
62 #include "MueLu_TentativePFactory.hpp"
63 #include "MueLu_Utilities_kokkos.hpp"
64 
65 #include <sstream>
66 
67 namespace MueLu {
68 
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());
72 
73 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
74  SET_VALID_ENTRY("sa: damping factor");
75  SET_VALID_ENTRY("sa: calculate eigenvalue estimate");
76  SET_VALID_ENTRY("sa: eigenvalue estimate num iterations");
77 #undef SET_VALID_ENTRY
78 
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");
81 
82  // Make sure we don't recursively validate options for the matrixmatrix kernels
83  ParameterList norecurse;
84  norecurse.disableRecursiveValidation();
85  validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
86 
87 
88  return validParamList;
89  }
90 
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");
94 
95  // Get default tentative prolongator factory
96  // Getting it that way ensure that the same factory instance will be used for both SaPFactory_kokkos and NullspaceFactory.
97  RCP<const FactoryBase> initialPFact = GetFactory("P");
98  if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent"); }
99  coarseLevel.DeclareInput("P", initialPFact.get(), this);
100  }
101 
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);
105  }
106 
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);
110 
111  // Add debugging information
112  DeviceType::execution_space::print_configuration(GetOStream(Runtime1));
113 
114  typedef typename Teuchos::ScalarTraits<SC>::magnitudeType Magnitude;
115 
116  // Get default tentative prolongator factory
117  // Getting it that way ensure that the same factory instance will be used for both SaPFactory_kokkos and NullspaceFactory.
118  // -- Warning: Do not use directly initialPFact_. Use initialPFact instead everywhere!
119  RCP<const FactoryBase> initialPFact = GetFactory("P");
120  if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent"); }
121  const ParameterList& pL = GetParameterList();
122 
123  // Level Get
124  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
125  RCP<Matrix> Ptent = coarseLevel.Get< RCP<Matrix> >("P", initialPFact.get());
126 
127  if(restrictionMode_) {
128  SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
129 
130  A = Utilities_kokkos::Transpose(*A, true); // build transpose of A explicitly
131  }
132 
133  //Build final prolongator
134  RCP<Matrix> finalP; // output
135 
136  // Reuse pattern if available
137  RCP<ParameterList> APparams;
138  if(pL.isSublist("matrixmatrix: kernel params"))
139  APparams=rcp(new ParameterList(pL.sublist("matrixmatrix: kernel params")));
140  else
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;
144 
145  APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
146 
147  if (APparams->isParameter("graph"))
148  finalP = APparams->get< RCP<Matrix> >("graph");
149  }
150  // By default, we don't need global constants for SaP
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));
153 
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");
157  if (dampingFactor != Teuchos::ScalarTraits<SC>::zero()) {
158 
159  SC lambdaMax;
160  {
161  SubFactoryMonitor m2(*this, "Eigenvalue estimate", coarseLevel);
162  lambdaMax = A->GetMaxEigenvalueEstimate();
163  if (lambdaMax == -Teuchos::ScalarTraits<SC>::one() || estimateMaxEigen) {
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);
168  } else {
169  GetOStream(Statistics1) << "Using cached max eigenvalue estimate" << std::endl;
170  }
171  GetOStream(Statistics0) << "Prolongator damping factor = " << dampingFactor/lambdaMax << " (" << dampingFactor << " / " << lambdaMax << ")" << std::endl;
172  }
173 
174  {
175  SubFactoryMonitor m2(*this, "Fused (I-omega*D^{-1} A)*Ptent", coarseLevel);
176  RCP<Vector> invDiag = Utilities_kokkos::GetMatrixDiagonalInverse(*A);
177 
178  SC omega = dampingFactor / lambdaMax;
179 
180  // finalP = Ptent + (I - \omega D^{-1}A) Ptent
181  finalP = Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(Statistics2), std::string("MueLu::SaP-") + toString(coarseLevel.GetLevelID()), APparams);
182  }
183 
184  } else {
185  finalP = Ptent;
186  }
187 
188  // Level Set
189  if (!restrictionMode_) {
190  // prolongation factory is in prolongation mode
191  Set(coarseLevel, "P", finalP);
192 
193  // NOTE: EXPERIMENTAL
194  if (Ptent->IsView("stridedMaps"))
195  finalP->CreateView("stridedMaps", Ptent);
196 
197  } else {
198  // prolongation factory is in restriction mode
199  RCP<Matrix> R = Utilities_kokkos::Transpose(*finalP, true);
200  Set(coarseLevel, "R", R);
201 
202  // NOTE: EXPERIMENTAL
203  if (Ptent->IsView("stridedMaps"))
204  R->CreateView("stridedMaps", Ptent, true);
205  }
206 
207  if (IsPrint(Statistics2)) {
208  RCP<ParameterList> params = rcp(new ParameterList());
209  params->set("printLoadBalancingInfo", true);
210  params->set("printCommInfo", true);
211  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*finalP, (!restrictionMode_ ? "P" : "R"), params);
212  }
213 
214  } //Build()
215 
216 } //namespace MueLu
217 
218 #endif // HAVE_MUELU_KOKKOS_REFACTOR
219 #endif // MUELU_SAPFACTORY_KOKKOS_DEF_HPP
220 
221 //TODO: restrictionMode_ should use the parameter list.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
scalar_type dampingFactor
Print more statistics.
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)