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  return validParamList;
83  }
84 
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");
88 
89  // Get default tentative prolongator factory
90  // Getting it that way ensure that the same factory instance will be used for both SaPFactory_kokkos and NullspaceFactory.
91  RCP<const FactoryBase> initialPFact = GetFactory("P");
92  if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent"); }
93  coarseLevel.DeclareInput("P", initialPFact.get(), this);
94  }
95 
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);
99  }
100 
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);
104 
105  // Add debugging information
106  DeviceType::execution_space::print_configuration(GetOStream(Runtime1));
107 
108  typedef typename Teuchos::ScalarTraits<SC>::magnitudeType Magnitude;
109 
110  // Get default tentative prolongator factory
111  // Getting it that way ensure that the same factory instance will be used for both SaPFactory_kokkos and NullspaceFactory.
112  // -- Warning: Do not use directly initialPFact_. Use initialPFact instead everywhere!
113  RCP<const FactoryBase> initialPFact = GetFactory("P");
114  if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent"); }
115 
116  // Level Get
117  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
118  RCP<Matrix> Ptent = coarseLevel.Get< RCP<Matrix> >("P", initialPFact.get());
119 
120  if(restrictionMode_) {
121  SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
122 
123  A = Utilities_kokkos::Transpose(*A, true); // build transpose of A explicitly
124  }
125 
126  //Build final prolongator
127  RCP<Matrix> finalP; // output
128 
129  // Reuse pattern if available
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;
133 
134  APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
135 
136  if (APparams->isParameter("graph"))
137  finalP = APparams->get< RCP<Matrix> >("graph");
138  }
139 
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()) {
145 
146  SC lambdaMax;
147  {
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);
155  } else {
156  GetOStream(Statistics1) << "Using cached max eigenvalue estimate" << std::endl;
157  }
158  GetOStream(Statistics0) << "Prolongator damping factor = " << dampingFactor/lambdaMax << " (" << dampingFactor << " / " << lambdaMax << ")" << std::endl;
159  }
160 
161  {
162  SubFactoryMonitor m2(*this, "Fused (I-omega*D^{-1} A)*Ptent", coarseLevel);
163  RCP<Vector> invDiag = Utilities_kokkos::GetMatrixDiagonalInverse(*A);
164 
165  SC omega = dampingFactor / lambdaMax;
166 
167  // finalP = Ptent + (I - \omega D^{-1}A) Ptent
168  finalP = Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(Statistics2), std::string("MueLu::SaP-") + toString(coarseLevel.GetLevelID()), APparams);
169  }
170 
171  } else {
172  finalP = Ptent;
173  }
174 
175  // Level Set
176  if (!restrictionMode_) {
177  // prolongation factory is in prolongation mode
178  Set(coarseLevel, "P", finalP);
179 
180  // NOTE: EXPERIMENTAL
181  if (Ptent->IsView("stridedMaps"))
182  finalP->CreateView("stridedMaps", Ptent);
183 
184  } else {
185  // prolongation factory is in restriction mode
186  RCP<Matrix> R = Utilities_kokkos::Transpose(*finalP, true);
187  Set(coarseLevel, "R", R);
188 
189  // NOTE: EXPERIMENTAL
190  if (Ptent->IsView("stridedMaps"))
191  R->CreateView("stridedMaps", Ptent, true);
192  }
193 
194  if (IsPrint(Statistics1)) {
195  RCP<ParameterList> params = rcp(new ParameterList());
196  params->set("printLoadBalancingInfo", true);
197  params->set("printCommInfo", true);
198  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*finalP, (!restrictionMode_ ? "P" : "R"), params);
199  }
200 
201  } //Build()
202 
203 } //namespace MueLu
204 
205 #endif // HAVE_MUELU_KOKKOS_REFACTOR
206 #endif // MUELU_SAPFACTORY_KOKKOS_DEF_HPP
207 
208 //TODO: restrictionMode_ should use the parameter list.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Print more statistics.
LocalOrdinal LO
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 > &params)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Scalar SC
Description of what is happening (more verbose)
#define SET_VALID_ENTRY(name)