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  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()) {
134 
135  SC lambdaMax;
136  {
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);
144  } else {
145  GetOStream(Statistics1) << "Using cached max eigenvalue estimate" << std::endl;
146  }
147  GetOStream(Statistics0) << "Prolongator damping factor = " << dampingFactor/lambdaMax << " (" << dampingFactor << " / " << lambdaMax << ")" << std::endl;
148  }
149 
150  {
151  SubFactoryMonitor m2(*this, "Fused (I-omega*D^{-1} A)*Ptent", coarseLevel);
152  RCP<Vector> invDiag = Utilities_kokkos::GetMatrixDiagonalInverse(*A);
153 
154  SC omega = dampingFactor / lambdaMax;
155 
156  // finalP = Ptent + (I - \omega D^{-1}A) Ptent
157  finalP = Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(Statistics2), std::string("MueLu::SaP-") + toString(coarseLevel.GetLevelID()));
158  }
159 
160  } else {
161  finalP = Ptent;
162  }
163 
164  // Level Set
165  if (!restrictionMode_) {
166  // prolongation factory is in prolongation mode
167  Set(coarseLevel, "P", finalP);
168 
169  // NOTE: EXPERIMENTAL
170  if (Ptent->IsView("stridedMaps"))
171  finalP->CreateView("stridedMaps", Ptent);
172 
173  } else {
174  // prolongation factory is in restriction mode
175  RCP<Matrix> R = Utilities_kokkos::Transpose(*finalP, true);
176  Set(coarseLevel, "R", R);
177 
178  // NOTE: EXPERIMENTAL
179  if (Ptent->IsView("stridedMaps"))
180  R->CreateView("stridedMaps", Ptent, true);
181  }
182 
183  if (IsPrint(Statistics1)) {
184  RCP<ParameterList> params = rcp(new ParameterList());
185  params->set("printLoadBalancingInfo", true);
186  params->set("printCommInfo", true);
187  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*finalP, (!restrictionMode_ ? "P" : "R"), params);
188  }
189 
190  } //Build()
191 
192 } //namespace MueLu
193 
194 #endif // HAVE_MUELU_KOKKOS_REFACTOR
195 #endif // MUELU_SAPFACTORY_KOKKOS_DEF_HPP
196 
197 //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.
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)