MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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;
177  {
178  SubFactoryMonitor m3(*this, "Diagonal Extraction", coarseLevel);
179  invDiag = Utilities_kokkos::GetMatrixDiagonalInverse(*A);
180  }
181  SC omega = dampingFactor / lambdaMax;
182  TEUCHOS_TEST_FOR_EXCEPTION(!std::isfinite(Teuchos::ScalarTraits<SC>::magnitude(omega)), Exceptions::RuntimeError, "Prolongator damping factor needs to be finite.");
183 
184  {
185  SubFactoryMonitor m3(*this, "Xpetra::IteratorOps::Jacobi", coarseLevel);
186  // finalP = Ptent + (I - \omega D^{-1}A) Ptent
187  finalP = Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(Statistics2), std::string("MueLu::SaP-") + toString(coarseLevel.GetLevelID()), APparams);
188  }
189  }
190 
191  } else {
192  finalP = Ptent;
193  }
194 
195  // Level Set
196  if (!restrictionMode_) {
197  // prolongation factory is in prolongation mode
198  if(!finalP.is_null()) {std::ostringstream oss; oss << "P_" << coarseLevel.GetLevelID(); finalP->setObjectLabel(oss.str());}
199  Set(coarseLevel, "P", finalP);
200 
201  // NOTE: EXPERIMENTAL
202  if (Ptent->IsView("stridedMaps"))
203  finalP->CreateView("stridedMaps", Ptent);
204 
205  } else {
206  // prolongation factory is in restriction mode
207  RCP<Matrix> R = Utilities_kokkos::Transpose(*finalP, true);
208  Set(coarseLevel, "R", R);
209  if(!R.is_null()) {std::ostringstream oss; oss << "R_" << coarseLevel.GetLevelID(); R->setObjectLabel(oss.str());}
210 
211  // NOTE: EXPERIMENTAL
212  if (Ptent->IsView("stridedMaps"))
213  R->CreateView("stridedMaps", Ptent, true);
214  }
215 
216  if (IsPrint(Statistics2)) {
217  RCP<ParameterList> params = rcp(new ParameterList());
218  params->set("printLoadBalancingInfo", true);
219  params->set("printCommInfo", true);
220  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*finalP, (!restrictionMode_ ? "P" : "R"), params);
221  }
222 
223  } //Build()
224 
225 } //namespace MueLu
226 
227 #endif // HAVE_MUELU_KOKKOS_REFACTOR
228 #endif // MUELU_SAPFACTORY_KOKKOS_DEF_HPP
229 
230 //TODO: restrictionMode_ should use the parameter list.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print more statistics.
One-liner description of what is happening.
Print even more statistics.
Print statistics that do not involve significant additional computation.
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)