MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_RAPFactory_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_RAPFACTORY_DEF_HPP
47#define MUELU_RAPFACTORY_DEF_HPP
48
49
50#include <sstream>
51
52#include <Xpetra_Matrix.hpp>
53#include <Xpetra_MatrixFactory.hpp>
54#include <Xpetra_MatrixMatrix.hpp>
55#include <Xpetra_MatrixUtils.hpp>
56#include <Xpetra_TripleMatrixMultiply.hpp>
57#include <Xpetra_Vector.hpp>
58#include <Xpetra_VectorFactory.hpp>
59#include <Xpetra_IO.hpp>
60
62
63#include "MueLu_MasterList.hpp"
64#include "MueLu_Monitor.hpp"
65#include "MueLu_PerfUtils.hpp"
67//#include "MueLu_Utilities.hpp"
68
69namespace MueLu {
70
71 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74
75 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77 RCP<ParameterList> validParamList = rcp(new ParameterList());
78
79#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
80 SET_VALID_ENTRY("transpose: use implicit");
81 SET_VALID_ENTRY("rap: triple product");
82 SET_VALID_ENTRY("rap: fix zero diagonals");
83 SET_VALID_ENTRY("rap: fix zero diagonals threshold");
84 SET_VALID_ENTRY("rap: fix zero diagonals replacement");
85 SET_VALID_ENTRY("rap: relative diagonal floor");
86#undef SET_VALID_ENTRY
87 validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A used during the prolongator smoothing process");
88 validParamList->set< RCP<const FactoryBase> >("P", null, "Prolongator factory");
89 validParamList->set< RCP<const FactoryBase> >("R", null, "Restrictor factory");
90
91 validParamList->set< bool > ("CheckMainDiagonal", false, "Check main diagonal for zeros");
92 validParamList->set< bool > ("RepairMainDiagonal", false, "Repair zeros on main diagonal");
93
94 // Make sure we don't recursively validate options for the matrixmatrix kernels
95 ParameterList norecurse;
96 norecurse.disableRecursiveValidation();
97 validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
98
99 return validParamList;
100 }
101
102 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
104 const Teuchos::ParameterList& pL = GetParameterList();
105 if (pL.get<bool>("transpose: use implicit") == false)
106 Input(coarseLevel, "R");
107
108 Input(fineLevel, "A");
109 Input(coarseLevel, "P");
110
111 // call DeclareInput of all user-given transfer factories
112 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
113 (*it)->CallDeclareInput(coarseLevel);
114
115 hasDeclaredInput_ = true;
116 }
117
118 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
120 const bool doTranspose = true;
121 const bool doFillComplete = true;
122 const bool doOptimizeStorage = true;
123 RCP<Matrix> Ac;
124 {
125 FactoryMonitor m(*this, "Computing Ac", coarseLevel);
126 std::ostringstream levelstr;
127 levelstr << coarseLevel.GetLevelID();
128 std::string labelstr = FormattingHelper::getColonLabel(coarseLevel.getObjectLabel());
129
130 TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_ == false, Exceptions::RuntimeError,
131 "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
132
133 const Teuchos::ParameterList& pL = GetParameterList();
134 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
135 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P"), AP;
136 // We don't have a valid P (e.g., # global aggregates = 0) so we bail.
137 // This level will ultimately be removed in MueLu_Hierarchy_defs.h via a resize()
138 if (P == Teuchos::null) {
139 Ac = Teuchos::null;
140 Set(coarseLevel, "A", Ac);
141 return;
142 }
143
144 bool isEpetra = A->getRowMap()->lib() == Xpetra::UseEpetra;
145 bool isGPU =
146#ifdef KOKKOS_ENABLE_CUDA
147 (typeid(Node).name() == typeid(Kokkos::Compat::KokkosCudaWrapperNode).name()) ||
148#endif
149#ifdef KOKKOS_ENABLE_HIP
150 (typeid(Node).name() == typeid(Kokkos::Compat::KokkosHIPWrapperNode).name()) ||
151#endif
152 false;
153
154 if (pL.get<bool>("rap: triple product") == false || isEpetra || isGPU) {
155 if (pL.get<bool>("rap: triple product") && isEpetra)
156 GetOStream(Warnings1) << "Switching from triple product to R x (A x P) since triple product has not been implemented for Epetra.\n";
157#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
158 if (pL.get<bool>("rap: triple product") && isGPU)
159 GetOStream(Warnings1) << "Switching from triple product to R x (A x P) since triple product has not been implemented for "
160 << Node::execution_space::name() << std::endl;
161#endif
162
163 // Reuse pattern if available (multiple solve)
164 RCP<ParameterList> APparams = rcp(new ParameterList);
165 if(pL.isSublist("matrixmatrix: kernel params"))
166 APparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
167
168 // By default, we don't need global constants for A*P
169 APparams->set("compute global constants: temporaries",APparams->get("compute global constants: temporaries",false));
170 APparams->set("compute global constants",APparams->get("compute global constants",false));
171
172 if (coarseLevel.IsAvailable("AP reuse data", this)) {
173 GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
174
175 APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
176
177 if (APparams->isParameter("graph"))
178 AP = APparams->get< RCP<Matrix> >("graph");
179 }
180
181 {
182 SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
183
184 AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(Statistics2),
185 doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::A*P-")+levelstr.str(), APparams);
186 }
187
188 // Reuse coarse matrix memory if available (multiple solve)
189 RCP<ParameterList> RAPparams = rcp(new ParameterList);
190 if(pL.isSublist("matrixmatrix: kernel params"))
191 RAPparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
192
193 if (coarseLevel.IsAvailable("RAP reuse data", this)) {
194 GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous RAP data" << std::endl;
195
196 RAPparams = coarseLevel.Get< RCP<ParameterList> >("RAP reuse data", this);
197
198 if (RAPparams->isParameter("graph"))
199 Ac = RAPparams->get< RCP<Matrix> >("graph");
200
201 // Some eigenvalue may have been cached with the matrix in the previous run.
202 // As the matrix values will be updated, we need to reset the eigenvalue.
203 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
204 }
205
206 // We *always* need global constants for the RAP, but not for the temps
207 RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
208 RAPparams->set("compute global constants",true);
209
210 // Allow optimization of storage.
211 // This is necessary for new faster Epetra MM kernels.
212 // Seems to work with matrix modifications to repair diagonal entries.
213
214 if (pL.get<bool>("transpose: use implicit") == true) {
215 SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
216
217 Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
218 doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
219
220 } else {
221 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
222
223 SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
224
225 Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
226 doFillComplete, doOptimizeStorage, labelstr+std::string("MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
227 }
228
229 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >("rap: relative diagonal floor")();
230 if(relativeFloor.size() > 0) {
231 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(Statistics2));
232 }
233
234 bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
235 bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
236 if (checkAc || repairZeroDiagonals) {
237 using magnitudeType = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
238 magnitudeType threshold;
239 if (pL.isType<magnitudeType>("rap: fix zero diagonals threshold"))
240 threshold = pL.get<magnitudeType>("rap: fix zero diagonals threshold");
241 else
242 threshold = Teuchos::as<magnitudeType>(pL.get<double>("rap: fix zero diagonals threshold"));
243 Scalar replacement = Teuchos::as<Scalar>(pL.get<double>("rap: fix zero diagonals replacement"));
244 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1), threshold, replacement);
245 }
246
247 if (IsPrint(Statistics2)) {
248 RCP<ParameterList> params = rcp(new ParameterList());;
249 params->set("printLoadBalancingInfo", true);
250 params->set("printCommInfo", true);
251 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
252 }
253
254 if(!Ac.is_null()) {std::ostringstream oss; oss << "A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
255 Set(coarseLevel, "A", Ac);
256
257 if (!isGPU) {
258 APparams->set("graph", AP);
259 Set(coarseLevel, "AP reuse data", APparams);
260 }
261 if (!isGPU) {
262 RAPparams->set("graph", Ac);
263 Set(coarseLevel, "RAP reuse data", RAPparams);
264 }
265 } else {
266 RCP<ParameterList> RAPparams = rcp(new ParameterList);
267 if(pL.isSublist("matrixmatrix: kernel params"))
268 RAPparams->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
269
270 if (coarseLevel.IsAvailable("RAP reuse data", this)) {
271 GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous RAP data" << std::endl;
272
273 RAPparams = coarseLevel.Get< RCP<ParameterList> >("RAP reuse data", this);
274
275 if (RAPparams->isParameter("graph"))
276 Ac = RAPparams->get< RCP<Matrix> >("graph");
277
278 // Some eigenvalue may have been cached with the matrix in the previous run.
279 // As the matrix values will be updated, we need to reset the eigenvalue.
280 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
281 }
282
283 // We *always* need global constants for the RAP, but not for the temps
284 RAPparams->set("compute global constants: temporaries",RAPparams->get("compute global constants: temporaries",false));
285 RAPparams->set("compute global constants",true);
286
287 if (pL.get<bool>("transpose: use implicit") == true) {
288
289 Ac = MatrixFactory::Build(P->getDomainMap(), Teuchos::as<LO>(0));
290
291 SubFactoryMonitor m2(*this, "MxMxM: R x A x P (implicit)", coarseLevel);
292
293 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
294 MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
295 doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-implicit-")+levelstr.str(),
296 RAPparams);
297 } else {
298 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
299 Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));
300
301 SubFactoryMonitor m2(*this, "MxMxM: R x A x P (explicit)", coarseLevel);
302
303 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
304 MultiplyRAP(*R, !doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
305 doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-explicit-")+levelstr.str(),
306 RAPparams);
307 }
308
309 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >("rap: relative diagonal floor")();
310 if(relativeFloor.size() > 0) {
311 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(Statistics2));
312 }
313
314 bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
315 bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
316 if (checkAc || repairZeroDiagonals) {
317 using magnitudeType = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
318 magnitudeType threshold;
319 if (pL.isType<magnitudeType>("rap: fix zero diagonals threshold"))
320 threshold = pL.get<magnitudeType>("rap: fix zero diagonals threshold");
321 else
322 threshold = Teuchos::as<magnitudeType>(pL.get<double>("rap: fix zero diagonals threshold"));
323 Scalar replacement = Teuchos::as<Scalar>(pL.get<double>("rap: fix zero diagonals replacement"));
324 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1), threshold, replacement);
325 }
326
327
328 if (IsPrint(Statistics2)) {
329 RCP<ParameterList> params = rcp(new ParameterList());;
330 params->set("printLoadBalancingInfo", true);
331 params->set("printCommInfo", true);
332 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
333 }
334
335 if(!Ac.is_null()) {std::ostringstream oss; oss << "A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
336 Set(coarseLevel, "A", Ac);
337
338 if (!isGPU) {
339 RAPparams->set("graph", Ac);
340 Set(coarseLevel, "RAP reuse data", RAPparams);
341 }
342 }
343
344
345 }
346
347#ifdef HAVE_MUELU_DEBUG
348 MatrixUtils::checkLocalRowMapMatchesColMap(*Ac);
349#endif // HAVE_MUELU_DEBUG
350
351 if (transferFacts_.begin() != transferFacts_.end()) {
352 SubFactoryMonitor m(*this, "Projections", coarseLevel);
353
354 // call Build of all user-given transfer factories
355 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
356 RCP<const FactoryBase> fac = *it;
357 GetOStream(Runtime0) << "RAPFactory: call transfer factory: " << fac->description() << std::endl;
358 fac->CallBuild(coarseLevel);
359 // Coordinates transfer is marginally different from all other operations
360 // because it is *optional*, and not required. For instance, we may need
361 // coordinates only on level 4 if we start repartitioning from that level,
362 // but we don't need them on level 1,2,3. As our current Hierarchy setup
363 // assumes propagation of dependencies only through three levels, this
364 // means that we need to rely on other methods to propagate optional data.
365 //
366 // The method currently used is through RAP transfer factories, which are
367 // simply factories which are called at the end of RAP with a single goal:
368 // transfer some fine data to coarser level. Because these factories are
369 // kind of outside of the mainline factories, they behave different. In
370 // particular, we call their Build method explicitly, rather than through
371 // Get calls. This difference is significant, as the Get call is smart
372 // enough to know when to release all factory dependencies, and Build is
373 // dumb. This led to the following CoordinatesTransferFactory sequence:
374 // 1. Request level 0
375 // 2. Request level 1
376 // 3. Request level 0
377 // 4. Release level 0
378 // 5. Release level 1
379 //
380 // The problem is missing "6. Release level 0". Because it was missing,
381 // we had outstanding request on "Coordinates", "Aggregates" and
382 // "CoarseMap" on level 0.
383 //
384 // This was fixed by explicitly calling Release on transfer factories in
385 // RAPFactory. I am still unsure how exactly it works, but now we have
386 // clear data requests for all levels.
387 coarseLevel.Release(*fac);
388 }
389 }
390
391 }
392
393 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
395 // check if it's a TwoLevelFactoryBase based transfer factory
396 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast,
397 "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. "
398 "This is very strange. (Note: you can remove this exception if there's a good reason for)");
399 TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_, Exceptions::RuntimeError, "MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
400 transferFacts_.push_back(factory);
401 }
402
403} //namespace MueLu
404
405#define MUELU_RAPFACTORY_SHORT
406#endif // MUELU_RAPFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
MueLu::DefaultScalar Scalar
MueLu::DefaultNode Node
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
int GetLevelID() const
Return level number.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Runtime0
One-liner description of what is happening.
@ Warnings1
Additional warnings.
static std::string getColonLabel(const std::string &label)
Helper function for object label.