MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_ParameterListInterpreter_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_PARAMETERLISTINTERPRETER_DEF_HPP
11 #define MUELU_PARAMETERLISTINTERPRETER_DEF_HPP
12 
14 
15 #include <Xpetra_Matrix.hpp>
16 #include <Xpetra_MatrixUtils.hpp>
17 
18 #include "MueLu_ConfigDefs.hpp"
19 
21 
22 #include "MueLu_MasterList.hpp"
23 #include "MueLu_Level.hpp"
24 #include "MueLu_Hierarchy.hpp"
25 #include "MueLu_FactoryManager.hpp"
26 
27 #include "MueLu_AggregationExportFactory.hpp"
28 #include "MueLu_AggregateQualityEstimateFactory.hpp"
29 #include "MueLu_AmalgamationFactory.hpp"
30 #include "MueLu_BrickAggregationFactory.hpp"
31 #include "MueLu_ClassicalMapFactory.hpp"
32 #include "MueLu_ClassicalPFactory.hpp"
33 #include "MueLu_CoalesceDropFactory.hpp"
34 #include "MueLu_CoarseMapFactory.hpp"
35 #include "MueLu_ConstraintFactory.hpp"
36 #include "MueLu_CoordinatesTransferFactory.hpp"
37 #include "MueLu_DirectSolver.hpp"
38 #include "MueLu_EminPFactory.hpp"
39 #include "MueLu_Exceptions.hpp"
40 #include "MueLu_FacadeClassFactory.hpp"
41 #include "MueLu_FactoryFactory.hpp"
42 #include "MueLu_FilteredAFactory.hpp"
43 #include "MueLu_GenericRFactory.hpp"
44 #include "MueLu_InitialBlockNumberFactory.hpp"
45 #include "MueLu_LineDetectionFactory.hpp"
46 #include "MueLu_LocalOrdinalTransferFactory.hpp"
47 #include "MueLu_MatrixAnalysisFactory.hpp"
48 #include "MueLu_MultiVectorTransferFactory.hpp"
49 #include "MueLu_NotayAggregationFactory.hpp"
50 #include "MueLu_NullspaceFactory.hpp"
51 #include "MueLu_PatternFactory.hpp"
52 #include "MueLu_ReplicatePFactory.hpp"
53 #include "MueLu_CombinePFactory.hpp"
54 #include "MueLu_PgPFactory.hpp"
55 #include "MueLu_RAPFactory.hpp"
56 #include "MueLu_RAPShiftFactory.hpp"
57 #include "MueLu_RebalanceAcFactory.hpp"
58 #include "MueLu_RebalanceTransferFactory.hpp"
59 #include "MueLu_RepartitionFactory.hpp"
60 #include "MueLu_RepartitionHeuristicFactory.hpp"
61 #include "MueLu_ReitzingerPFactory.hpp"
62 #include "MueLu_SaPFactory.hpp"
63 #include "MueLu_ScaledNullspaceFactory.hpp"
64 #include "MueLu_SemiCoarsenPFactory.hpp"
65 #include "MueLu_SmootherFactory.hpp"
66 #include "MueLu_SmooVecCoalesceDropFactory.hpp"
67 #include "MueLu_TentativePFactory.hpp"
68 #include "MueLu_TogglePFactory.hpp"
69 #include "MueLu_ToggleCoordinatesTransferFactory.hpp"
70 #include "MueLu_TransPFactory.hpp"
71 #include "MueLu_UncoupledAggregationFactory.hpp"
72 #include "MueLu_ZoltanInterface.hpp"
73 #include "MueLu_Zoltan2Interface.hpp"
74 #include "MueLu_NodePartitionInterface.hpp"
75 #include "MueLu_LowPrecisionFactory.hpp"
76 
77 #include "MueLu_CoalesceDropFactory_kokkos.hpp"
78 #include "MueLu_SemiCoarsenPFactory_kokkos.hpp"
79 #include "MueLu_TentativePFactory_kokkos.hpp"
80 
81 #ifdef HAVE_MUELU_MATLAB
82 #include "../matlab/src/MueLu_MatlabSmoother_decl.hpp"
83 #include "../matlab/src/MueLu_MatlabSmoother_def.hpp"
84 #include "../matlab/src/MueLu_TwoLevelMatlabFactory_decl.hpp"
85 #include "../matlab/src/MueLu_TwoLevelMatlabFactory_def.hpp"
86 #include "../matlab/src/MueLu_SingleLevelMatlabFactory_decl.hpp"
87 #include "../matlab/src/MueLu_SingleLevelMatlabFactory_def.hpp"
88 #endif
89 
90 #ifdef HAVE_MUELU_INTREPID2
91 #include "MueLu_IntrepidPCoarsenFactory.hpp"
92 #endif
93 
94 #include <unordered_set>
95 
96 namespace MueLu {
97 
98 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100  : factFact_(factFact) {
101  RCP<Teuchos::TimeMonitor> tM = rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(std::string("MueLu: ParameterListInterpreter (ParameterList)"))));
102  if (facadeFact == Teuchos::null)
104  else
105  facadeFact_ = facadeFact;
106 
107  if (paramList.isParameter("xml parameter file")) {
108  std::string filename = paramList.get("xml parameter file", "");
109  if (filename.length() != 0) {
110  TEUCHOS_TEST_FOR_EXCEPTION(comm.is_null(), Exceptions::RuntimeError, "xml parameter file requires a valid comm");
111 
112  ParameterList paramList2 = paramList;
113  Teuchos::updateParametersFromXmlFileAndBroadcast(filename, Teuchos::Ptr<Teuchos::ParameterList>(&paramList2), *comm);
114  SetParameterList(paramList2);
115 
116  } else {
117  SetParameterList(paramList);
118  }
119 
120  } else {
121  SetParameterList(paramList);
122  }
123 }
124 
125 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
127  : factFact_(factFact) {
128  RCP<Teuchos::TimeMonitor> tM = rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(std::string("MueLu: ParameterListInterpreter (XML)"))));
129  if (facadeFact == Teuchos::null)
131  else
132  facadeFact_ = facadeFact;
133 
134  ParameterList paramList;
135  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<ParameterList>(&paramList), comm);
136  SetParameterList(paramList);
137 }
138 
139 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141 
142 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
144  Cycle_ = Hierarchy::GetDefaultCycle();
145  WCycleStartLevel_ = Hierarchy::GetDefaultCycleStartLevel();
146  scalingFactor_ = Teuchos::ScalarTraits<double>::one();
147  blockSize_ = 1;
148  dofOffset_ = 0;
149 
150  if (paramList.isSublist("Hierarchy")) {
151  SetFactoryParameterList(paramList);
152 
153  } else if (paramList.isParameter("MueLu preconditioner") == true) {
154  this->GetOStream(Runtime0) << "Use facade class: " << paramList.get<std::string>("MueLu preconditioner") << std::endl;
155  Teuchos::RCP<ParameterList> pp = facadeFact_->SetParameterList(paramList);
156  SetFactoryParameterList(*pp);
157 
158  } else {
159  // The validator doesn't work correctly for non-serializable data (Hint: template parameters), so strip it out
160  ParameterList serialList, nonSerialList;
161 
162  ExtractNonSerializableData(paramList, serialList, nonSerialList);
163  Validate(serialList);
164  SetEasyParameterList(paramList);
165  }
166 }
167 
168 // =====================================================================================================
169 // ====================================== EASY interpreter =============================================
170 // =====================================================================================================
172 static inline bool areSame(const ParameterList& list1, const ParameterList& list2);
173 
174 // Get value from one of the lists, or set it to default
175 // Use case: check for a parameter value in a level-specific sublist, then in a root level list;
176 // if it is absent from both, set it to default
177 #define MUELU_SET_VAR_2LIST(paramList, defaultList, paramName, paramType, varName) \
178  paramType varName; \
179  if (paramList.isParameter(paramName)) \
180  varName = paramList.get<paramType>(paramName); \
181  else if (defaultList.isParameter(paramName)) \
182  varName = defaultList.get<paramType>(paramName); \
183  else \
184  varName = MasterList::getDefault<paramType>(paramName);
185 
186 #define MUELU_TEST_AND_SET_VAR(paramList, paramName, paramType, varName) \
187  (paramList.isParameter(paramName) ? varName = paramList.get<paramType>(paramName), true : false)
188 
189 // Set parameter in a list if it is present in any of two lists
190 // User case: set factory specific parameter, first checking for a level-specific value, then cheking root level value
191 #define MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, paramName, paramType, listWrite) \
192  try { \
193  if (paramList.isParameter(paramName)) \
194  listWrite.set(paramName, paramList.get<paramType>(paramName)); \
195  else if (defaultList.isParameter(paramName)) \
196  listWrite.set(paramName, defaultList.get<paramType>(paramName)); \
197  } catch (Teuchos::Exceptions::InvalidParameterType&) { \
198  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true, Teuchos::Exceptions::InvalidParameterType, \
199  "Error: parameter \"" << paramName << "\" must be of type " << Teuchos::TypeNameTraits<paramType>::name()); \
200  }
201 
202 #define MUELU_TEST_PARAM_2LIST(paramList, defaultList, paramName, paramType, cmpValue) \
203  (cmpValue == (paramList.isParameter(paramName) ? paramList.get<paramType>(paramName) : (defaultList.isParameter(paramName) ? defaultList.get<paramType>(paramName) : MasterList::getDefault<paramType>(paramName))))
204 
205 #define MUELU_KOKKOS_FACTORY(varName, oldFactory, newFactory) \
206  RCP<Factory> varName; \
207  if (!useKokkos_) \
208  varName = rcp(new oldFactory()); \
209  else \
210  varName = rcp(new newFactory());
211 #define MUELU_KOKKOS_FACTORY_NO_DECL(varName, oldFactory, newFactory) \
212  if (!useKokkos_) \
213  varName = rcp(new oldFactory()); \
214  else \
215  varName = rcp(new newFactory());
216 
217 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
219  SetEasyParameterList(const ParameterList& constParamList) {
220  ParameterList paramList;
221 
222  MUELU_SET_VAR_2LIST(constParamList, constParamList, "problem: type", std::string, problemType);
223  if (problemType != "unknown") {
224  paramList = *MasterList::GetProblemSpecificList(problemType);
225  paramList.setParameters(constParamList);
226  } else {
227  // Create a non const copy of the parameter list
228  // Working with a modifiable list is much much easier than with original one
229  paramList = constParamList;
230  }
231 
232  // Check for Kokkos
233  useKokkos_ = !Node::is_serial;
234  (void)MUELU_TEST_AND_SET_VAR(paramList, "use kokkos refactor", bool, useKokkos_);
235 
236  // Check for timer synchronization
237  MUELU_SET_VAR_2LIST(paramList, paramList, "synchronize factory timers", bool, syncTimers);
238  if (syncTimers)
240 
241  // Translate cycle type parameter
242  if (paramList.isParameter("cycle type")) {
243  std::map<std::string, CycleType> cycleMap;
244  cycleMap["V"] = VCYCLE;
245  cycleMap["W"] = WCYCLE;
246 
247  auto cycleType = paramList.get<std::string>("cycle type");
248  TEUCHOS_TEST_FOR_EXCEPTION(cycleMap.count(cycleType) == 0, Exceptions::RuntimeError,
249  "Invalid cycle type: \"" << cycleType << "\"");
250  Cycle_ = cycleMap[cycleType];
251  }
252 
253  if (paramList.isParameter("W cycle start level")) {
254  WCycleStartLevel_ = paramList.get<int>("W cycle start level");
255  }
256 
257  if (paramList.isParameter("coarse grid correction scaling factor"))
258  scalingFactor_ = paramList.get<double>("coarse grid correction scaling factor");
259 
260  this->maxCoarseSize_ = paramList.get<int>("coarse: max size", MasterList::getDefault<int>("coarse: max size"));
261  this->numDesiredLevel_ = paramList.get<int>("max levels", MasterList::getDefault<int>("max levels"));
262  blockSize_ = paramList.get<int>("number of equations", MasterList::getDefault<int>("number of equations"));
263 
264  (void)MUELU_TEST_AND_SET_VAR(paramList, "debug: graph level", int, this->graphOutputLevel_);
265 
266  // Generic data keeping (this keeps the data on all levels)
267  if (paramList.isParameter("keep data"))
268  this->dataToKeep_ = Teuchos::getArrayFromStringParameter<std::string>(paramList, "keep data");
269 
270  // Export level data
271  if (paramList.isSublist("export data")) {
272  ParameterList printList = paramList.sublist("export data");
273 
274  // Vectors, aggregates and other things that need special handling
275  if (printList.isParameter("Nullspace"))
276  this->nullspaceToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "Nullspace");
277  if (printList.isParameter("Coordinates"))
278  this->coordinatesToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "Coordinates");
279  if (printList.isParameter("Material"))
280  this->materialToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "Material");
281  if (printList.isParameter("Aggregates"))
282  this->aggregatesToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "Aggregates");
283  if (printList.isParameter("pcoarsen: element to node map"))
284  this->elementToNodeMapsToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "pcoarsen: element to node map");
285 
286  // If we asked for an arbitrary matrix to be printed, we do that here
287  for (auto iter = printList.begin(); iter != printList.end(); iter++) {
288  const std::string& name = printList.name(iter);
289  // Ignore the special cases
290  if (name == "Nullspace" || name == "Coordinates" || name == "Material" || name == "Aggregates" || name == "pcoarsen: element to node map")
291  continue;
292 
293  this->matricesToPrint_[name] = Teuchos::getArrayFromStringParameter<int>(printList, name);
294  }
295  }
296 
297  // Set verbosity parameter
299  {
300  MUELU_SET_VAR_2LIST(paramList, paramList, "verbosity", std::string, verbosityLevel);
301  this->verbosity_ = toVerbLevel(verbosityLevel);
302  VerboseObject::SetDefaultVerbLevel(this->verbosity_);
303  }
304 
305  MUELU_SET_VAR_2LIST(paramList, paramList, "output filename", std::string, outputFilename);
306  if (outputFilename != "")
307  VerboseObject::SetMueLuOFileStream(outputFilename);
308 
309  // Detect if we need to transfer coordinates to coarse levels. We do that iff
310  // - we use "distance laplacian" dropping on some level, or
311  // - we use a repartitioner on some level that needs coordinates
312  // - we use brick aggregation
313  // - we use Ifpack2 line partitioner
314  // This is not ideal, as we may have "repartition: enable" turned on by default
315  // and not present in the list, but it is better than nothing.
316  useCoordinates_ = false;
317  useBlockNumber_ = false;
318  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: strength-of-connection: matrix", std::string, "distance laplacian"))
319  useCoordinates_ = true;
320  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: use blocking", bool, true))
321  useBlockNumber_ = true;
322  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "distance laplacian") ||
323  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: type", std::string, "brick") ||
324  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: export visualization data", bool, true)) {
325  useCoordinates_ = true;
326  } else if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "block diagonal distance laplacian")) {
327  useCoordinates_ = true;
328  useBlockNumber_ = true;
329  } else if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "block diagonal") ||
330  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "block diagonal classical") ||
331  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "block diagonal signed classical") ||
332  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "block diagonal colored signed classical") ||
333  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "signed classical")) {
334  useBlockNumber_ = true;
335  } else if (paramList.isSublist("smoother: params")) {
336  const auto smooParamList = paramList.sublist("smoother: params");
337  if (smooParamList.isParameter("partitioner: type") &&
338  (smooParamList.get<std::string>("partitioner: type") == "line")) {
339  useCoordinates_ = true;
340  }
341  } else {
342  for (int levelID = 0; levelID < this->numDesiredLevel_; levelID++) {
343  std::string levelStr = "level " + toString(levelID);
344 
345  if (paramList.isSublist(levelStr)) {
346  const ParameterList& levelList = paramList.sublist(levelStr);
347 
348  if (MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: drop scheme", std::string, "distance laplacian") ||
349  MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: type", std::string, "brick") ||
350  MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: export visualization data", bool, true)) {
351  useCoordinates_ = true;
352  } else if (MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: drop scheme", std::string, "block diagonal distance laplacian")) {
353  useCoordinates_ = true;
354  useBlockNumber_ = true;
355  } else if (MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: drop scheme", std::string, "block diagonal") ||
356  MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: drop scheme", std::string, "block diagonal classical") ||
357  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "block diagonal signed classical") ||
358  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "block diagonal colored signed classical") ||
359  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "signed classical")) {
360  useBlockNumber_ = true;
361  }
362  }
363  }
364  }
365 
366  useMaterial_ = false;
367  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: distance laplacian metric", std::string, "material")) {
368  useMaterial_ = true;
369  }
370 
371  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "repartition: enable", bool, true)) {
372  // We don't need coordinates if we're doing the in-place restriction
373  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "repartition: use subcommunicators", bool, true) &&
374  MUELU_TEST_PARAM_2LIST(paramList, paramList, "repartition: use subcommunicators in place", bool, true)) {
375  // do nothing --- these don't need coordinates
376  } else if (!paramList.isSublist("repartition: params")) {
377  useCoordinates_ = true;
378  } else {
379  const ParameterList& repParams = paramList.sublist("repartition: params");
380  if (repParams.isType<std::string>("algorithm")) {
381  const std::string algo = repParams.get<std::string>("algorithm");
382  if (algo == "multijagged" || algo == "rcb") {
383  useCoordinates_ = true;
384  }
385  } else {
386  useCoordinates_ = true;
387  }
388  }
389  }
390  for (int levelID = 0; levelID < this->numDesiredLevel_; levelID++) {
391  std::string levelStr = "level " + toString(levelID);
392 
393  if (paramList.isSublist(levelStr)) {
394  const ParameterList& levelList = paramList.sublist(levelStr);
395 
396  if (MUELU_TEST_PARAM_2LIST(levelList, paramList, "repartition: enable", bool, true)) {
397  if (!levelList.isSublist("repartition: params")) {
398  useCoordinates_ = true;
399  break;
400  } else {
401  const ParameterList& repParams = levelList.sublist("repartition: params");
402  if (repParams.isType<std::string>("algorithm")) {
403  const std::string algo = repParams.get<std::string>("algorithm");
404  if (algo == "multijagged" || algo == "rcb") {
405  useCoordinates_ = true;
406  break;
407  }
408  } else {
409  useCoordinates_ = true;
410  break;
411  }
412  }
413  }
414  }
415  }
416 
417  // Detect if we do implicit P and R rebalance
418  changedPRrebalance_ = false;
419  changedPRViaCopyrebalance_ = false;
420  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "repartition: enable", bool, true)) {
421  changedPRrebalance_ = MUELU_TEST_AND_SET_VAR(paramList, "repartition: rebalance P and R", bool, this->doPRrebalance_);
422  changedPRViaCopyrebalance_ = MUELU_TEST_AND_SET_VAR(paramList, "repartition: explicit via new copy rebalance P and R", bool, this->doPRViaCopyrebalance_);
423  }
424 
425  // Detect if we use implicit transpose
426  changedImplicitTranspose_ = MUELU_TEST_AND_SET_VAR(paramList, "transpose: use implicit", bool, this->implicitTranspose_);
427 
428  // Detect if we use fuse prolongation and update
429  (void)MUELU_TEST_AND_SET_VAR(paramList, "fuse prolongation and update", bool, this->fuseProlongationAndUpdate_);
430 
431  // Detect if we suppress the dimension check of the user-given nullspace
432  (void)MUELU_TEST_AND_SET_VAR(paramList, "nullspace: suppress dimension check", bool, this->suppressNullspaceDimensionCheck_);
433 
434  if (paramList.isSublist("matvec params"))
435  this->matvecParams_ = Teuchos::parameterList(paramList.sublist("matvec params"));
436 
437  // Create default manager
438  // FIXME: should it be here, or higher up
439  RCP<FactoryManager> defaultManager = rcp(new FactoryManager());
440  defaultManager->SetVerbLevel(this->verbosity_);
441  defaultManager->SetKokkosRefactor(useKokkos_);
442 
443  // We will ignore keeps0
444  std::vector<keep_pair> keeps0;
445  UpdateFactoryManager(paramList, ParameterList(), *defaultManager, 0 /*levelID*/, keeps0);
446 
447  // std::cout<<"*** Default Manager ***"<<std::endl;
448  // defaultManager->Print();
449 
450  // Create level specific factory managers
451  for (int levelID = 0; levelID < this->numDesiredLevel_; levelID++) {
452  // Note, that originally if there were no level specific parameters, we
453  // simply copied the defaultManager However, with the introduction of
454  // levelID to UpdateFactoryManager (required for reuse), we can no longer
455  // guarantee that the kept variables are the same for each level even if
456  // dependency structure does not change.
457  RCP<FactoryManager> levelManager = rcp(new FactoryManager(*defaultManager));
458  levelManager->SetVerbLevel(defaultManager->GetVerbLevel());
459 
460  std::vector<keep_pair> keeps;
461  if (paramList.isSublist("level " + toString(levelID))) {
462  // We do this so the parameters on the level get flagged correctly as "used"
463  ParameterList& levelList = paramList.sublist("level " + toString(levelID), true /*mustAlreadyExist*/);
464  UpdateFactoryManager(levelList, paramList, *levelManager, levelID, keeps);
465 
466  } else {
467  ParameterList levelList;
468  UpdateFactoryManager(levelList, paramList, *levelManager, levelID, keeps);
469  }
470 
471  this->keep_[levelID] = keeps;
472  this->AddFactoryManager(levelID, 1, levelManager);
473 
474  // std::cout<<"*** Level "<<levelID<<" Manager ***"<<std::endl;
475  // levelManager->Print();
476  }
477 
478  // FIXME: parameters passed to packages, like Ifpack2, are not touched by us, resulting in "[unused]" flag
479  // being displayed. On the other hand, we don't want to simply iterate through them touching. I don't know
480  // what a good solution looks like
481  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "print initial parameters", bool, true))
482  this->GetOStream(static_cast<MsgType>(Runtime1), 0) << paramList << std::endl;
483 
484  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "print unused parameters", bool, true)) {
485  // Check unused parameters
486  ParameterList unusedParamList;
487 
488  // Check for unused parameters that aren't lists
489  for (ParameterList::ConstIterator it = paramList.begin(); it != paramList.end(); it++) {
490  const ParameterEntry& entry = paramList.entry(it);
491 
492  if (!entry.isList() && !entry.isUsed())
493  unusedParamList.setEntry(paramList.name(it), entry);
494  }
495 
496  // Check for unused parameters in level-specific sublists
497  for (int levelID = 0; levelID < this->numDesiredLevel_; levelID++) {
498  std::string levelStr = "level " + toString(levelID);
499 
500  if (paramList.isSublist(levelStr)) {
501  const ParameterList& levelList = paramList.sublist(levelStr);
502 
503  for (ParameterList::ConstIterator itr = levelList.begin(); itr != levelList.end(); ++itr) {
504  const ParameterEntry& entry = levelList.entry(itr);
505 
506  if (!entry.isList() && !entry.isUsed())
507  unusedParamList.sublist(levelStr).setEntry(levelList.name(itr), entry);
508  }
509  }
510  }
511 
512  if (unusedParamList.numParams() > 0) {
513  std::ostringstream unusedParamsStream;
514  int indent = 4;
515  unusedParamList.print(unusedParamsStream, indent);
516 
517  this->GetOStream(Warnings1) << "The following parameters were not used:\n"
518  << unusedParamsStream.str() << std::endl;
519  }
520  }
521 
523 }
524 
525 // =====================================================================================================
526 // ==================================== UpdateFactoryManager ===========================================
527 // =====================================================================================================
528 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
530  UpdateFactoryManager(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager,
531  int levelID, std::vector<keep_pair>& keeps) const {
532  // NOTE: Factory::SetParameterList must be called prior to Factory::SetFactory, as
533  // SetParameterList sets default values for non mentioned parameters, including factories
534 
535  using strings = std::unordered_set<std::string>;
536 
537  // shortcut
538  if (paramList.numParams() == 0 && defaultList.numParams() > 0)
539  paramList = ParameterList(defaultList);
540 
541  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
542  TEUCHOS_TEST_FOR_EXCEPTION(strings({"none", "tP", "RP", "emin", "RAP", "full", "S"}).count(reuseType) == 0,
543  Exceptions::RuntimeError, "Unknown \"reuse: type\" value: \"" << reuseType << "\". Please consult User's Guide.");
544 
545  MUELU_SET_VAR_2LIST(paramList, defaultList, "multigrid algorithm", std::string, multigridAlgo);
546  TEUCHOS_TEST_FOR_EXCEPTION(strings({"unsmoothed", "sa", "pg", "emin", "matlab", "pcoarsen", "classical", "smoothed reitzinger", "unsmoothed reitzinger", "replicate", "combine"}).count(multigridAlgo) == 0,
547  Exceptions::RuntimeError, "Unknown \"multigrid algorithm\" value: \"" << multigridAlgo << "\". Please consult User's Guide.");
548 #ifndef HAVE_MUELU_MATLAB
549  TEUCHOS_TEST_FOR_EXCEPTION(multigridAlgo == "matlab", Exceptions::RuntimeError,
550  "Cannot use matlab for multigrid algorithm - MueLu was not configured with MATLAB support.");
551 #endif
552 #ifndef HAVE_MUELU_INTREPID2
553  TEUCHOS_TEST_FOR_EXCEPTION(multigridAlgo == "pcoarsen", Exceptions::RuntimeError,
554  "Cannot use IntrepidPCoarsen prolongator factory - MueLu was not configured with Intrepid support.");
555 #endif
556 
557  // Only some combinations of reuse and multigrid algorithms are tested, all
558  // other are considered invalid at the moment
559  if (reuseType == "none" || reuseType == "S" || reuseType == "RP" || reuseType == "RAP") {
560  // This works for all kinds of multigrid algorithms
561 
562  } else if (reuseType == "tP" && (multigridAlgo != "sa" && multigridAlgo != "unsmoothed")) {
563  reuseType = "none";
564  this->GetOStream(Warnings0) << "Ignoring \"tP\" reuse option as it is only compatible with \"sa\", "
565  "or \"unsmoothed\" multigrid algorithms"
566  << std::endl;
567 
568  } else if (reuseType == "emin" && multigridAlgo != "emin") {
569  reuseType = "none";
570  this->GetOStream(Warnings0) << "Ignoring \"emin\" reuse option it is only compatible with "
571  "\"emin\" multigrid algorithm"
572  << std::endl;
573  }
574 
575  // == Non-serializable data ===
576  // Check both the parameter and the type
577  bool have_userP = false;
578  if (paramList.isParameter("P") && !paramList.get<RCP<Matrix> >("P").is_null())
579  have_userP = true;
580 
581  // === Coarse solver ===
582  UpdateFactoryManager_CoarseSolvers(paramList, defaultList, manager, levelID, keeps);
583 
584  // == Smoothers ==
585  UpdateFactoryManager_Smoothers(paramList, defaultList, manager, levelID, keeps);
586 
587  // === BlockNumber ===
588  if (levelID == 0)
589  UpdateFactoryManager_BlockNumber(paramList, defaultList, manager, levelID, keeps);
590 
591  // === Aggregation ===
592  if (multigridAlgo == "unsmoothed reitzinger" || multigridAlgo == "smoothed reitzinger")
593  UpdateFactoryManager_Reitzinger(paramList, defaultList, manager, levelID, keeps);
594  else
595  UpdateFactoryManager_Aggregation_TentativeP(paramList, defaultList, manager, levelID, keeps);
596 
597  // === Nullspace ===
598  RCP<Factory> nullSpaceFactory; // Cache thcAN is guy for the combination of semi-coarsening & repartitioning
599  UpdateFactoryManager_Nullspace(paramList, defaultList, manager, levelID, keeps, nullSpaceFactory);
600 
601  // === Prolongation ===
602  // NOTE: None of the UpdateFactoryManager routines called here check the
603  // multigridAlgo. This is intentional, to allow for reuse of components
604  // underneath. Thus, the multigridAlgo was checked in the beginning of the
605  // function.
606  if (have_userP) {
607  // User prolongator
608  manager.SetFactory("P", NoFactory::getRCP());
609 
610  } else if (multigridAlgo == "unsmoothed" || multigridAlgo == "unsmoothed reitzinger") {
611  // Unsmoothed aggregation
612  manager.SetFactory("P", manager.GetFactory("Ptent"));
613 
614  } else if (multigridAlgo == "classical") {
615  // Classical AMG
616  manager.SetFactory("P", manager.GetFactory("Ptent"));
617 
618  } else if (multigridAlgo == "sa" || multigridAlgo == "smoothed reitzinger") {
619  // Smoothed aggregation
620  UpdateFactoryManager_SA(paramList, defaultList, manager, levelID, keeps);
621 
622  } else if (multigridAlgo == "emin") {
623  // Energy minimization
624  UpdateFactoryManager_Emin(paramList, defaultList, manager, levelID, keeps);
625 
626  } else if (multigridAlgo == "replicate") {
627  UpdateFactoryManager_Replicate(paramList, defaultList, manager, levelID, keeps);
628 
629  } else if (multigridAlgo == "combine") {
630  UpdateFactoryManager_Combine(paramList, defaultList, manager, levelID, keeps);
631 
632  } else if (multigridAlgo == "pg") {
633  // Petrov-Galerkin
634  UpdateFactoryManager_PG(paramList, defaultList, manager, levelID, keeps);
635 
636  } else if (multigridAlgo == "matlab") {
637  // Matlab Coarsneing
638  UpdateFactoryManager_Matlab(paramList, defaultList, manager, levelID, keeps);
639 
640  } else if (multigridAlgo == "pcoarsen") {
641  // P-Coarsening
642  UpdateFactoryManager_PCoarsen(paramList, defaultList, manager, levelID, keeps);
643  }
644 
645  // === Semi-coarsening ===
646  UpdateFactoryManager_SemiCoarsen(paramList, defaultList, manager, levelID, keeps);
647 
648  // === Restriction ===
649  UpdateFactoryManager_Restriction(paramList, defaultList, manager, levelID, keeps);
650 
651  // === RAP ===
652  UpdateFactoryManager_RAP(paramList, defaultList, manager, levelID, keeps);
653 
654  // == BlockNumber Transfer ==
655  UpdateFactoryManager_LocalOrdinalTransfer("BlockNumber", multigridAlgo, paramList, defaultList, manager, levelID, keeps);
656 
657  // === Coordinates ===
658  UpdateFactoryManager_Coordinates(paramList, defaultList, manager, levelID, keeps);
659 
660  // === Material ===
661  UpdateFactoryManager_Material(paramList, defaultList, manager, levelID, keeps);
662 
663  // === Pre-Repartition Keeps for Reuse ===
664  if ((reuseType == "RP" || reuseType == "RAP" || reuseType == "full") && levelID)
665  keeps.push_back(keep_pair("Nullspace", manager.GetFactory("Nullspace").get()));
666 
667  if (reuseType == "RP" && levelID) {
668  keeps.push_back(keep_pair("P", manager.GetFactory("P").get()));
669  if (!this->implicitTranspose_)
670  keeps.push_back(keep_pair("R", manager.GetFactory("R").get()));
671  }
672  if ((reuseType == "tP" || reuseType == "RP" || reuseType == "emin") && useCoordinates_ && levelID)
673  keeps.push_back(keep_pair("Coordinates", manager.GetFactory("Coordinates").get()));
674 
675  // === Repartitioning ===
676  UpdateFactoryManager_Repartition(paramList, defaultList, manager, levelID, keeps, nullSpaceFactory);
677 
678  // === Lower precision transfers ===
679  UpdateFactoryManager_LowPrecision(paramList, defaultList, manager, levelID, keeps);
680 
681  // === Final Keeps for Reuse ===
682  if ((reuseType == "RAP" || reuseType == "full") && levelID) {
683  keeps.push_back(keep_pair("P", manager.GetFactory("P").get()));
684  if (!this->implicitTranspose_)
685  keeps.push_back(keep_pair("R", manager.GetFactory("R").get()));
686  keeps.push_back(keep_pair("A", manager.GetFactory("A").get()));
687  }
688 
689  // In case you ever want to inspect the FactoryManager as it is generated for each level
690  /*std::cout<<"*** Factory Manager on level "<<levelID<<" ***"<<std::endl;
691  manager.Print(); */
692 }
693 
694 // =====================================================================================================
695 // ========================================= Smoothers =================================================
696 // =====================================================================================================
697 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
700  FactoryManager& manager, int levelID, std::vector<keep_pair>& keeps) const {
701  MUELU_SET_VAR_2LIST(paramList, defaultList, "multigrid algorithm", std::string, multigridAlgo);
702  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
703  bool useMaxAbsDiagonalScaling = false;
704  if (defaultList.isParameter("sa: use rowsumabs diagonal scaling"))
705  useMaxAbsDiagonalScaling = defaultList.get<bool>("sa: use rowsumabs diagonal scaling");
706 
707  // === Smoothing ===
708  // FIXME: should custom smoother check default list too?
709  bool isCustomSmoother =
710  paramList.isParameter("smoother: pre or post") ||
711  paramList.isParameter("smoother: type") || paramList.isParameter("smoother: pre type") || paramList.isParameter("smoother: post type") ||
712  paramList.isSublist("smoother: params") || paramList.isSublist("smoother: pre params") || paramList.isSublist("smoother: post params") ||
713  paramList.isParameter("smoother: sweeps") || paramList.isParameter("smoother: pre sweeps") || paramList.isParameter("smoother: post sweeps") ||
714  paramList.isParameter("smoother: overlap") || paramList.isParameter("smoother: pre overlap") || paramList.isParameter("smoother: post overlap");
715 
716  MUELU_SET_VAR_2LIST(paramList, defaultList, "smoother: pre or post", std::string, PreOrPost);
717  if (PreOrPost == "none") {
718  manager.SetFactory("Smoother", Teuchos::null);
719 
720  } else if (isCustomSmoother) {
721  // FIXME: get default values from the factory
722  // NOTE: none of the smoothers at the moment use parameter validation framework, so we
723  // cannot get the default values from it.
724 #define TEST_MUTUALLY_EXCLUSIVE(arg1, arg2) \
725  TEUCHOS_TEST_FOR_EXCEPTION(paramList.isParameter(#arg1) && paramList.isParameter(#arg2), \
726  Exceptions::InvalidArgument, "You cannot specify both \"" #arg1 "\" and \"" #arg2 "\"");
727 #define TEST_MUTUALLY_EXCLUSIVE_S(arg1, arg2) \
728  TEUCHOS_TEST_FOR_EXCEPTION(paramList.isSublist(#arg1) && paramList.isSublist(#arg2), \
729  Exceptions::InvalidArgument, "You cannot specify both \"" #arg1 "\" and \"" #arg2 "\"");
730 
731  TEST_MUTUALLY_EXCLUSIVE("smoother: type", "smoother: pre type");
732  TEST_MUTUALLY_EXCLUSIVE("smoother: type", "smoother: post type");
733  TEST_MUTUALLY_EXCLUSIVE("smoother: sweeps", "smoother: pre sweeps");
734  TEST_MUTUALLY_EXCLUSIVE("smoother: sweeps", "smoother: post sweeps");
735  TEST_MUTUALLY_EXCLUSIVE("smoother: overlap", "smoother: pre overlap");
736  TEST_MUTUALLY_EXCLUSIVE("smoother: overlap", "smoother: post overlap");
737  TEST_MUTUALLY_EXCLUSIVE_S("smoother: params", "smoother: pre params");
738  TEST_MUTUALLY_EXCLUSIVE_S("smoother: params", "smoother: post params");
739  TEUCHOS_TEST_FOR_EXCEPTION(PreOrPost == "both" && (paramList.isParameter("smoother: pre type") != paramList.isParameter("smoother: post type")),
740  Exceptions::InvalidArgument, "You must specify both \"smoother: pre type\" and \"smoother: post type\"");
741 
742  // Default values
743  int overlap = 0;
744  ParameterList defaultSmootherParams;
745  defaultSmootherParams.set("relaxation: type", "Symmetric Gauss-Seidel");
746  defaultSmootherParams.set("relaxation: sweeps", Teuchos::OrdinalTraits<LO>::one());
747  defaultSmootherParams.set("relaxation: damping factor", Teuchos::ScalarTraits<Scalar>::one());
748 
749  RCP<SmootherFactory> preSmoother = Teuchos::null, postSmoother = Teuchos::null;
750  std::string preSmootherType, postSmootherType;
751  ParameterList preSmootherParams, postSmootherParams;
752 
753  if (paramList.isParameter("smoother: overlap"))
754  overlap = paramList.get<int>("smoother: overlap");
755 
756  if (PreOrPost == "pre" || PreOrPost == "both") {
757  if (paramList.isParameter("smoother: pre type")) {
758  preSmootherType = paramList.get<std::string>("smoother: pre type");
759  } else {
760  MUELU_SET_VAR_2LIST(paramList, defaultList, "smoother: type", std::string, preSmootherTypeTmp);
761  preSmootherType = preSmootherTypeTmp;
762  }
763  if (paramList.isParameter("smoother: pre overlap"))
764  overlap = paramList.get<int>("smoother: pre overlap");
765 
766  if (paramList.isSublist("smoother: pre params"))
767  preSmootherParams = paramList.sublist("smoother: pre params");
768  else if (paramList.isSublist("smoother: params"))
769  preSmootherParams = paramList.sublist("smoother: params");
770  else if (defaultList.isSublist("smoother: params"))
771  preSmootherParams = defaultList.sublist("smoother: params");
772  else if (preSmootherType == "RELAXATION")
773  preSmootherParams = defaultSmootherParams;
774 
775  if (preSmootherType == "CHEBYSHEV" && useMaxAbsDiagonalScaling)
776  preSmootherParams.set("chebyshev: use rowsumabs diagonal scaling", true);
777 
778 #ifdef HAVE_MUELU_INTREPID2
779  // Propagate P-coarsening for Topo smoothing
780  if (multigridAlgo == "pcoarsen" && preSmootherType == "TOPOLOGICAL" &&
781  defaultList.isParameter("pcoarsen: schedule") && defaultList.isParameter("pcoarsen: element")) {
782  // P-Coarsening by schedule (new interface)
783  // NOTE: levelID represents the *coarse* level in this case
784  auto pcoarsen_schedule = Teuchos::getArrayFromStringParameter<int>(defaultList, "pcoarsen: schedule");
785  auto pcoarsen_element = defaultList.get<std::string>("pcoarsen: element");
786 
787  if (levelID < (int)pcoarsen_schedule.size()) {
788  // Topo info for P-Coarsening
789  auto lo = pcoarsen_element + std::to_string(pcoarsen_schedule[levelID]);
790  preSmootherParams.set("pcoarsen: hi basis", lo);
791  }
792  }
793 #endif
794 
795 #ifdef HAVE_MUELU_MATLAB
796  if (preSmootherType == "matlab")
797  preSmoother = rcp(new SmootherFactory(rcp(new MatlabSmoother(preSmootherParams))));
798  else
799 #endif
800  preSmoother = rcp(new SmootherFactory(rcp(new TrilinosSmoother(preSmootherType, preSmootherParams, overlap))));
801  }
802 
803  if (PreOrPost == "post" || PreOrPost == "both") {
804  if (paramList.isParameter("smoother: post type"))
805  postSmootherType = paramList.get<std::string>("smoother: post type");
806  else {
807  MUELU_SET_VAR_2LIST(paramList, defaultList, "smoother: type", std::string, postSmootherTypeTmp);
808  postSmootherType = postSmootherTypeTmp;
809  }
810 
811  if (paramList.isSublist("smoother: post params"))
812  postSmootherParams = paramList.sublist("smoother: post params");
813  else if (paramList.isSublist("smoother: params"))
814  postSmootherParams = paramList.sublist("smoother: params");
815  else if (defaultList.isSublist("smoother: params"))
816  postSmootherParams = defaultList.sublist("smoother: params");
817  else if (postSmootherType == "RELAXATION")
818  postSmootherParams = defaultSmootherParams;
819  if (paramList.isParameter("smoother: post overlap"))
820  overlap = paramList.get<int>("smoother: post overlap");
821 
822  if (postSmootherType == "CHEBYSHEV" && useMaxAbsDiagonalScaling)
823  postSmootherParams.set("chebyshev: use rowsumabs diagonal scaling", true);
824 
825  if (postSmootherType == preSmootherType && areSame(preSmootherParams, postSmootherParams))
826  postSmoother = preSmoother;
827  else {
828 #ifdef HAVE_MUELU_INTREPID2
829  // Propagate P-coarsening for Topo smoothing
830  if (multigridAlgo == "pcoarsen" && preSmootherType == "TOPOLOGICAL" &&
831  defaultList.isParameter("pcoarsen: schedule") && defaultList.isParameter("pcoarsen: element")) {
832  // P-Coarsening by schedule (new interface)
833  // NOTE: levelID represents the *coarse* level in this case
834  auto pcoarsen_schedule = Teuchos::getArrayFromStringParameter<int>(defaultList, "pcoarsen: schedule");
835  auto pcoarsen_element = defaultList.get<std::string>("pcoarsen: element");
836 
837  if (levelID < (int)pcoarsen_schedule.size()) {
838  // Topo info for P-Coarsening
839  auto lo = pcoarsen_element + std::to_string(pcoarsen_schedule[levelID]);
840  postSmootherParams.set("pcoarsen: hi basis", lo);
841  }
842  }
843 #endif
844 
845 #ifdef HAVE_MUELU_MATLAB
846  if (postSmootherType == "matlab")
847  postSmoother = rcp(new SmootherFactory(rcp(new MatlabSmoother(postSmootherParams))));
848  else
849 #endif
850  postSmoother = rcp(new SmootherFactory(rcp(new TrilinosSmoother(postSmootherType, postSmootherParams, overlap))));
851  }
852  }
853 
854  if (preSmoother == postSmoother)
855  manager.SetFactory("Smoother", preSmoother);
856  else {
857  manager.SetFactory("PreSmoother", preSmoother);
858  manager.SetFactory("PostSmoother", postSmoother);
859  }
860  }
861 
862  // The first clause is not necessary, but it is here for clarity Smoothers
863  // are reused if smoother explicitly said to reuse them, or if any other
864  // reuse option is enabled
865  bool reuseSmoothers = (reuseType == "S" || reuseType != "none");
866  if (reuseSmoothers) {
867  auto preSmootherFactory = rcp_const_cast<Factory>(rcp_dynamic_cast<const Factory>(manager.GetFactory("PreSmoother")));
868 
869  if (preSmootherFactory != Teuchos::null) {
870  ParameterList postSmootherFactoryParams;
871  postSmootherFactoryParams.set("keep smoother data", true);
872  preSmootherFactory->SetParameterList(postSmootherFactoryParams);
873 
874  keeps.push_back(keep_pair("PreSmoother data", preSmootherFactory.get()));
875  }
876 
877  auto postSmootherFactory = rcp_const_cast<Factory>(rcp_dynamic_cast<const Factory>(manager.GetFactory("PostSmoother")));
878  if (postSmootherFactory != Teuchos::null) {
879  ParameterList postSmootherFactoryParams;
880  postSmootherFactoryParams.set("keep smoother data", true);
881  postSmootherFactory->SetParameterList(postSmootherFactoryParams);
882 
883  keeps.push_back(keep_pair("PostSmoother data", postSmootherFactory.get()));
884  }
885 
886  auto coarseFactory = rcp_const_cast<Factory>(rcp_dynamic_cast<const Factory>(manager.GetFactory("CoarseSolver")));
887  if (coarseFactory != Teuchos::null) {
888  ParameterList coarseFactoryParams;
889  coarseFactoryParams.set("keep smoother data", true);
890  coarseFactory->SetParameterList(coarseFactoryParams);
891 
892  keeps.push_back(keep_pair("PreSmoother data", coarseFactory.get()));
893  }
894  }
895 
896  if ((reuseType == "RAP" && levelID) || (reuseType == "full")) {
897  // The difference between "RAP" and "full" is keeping smoothers. However,
898  // as in both cases we keep coarse matrices, we do not need to update
899  // coarse smoothers. On the other hand, if a user changes fine level
900  // matrix, "RAP" would update the fine level smoother, while "full" would
901  // not
902  keeps.push_back(keep_pair("PreSmoother", manager.GetFactory("PreSmoother").get()));
903  keeps.push_back(keep_pair("PostSmoother", manager.GetFactory("PostSmoother").get()));
904 
905  // We do keep_pair("PreSmoother", manager.GetFactory("CoarseSolver").get())
906  // as the coarse solver factory is in fact a smoothing factory, so the
907  // only pieces of data it generates are PreSmoother and PostSmoother
908  keeps.push_back(keep_pair("PreSmoother", manager.GetFactory("CoarseSolver").get()));
909  }
910 }
911 
912 // =====================================================================================================
913 // ====================================== Coarse Solvers ===============================================
914 // =====================================================================================================
915 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
918  FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
919  // FIXME: should custom coarse solver check default list too?
920  bool isCustomCoarseSolver =
921  paramList.isParameter("coarse: type") ||
922  paramList.isParameter("coarse: params");
923  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "coarse: type", std::string, "none")) {
924  manager.SetFactory("CoarseSolver", Teuchos::null);
925 
926  } else if (isCustomCoarseSolver) {
927  // FIXME: get default values from the factory
928  // NOTE: none of the smoothers at the moment use parameter validation framework, so we
929  // cannot get the default values from it.
930  MUELU_SET_VAR_2LIST(paramList, defaultList, "coarse: type", std::string, coarseType);
931 
932  int overlap = 0;
933  if (paramList.isParameter("coarse: overlap"))
934  overlap = paramList.get<int>("coarse: overlap");
935 
936  ParameterList coarseParams;
937  if (paramList.isSublist("coarse: params"))
938  coarseParams = paramList.sublist("coarse: params");
939  else if (defaultList.isSublist("coarse: params"))
940  coarseParams = defaultList.sublist("coarse: params");
941 
942  using strings = std::unordered_set<std::string>;
943 
944  RCP<SmootherPrototype> coarseSmoother;
945  // TODO: this is not a proper place to check. If we consider direct solver to be a special
946  // case of smoother, we would like to unify Amesos and Ifpack2 smoothers in src/Smoothers, and
947  // have a single factory responsible for those. Then, this check would belong there.
948  if (strings({"RELAXATION", "CHEBYSHEV", "ILUT", "ILU", "RILUK", "SCHWARZ", "Amesos",
949  "BLOCK RELAXATION", "BLOCK_RELAXATION", "BLOCKRELAXATION",
950  "SPARSE BLOCK RELAXATION", "SPARSE_BLOCK_RELAXATION", "SPARSEBLOCKRELAXATION",
951  "LINESMOOTHING_BANDEDRELAXATION", "LINESMOOTHING_BANDED_RELAXATION", "LINESMOOTHING_BANDED RELAXATION",
952  "LINESMOOTHING_TRIDIRELAXATION", "LINESMOOTHING_TRIDI_RELAXATION", "LINESMOOTHING_TRIDI RELAXATION",
953  "LINESMOOTHING_TRIDIAGONALRELAXATION", "LINESMOOTHING_TRIDIAGONAL_RELAXATION", "LINESMOOTHING_TRIDIAGONAL RELAXATION",
954  "TOPOLOGICAL", "FAST_ILU", "FAST_IC", "FAST_ILDL", "HIPTMAIR"})
955  .count(coarseType)) {
956  coarseSmoother = rcp(new TrilinosSmoother(coarseType, coarseParams, overlap));
957  } else {
958 #ifdef HAVE_MUELU_MATLAB
959  if (coarseType == "matlab")
960  coarseSmoother = rcp(new MatlabSmoother(coarseParams));
961  else
962 #endif
963  coarseSmoother = rcp(new DirectSolver(coarseType, coarseParams));
964  }
965 
966  manager.SetFactory("CoarseSolver", rcp(new SmootherFactory(coarseSmoother)));
967  }
968 }
969 
970 // =====================================================================================================
971 // ========================================= TentativeP=================================================
972 // =====================================================================================================
973 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
976  FactoryManager& manager, int levelID, std::vector<keep_pair>& keeps) const {
977  ParameterList rParams;
978  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: enable", bool, rParams);
979  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, rParams);
980  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: constant column sums", bool, rParams);
981  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: calculate qr", bool, rParams);
982 
983  RCP<Factory> rFactory = rcp(new ReitzingerPFactory());
984  rFactory->SetParameterList(rParams);
985 
986  // These are all going to be user provided, so NoFactory
987  rFactory->SetFactory("Pnodal", NoFactory::getRCP());
988  rFactory->SetFactory("NodeAggMatrix", NoFactory::getRCP());
989  // rFactory->SetFactory("NodeMatrix", NoFactory::getRCP());
990 
991  if (levelID > 1)
992  rFactory->SetFactory("D0", this->GetFactoryManager(levelID - 1)->GetFactory("D0"));
993  else
994  rFactory->SetFactory("D0", NoFactory::getRCP());
995 
996  manager.SetFactory("Ptent", rFactory);
997  manager.SetFactory("D0", rFactory);
998  manager.SetFactory("InPlaceMap", rFactory);
999 }
1000 
1001 // =====================================================================================================
1002 // ========================================= TentativeP=================================================
1003 // =====================================================================================================
1004 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1007  FactoryManager& manager, int levelID, std::vector<keep_pair>& keeps) const {
1008  using strings = std::unordered_set<std::string>;
1009 
1010  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
1011 
1012  MUELU_SET_VAR_2LIST(paramList, defaultList, "aggregation: type", std::string, aggType);
1013  TEUCHOS_TEST_FOR_EXCEPTION(!strings({"uncoupled", "coupled", "brick", "matlab", "notay", "classical"}).count(aggType),
1014  Exceptions::RuntimeError, "Unknown aggregation algorithm: \"" << aggType << "\". Please consult User's Guide.");
1015 
1016  // Only doing this for classical because otherwise, the gold tests get broken badly
1017  RCP<AmalgamationFactory> amalgFact;
1018  if (aggType == "classical") {
1019  amalgFact = rcp(new AmalgamationFactory());
1020  manager.SetFactory("UnAmalgamationInfo", amalgFact);
1021  }
1022 
1023  // Aggregation graph
1024  RCP<Factory> dropFactory;
1025 
1026  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "matlab")) {
1027 #ifdef HAVE_MUELU_MATLAB
1028  dropFactory = rcp(new SingleLevelMatlabFactory());
1029  ParameterList socParams = paramList.sublist("strength-of-connection: params");
1030  dropFactory->SetParameterList(socParams);
1031 #else
1032  throw std::runtime_error("Cannot use MATLAB evolutionary strength-of-connection - MueLu was not configured with MATLAB support.");
1033 #endif
1034  } else if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "unsupported vector smoothing")) {
1036  ParameterList dropParams;
1037  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop scheme", std::string, dropParams);
1038  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: block diagonal: interleaved blocksize", int, dropParams);
1039  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: number of random vectors", int, dropParams);
1040  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: number of times to pre or post smooth", int, dropParams);
1041  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: penalty parameters", Teuchos::Array<double>, dropParams);
1042  dropFactory->SetParameterList(dropParams);
1043  } else {
1045  ParameterList dropParams;
1046  if (!rcp_dynamic_cast<CoalesceDropFactory>(dropFactory).is_null())
1047  dropParams.set("lightweight wrap", true);
1048  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop scheme", std::string, dropParams);
1049  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: row sum drop tol", double, dropParams);
1050  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: block diagonal: interleaved blocksize", int, dropParams);
1051  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol", double, dropParams);
1052  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: use ml scaling of drop tol", bool, dropParams);
1053 
1054  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: Dirichlet threshold", double, dropParams);
1055  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: greedy Dirichlet", bool, dropParams);
1056  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: distance laplacian metric", std::string, dropParams);
1057 #ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
1058  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: distance laplacian algo", std::string, dropParams);
1059  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: classical algo", std::string, dropParams);
1060 #endif
1061  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: distance laplacian directional weights", Teuchos::Array<double>, dropParams);
1062  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: coloring: localize color graph", bool, dropParams);
1063  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: dropping may create Dirichlet", bool, dropParams);
1064  if (useKokkos_) {
1065  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: use blocking", bool, dropParams);
1066  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: symmetrize graph after dropping", bool, dropParams);
1067  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: strength-of-connection: matrix", std::string, dropParams);
1068  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: strength-of-connection: measure", std::string, dropParams);
1069  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use lumping", bool, dropParams);
1070  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse graph", bool, dropParams);
1071  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse eigenvalue", bool, dropParams);
1072  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use root stencil", bool, dropParams);
1073  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: Dirichlet threshold", double, dropParams);
1074  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use spread lumping", bool, dropParams);
1075  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom growth factor", double, dropParams);
1076  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom cap", double, dropParams);
1077  }
1078 
1079 #ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
1080  if (!dropParams.isParameter("aggregation: drop scheme") ||
1081  (dropParams.isParameter("aggregation: drop scheme") &&
1082  ((dropParams.get<std::string>("aggregation: drop scheme") != "point-wise") && (dropParams.get<std::string>("aggregation: drop scheme") != "cut-drop")))) {
1083  Teuchos::ParameterList dropParamsWithDefaults(dropParams);
1084 
1085 #define MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(paramList, paramName, paramType) \
1086  if (!paramList.isParameter(paramName)) { \
1087  paramList.set(paramName, MasterList::getDefault<paramType>(paramName)); \
1088  }
1089 
1090  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: drop scheme", std::string);
1091  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: strength-of-connection: matrix", std::string);
1092  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: strength-of-connection: measure", std::string);
1093  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: use blocking", bool);
1094 
1095 #undef MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST
1096 
1097  // We are using the old style of dropping params
1098  TEUCHOS_TEST_FOR_EXCEPTION(dropParams.isParameter("aggregation: strength-of-connection: matrix") ||
1099  dropParams.isParameter("aggregation: strength-of-connection: measure") ||
1100  dropParams.isParameter("aggregation: use blocking"),
1102  "The inputs contain a mix of old and new dropping parameters:\n\n"
1103  << dropParams << "\n\nKeep in mind that defaults are set for old parameters, so this gets interpreted as\n\n"
1104  << dropParamsWithDefaults);
1105  }
1106 #endif
1107 
1108  if (!amalgFact.is_null())
1109  dropFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
1110 
1111  if (dropParams.isParameter("aggregation: drop scheme")) {
1112  std::string drop_scheme = dropParams.get<std::string>("aggregation: drop scheme");
1113  if (drop_scheme == "block diagonal colored signed classical")
1114  manager.SetFactory("Coloring Graph", dropFactory);
1115  if ((MUELU_TEST_PARAM_2LIST(dropParams, defaultList, "aggregation: use blocking", bool, true)) ||
1116  (drop_scheme.find("block diagonal") != std::string::npos || drop_scheme == "signed classical")) {
1117  if (levelID > 0)
1118  dropFactory->SetFactory("BlockNumber", this->GetFactoryManager(levelID - 1)->GetFactory("BlockNumber"));
1119  else
1120  dropFactory->SetFactory("BlockNumber", manager.GetFactory("BlockNumber"));
1121  }
1122  }
1123 
1124  dropFactory->SetParameterList(dropParams);
1125  }
1126  manager.SetFactory("Graph", dropFactory);
1127 
1128 // Aggregation scheme
1129 #ifndef HAVE_MUELU_MATLAB
1130  if (aggType == "matlab")
1131  throw std::runtime_error("Cannot use MATLAB aggregation - MueLu was not configured with MATLAB support.");
1132 #endif
1133  RCP<Factory> aggFactory;
1134  if (aggType == "uncoupled") {
1135  aggFactory = rcp(new UncoupledAggregationFactory());
1136  ParameterList aggParams;
1137  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: mode", std::string, aggParams);
1138  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: ordering", std::string, aggParams);
1139  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: min agg size", int, aggParams);
1140  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: max agg size", int, aggParams);
1141  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: max selected neighbors", int, aggParams);
1142  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: backend", std::string, aggParams);
1143  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: phase 1 algorithm", std::string, aggParams);
1144  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: deterministic", bool, aggParams);
1145  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: coloring algorithm", std::string, aggParams);
1146  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 1", bool, aggParams);
1147  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 2a", bool, aggParams);
1148  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 2b", bool, aggParams);
1149  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 3", bool, aggParams);
1150  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: match ML phase1", bool, aggParams);
1151  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: match ML phase2a", bool, aggParams);
1152  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: match ML phase2b", bool, aggParams);
1153  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: phase2a agg factor", double, aggParams);
1154  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: preserve Dirichlet points", bool, aggParams);
1155  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: error on nodes with no on-rank neighbors", bool, aggParams);
1156  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: phase3 avoid singletons", bool, aggParams);
1157  aggFactory->SetParameterList(aggParams);
1158  // make sure that the aggregation factory has all necessary data
1159  aggFactory->SetFactory("DofsPerNode", manager.GetFactory("Graph"));
1160  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
1161  // aggFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
1162 
1163  } else if (aggType == "brick") {
1164  aggFactory = rcp(new BrickAggregationFactory());
1165  ParameterList aggParams;
1166  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick x size", int, aggParams);
1167  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick y size", int, aggParams);
1168  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick z size", int, aggParams);
1169  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick x Dirichlet", bool, aggParams);
1170  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick y Dirichlet", bool, aggParams);
1171  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick z Dirichlet", bool, aggParams);
1172  aggFactory->SetParameterList(aggParams);
1173 
1174  // Unlike other factories, BrickAggregationFactory makes the Graph/DofsPerNode itself
1175  manager.SetFactory("Graph", aggFactory);
1176  manager.SetFactory("DofsPerNode", aggFactory);
1177  manager.SetFactory("Filtering", aggFactory);
1178  if (levelID > 1) {
1179  // We check for levelID > 0, as in the interpreter aggFactory for
1180  // levelID really corresponds to level 0. Managers are clunky, as they
1181  // contain factories for two different levels
1182  aggFactory->SetFactory("Coordinates", this->GetFactoryManager(levelID - 1)->GetFactory("Coordinates"));
1183  }
1184  } else if (aggType == "classical") {
1185  // Map and coloring
1186  RCP<Factory> mapFact = rcp(new ClassicalMapFactory());
1187  ParameterList mapParams;
1188  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: deterministic", bool, mapParams);
1189  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: coloring algorithm", std::string, mapParams);
1190 
1191  ParameterList tempParams;
1192  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop scheme", std::string, tempParams);
1193  std::string drop_algo = tempParams.get<std::string>("aggregation: drop scheme");
1194  if (drop_algo == "block diagonal colored signed classical") {
1195  mapParams.set("aggregation: coloring: use color graph", true);
1196  mapFact->SetFactory("Coloring Graph", manager.GetFactory("Coloring Graph"));
1197  }
1198  mapFact->SetParameterList(mapParams);
1199  mapFact->SetFactory("Graph", manager.GetFactory("Graph"));
1200  mapFact->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
1201 
1202  manager.SetFactory("FC Splitting", mapFact);
1203  manager.SetFactory("CoarseMap", mapFact);
1204 
1205  aggFactory = rcp(new ClassicalPFactory());
1206  ParameterList aggParams;
1207  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: classical scheme", std::string, aggParams);
1208  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop scheme", std::string, aggParams);
1209  aggFactory->SetParameterList(aggParams);
1210  aggFactory->SetFactory("FC Splitting", manager.GetFactory("FC Splitting"));
1211  aggFactory->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1212  aggFactory->SetFactory("DofsPerNode", manager.GetFactory("Graph"));
1213  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
1214 
1215  if (drop_algo.find("block diagonal") != std::string::npos || drop_algo == "signed classical") {
1216  if (levelID > 0)
1217  aggFactory->SetFactory("BlockNumber", this->GetFactoryManager(levelID - 1)->GetFactory("BlockNumber"));
1218  else
1219  aggFactory->SetFactory("BlockNumber", manager.GetFactory("BlockNumber"));
1220  }
1221 
1222  // Now we short-circuit, because we neither need nor want TentativePFactory here
1223  manager.SetFactory("Ptent", aggFactory);
1224  manager.SetFactory("P Graph", aggFactory);
1225 
1226  if (reuseType == "tP" && levelID) {
1227  // keeps.push_back(keep_pair("Nullspace", Ptent.get()));
1228  keeps.push_back(keep_pair("Ptent", aggFactory.get()));
1229  }
1230  return;
1231  } else if (aggType == "notay") {
1232  aggFactory = rcp(new NotayAggregationFactory());
1233  ParameterList aggParams;
1234  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: pairwise: size", int, aggParams);
1235  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: pairwise: tie threshold", double, aggParams);
1236  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: Dirichlet threshold", double, aggParams);
1237  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: ordering", std::string, aggParams);
1238  aggFactory->SetParameterList(aggParams);
1239  aggFactory->SetFactory("DofsPerNode", manager.GetFactory("Graph"));
1240  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
1241  }
1242 #ifdef HAVE_MUELU_MATLAB
1243  else if (aggType == "matlab") {
1244  ParameterList aggParams = paramList.sublist("aggregation: params");
1245  aggFactory = rcp(new SingleLevelMatlabFactory());
1246  aggFactory->SetParameterList(aggParams);
1247  }
1248 #endif
1249 
1250  manager.SetFactory("Aggregates", aggFactory);
1251 
1252  // Coarse map
1253  RCP<Factory> coarseMap = rcp(new CoarseMapFactory());
1254  coarseMap->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1255  manager.SetFactory("CoarseMap", coarseMap);
1256 
1257  // Tentative P
1259  ParameterList ptentParams;
1260  if (paramList.isSublist("matrixmatrix: kernel params"))
1261  ptentParams.sublist("matrixmatrix: kernel params", false) = paramList.sublist("matrixmatrix: kernel params");
1262  if (defaultList.isSublist("matrixmatrix: kernel params"))
1263  ptentParams.sublist("matrixmatrix: kernel params", false) = defaultList.sublist("matrixmatrix: kernel params");
1264  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: calculate qr", bool, ptentParams);
1265  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: build coarse coordinates", bool, ptentParams);
1266  Ptent->SetParameterList(ptentParams);
1267  Ptent->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1268  Ptent->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1269  manager.SetFactory("Ptent", Ptent);
1270 
1271  if (reuseType == "tP" && levelID) {
1272  keeps.push_back(keep_pair("Nullspace", Ptent.get()));
1273  keeps.push_back(keep_pair("P", Ptent.get()));
1274  }
1275 }
1276 
1277 // =====================================================================================================
1278 // ============================================ RAP ====================================================
1279 // =====================================================================================================
1280 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1282  UpdateFactoryManager_RAP(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager,
1283  int levelID, std::vector<keep_pair>& keeps) const {
1284  if (paramList.isParameter("A") && !paramList.get<RCP<Matrix> >("A").is_null()) {
1285  // We have user matrix A
1286  manager.SetFactory("A", NoFactory::getRCP());
1287  return;
1288  }
1289 
1290  ParameterList RAPparams;
1291 
1292  RCP<RAPFactory> RAP;
1293  RCP<RAPShiftFactory> RAPs;
1294  // Allow for Galerkin or shifted RAP
1295  // FIXME: Should this not be some form of MUELU_SET_VAR_2LIST?
1296  std::string alg = paramList.get("rap: algorithm", "galerkin");
1297  if (alg == "shift" || alg == "non-galerkin") {
1298  RAPs = rcp(new RAPShiftFactory());
1299  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift", double, RAPparams);
1300  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift diagonal M", bool, RAPparams);
1301  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift low storage", bool, RAPparams);
1302  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift array", Teuchos::Array<double>, RAPparams);
1303  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: cfl array", Teuchos::Array<double>, RAPparams);
1304 
1305  } else {
1306  RAP = rcp(new RAPFactory());
1307  }
1308 
1309  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: relative diagonal floor", Teuchos::Array<double>, RAPparams);
1310 
1311  if (paramList.isSublist("matrixmatrix: kernel params"))
1312  RAPparams.sublist("matrixmatrix: kernel params", false) = paramList.sublist("matrixmatrix: kernel params");
1313  if (defaultList.isSublist("matrixmatrix: kernel params"))
1314  RAPparams.sublist("matrixmatrix: kernel params", false) = defaultList.sublist("matrixmatrix: kernel params");
1315  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "transpose: use implicit", bool, RAPparams);
1316  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: fix zero diagonals", bool, RAPparams);
1317  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: fix zero diagonals threshold", double, RAPparams);
1318  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: fix zero diagonals replacement", Scalar, RAPparams);
1319 
1320  // if "rap: triple product" has not been set and algorithm is "unsmoothed" switch triple product on
1321  if (!paramList.isParameter("rap: triple product") &&
1322  paramList.isType<std::string>("multigrid algorithm") &&
1323  paramList.get<std::string>("multigrid algorithm") == "unsmoothed")
1324  paramList.set("rap: triple product", true);
1325  else
1326  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: triple product", bool, RAPparams);
1327 
1328  try {
1329  if (paramList.isParameter("aggregation: allow empty prolongator columns")) {
1330  RAPparams.set("CheckMainDiagonal", paramList.get<bool>("aggregation: allow empty prolongator columns"));
1331  RAPparams.set("RepairMainDiagonal", paramList.get<bool>("aggregation: allow empty prolongator columns"));
1332  } else if (defaultList.isParameter("aggregation: allow empty prolongator columns")) {
1333  RAPparams.set("CheckMainDiagonal", defaultList.get<bool>("aggregation: allow empty prolongator columns"));
1334  RAPparams.set("RepairMainDiagonal", defaultList.get<bool>("aggregation: allow empty prolongator columns"));
1335  }
1336 
1339  "Error: parameter \"aggregation: allow empty prolongator columns\" must be of type " << Teuchos::TypeNameTraits<bool>::name());
1340  }
1341 
1342  if (!RAP.is_null()) {
1343  RAP->SetParameterList(RAPparams);
1344  RAP->SetFactory("P", manager.GetFactory("P"));
1345  } else {
1346  RAPs->SetParameterList(RAPparams);
1347  RAPs->SetFactory("P", manager.GetFactory("P"));
1348  }
1349 
1350  if (!this->implicitTranspose_) {
1351  if (!RAP.is_null())
1352  RAP->SetFactory("R", manager.GetFactory("R"));
1353  else
1354  RAPs->SetFactory("R", manager.GetFactory("R"));
1355  }
1356 
1357  // Matrix analysis
1358  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "matrix: compute analysis", bool, true)) {
1359  RCP<Factory> matrixAnalysisFact = rcp(new MatrixAnalysisFactory());
1360 
1361  if (!RAP.is_null())
1362  RAP->AddTransferFactory(matrixAnalysisFact);
1363  else
1364  RAPs->AddTransferFactory(matrixAnalysisFact);
1365  }
1366 
1367  // Aggregate qualities
1368  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: compute aggregate qualities", bool, true)) {
1369  RCP<Factory> aggQualityFact = rcp(new AggregateQualityEstimateFactory());
1370  ParameterList aggQualityParams;
1371  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: good aggregate threshold", double, aggQualityParams);
1372  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: file output", bool, aggQualityParams);
1373  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: file base", std::string, aggQualityParams);
1374  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: check symmetry", bool, aggQualityParams);
1375  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: algorithm", std::string, aggQualityParams);
1376  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: zero threshold", double, aggQualityParams);
1377  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: percentiles", Teuchos::Array<double>, aggQualityParams);
1378  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: mode", std::string, aggQualityParams);
1379  aggQualityFact->SetParameterList(aggQualityParams);
1380  aggQualityFact->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1381  aggQualityFact->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1382  manager.SetFactory("AggregateQualities", aggQualityFact);
1383 
1384  if (!RAP.is_null())
1385  RAP->AddTransferFactory(aggQualityFact);
1386  else
1387  RAPs->AddTransferFactory(aggQualityFact);
1388  }
1389 
1390  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: export visualization data", bool, true)) {
1392  ParameterList aggExportParams;
1393  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output filename", std::string, aggExportParams);
1394  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: agg style", std::string, aggExportParams);
1395  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: iter", int, aggExportParams);
1396  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: time step", int, aggExportParams);
1397  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: fine graph edges", bool, aggExportParams);
1398  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: coarse graph edges", bool, aggExportParams);
1399  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: build colormap", bool, aggExportParams);
1400  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: aggregate qualities", bool, aggExportParams);
1401  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: material", bool, aggExportParams);
1402  aggExport->SetParameterList(aggExportParams);
1403  aggExport->SetFactory("AggregateQualities", manager.GetFactory("AggregateQualities"));
1404  aggExport->SetFactory("DofsPerNode", manager.GetFactory("DofsPerNode"));
1405  aggExport->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1406  aggExport->SetFactory("Graph", manager.GetFactory("Graph"));
1407 
1408  if (!RAP.is_null())
1409  RAP->AddTransferFactory(aggExport);
1410  else
1411  RAPs->AddTransferFactory(aggExport);
1412  }
1413  if (!RAP.is_null())
1414  manager.SetFactory("A", RAP);
1415  else
1416  manager.SetFactory("A", RAPs);
1417 
1418  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
1419  MUELU_SET_VAR_2LIST(paramList, defaultList, "sa: use filtered matrix", bool, useFiltering);
1420  bool filteringChangesMatrix = useFiltering && !MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol", double, 0);
1421 
1422  if (reuseType == "RP" || (reuseType == "tP" && !filteringChangesMatrix)) {
1423  if (!RAP.is_null()) {
1424  keeps.push_back(keep_pair("AP reuse data", RAP.get()));
1425  keeps.push_back(keep_pair("RAP reuse data", RAP.get()));
1426 
1427  } else {
1428  keeps.push_back(keep_pair("AP reuse data", RAPs.get()));
1429  keeps.push_back(keep_pair("RAP reuse data", RAPs.get()));
1430  }
1431  }
1432 }
1433 
1434 // =====================================================================================================
1435 // ======================================= Coordinates =================================================
1436 // =====================================================================================================
1437 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1440  FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
1441  bool have_userCO = false;
1442  if (paramList.isParameter("Coordinates") && !paramList.get<RCP<MultiVector> >("Coordinates").is_null())
1443  have_userCO = true;
1444 
1445  if (useCoordinates_) {
1446  if (have_userCO) {
1447  manager.SetFactory("Coordinates", NoFactory::getRCP());
1448 
1449  } else {
1451  coords->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1452  coords->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1453  manager.SetFactory("Coordinates", coords);
1454 
1455  auto RAP = rcp_const_cast<RAPFactory>(rcp_dynamic_cast<const RAPFactory>(manager.GetFactory("A")));
1456  if (!RAP.is_null()) {
1457  RAP->AddTransferFactory(manager.GetFactory("Coordinates"));
1458  } else {
1459  auto RAPs = rcp_const_cast<RAPShiftFactory>(rcp_dynamic_cast<const RAPShiftFactory>(manager.GetFactory("A")));
1460  RAPs->AddTransferFactory(manager.GetFactory("Coordinates"));
1461  }
1462  }
1463  }
1464 }
1465 
1466 // ======================================================================================================
1467 // ======================================== Material ==================================================
1468 // =====================================================================================================
1469 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1471  UpdateFactoryManager_Material(ParameterList& paramList, const ParameterList& /* defaultList */,
1472  FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
1473  bool have_userMaterial = false;
1474  if (paramList.isParameter("Material") && !paramList.get<RCP<MultiVector> >("Material").is_null())
1475  have_userMaterial = true;
1476 
1477  if (useMaterial_) {
1478  if (have_userMaterial) {
1479  manager.SetFactory("Material", NoFactory::getRCP());
1480  } else {
1481  RCP<Factory> materialTransfer = rcp(new MultiVectorTransferFactory());
1482  ParameterList materialTransferParameters;
1483  materialTransferParameters.set("Vector name", "Material");
1484  materialTransferParameters.set("Transfer name", "Aggregates");
1485  materialTransferParameters.set("Normalize", true);
1486  materialTransfer->SetParameterList(materialTransferParameters);
1487  materialTransfer->SetFactory("Transfer factory", manager.GetFactory("Aggregates"));
1488  materialTransfer->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1489  manager.SetFactory("Material", materialTransfer);
1490 
1491  auto RAP = rcp_const_cast<RAPFactory>(rcp_dynamic_cast<const RAPFactory>(manager.GetFactory("A")));
1492  if (!RAP.is_null()) {
1493  RAP->AddTransferFactory(manager.GetFactory("Material"));
1494  } else {
1495  auto RAPs = rcp_const_cast<RAPShiftFactory>(rcp_dynamic_cast<const RAPShiftFactory>(manager.GetFactory("A")));
1496  RAPs->AddTransferFactory(manager.GetFactory("Material"));
1497  }
1498  }
1499  }
1500 }
1501 
1502 // =====================================================================================================
1503 // ================================= LocalOrdinalTransfer =============================================
1504 // =====================================================================================================
1505 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1507  UpdateFactoryManager_LocalOrdinalTransfer(const std::string& VarName, const std::string& multigridAlgo, ParameterList& paramList, const ParameterList& /* defaultList */,
1508  FactoryManager& manager, int levelID, std::vector<keep_pair>& /* keeps */) const {
1509  // NOTE: You would think this would be levelID > 0, but you'd be wrong, since the FactoryManager is basically
1510  // offset by a level from the things which actually do the work.
1511  if (useBlockNumber_ && (levelID > 0)) {
1512  auto RAP = rcp_const_cast<RAPFactory>(rcp_dynamic_cast<const RAPFactory>(manager.GetFactory("A")));
1513  auto RAPs = rcp_const_cast<RAPShiftFactory>(rcp_dynamic_cast<const RAPShiftFactory>(manager.GetFactory("A")));
1514  if (!RAP.is_null() || !RAPs.is_null()) {
1515  RCP<Factory> fact = rcp(new LocalOrdinalTransferFactory(VarName, multigridAlgo));
1516  if (multigridAlgo == "classical")
1517  fact->SetFactory("P Graph", manager.GetFactory("P Graph"));
1518  else
1519  fact->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1520  fact->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1521 
1522  fact->SetFactory(VarName, this->GetFactoryManager(levelID - 1)->GetFactory(VarName));
1523 
1524  manager.SetFactory(VarName, fact);
1525 
1526  if (!RAP.is_null())
1527  RAP->AddTransferFactory(manager.GetFactory(VarName));
1528  else
1529  RAPs->AddTransferFactory(manager.GetFactory(VarName));
1530  }
1531  }
1532 }
1533 
1534 // ======================================================================================================
1535 // ====================================== BlockNumber =================================================
1536 // =====================================================================================================
1537 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1540  FactoryManager& manager, int levelID, std::vector<keep_pair>& keeps) const {
1541  if (useBlockNumber_) {
1542  ParameterList myParams;
1544  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: block diagonal: interleaved blocksize", int, myParams);
1545  fact->SetParameterList(myParams);
1546  manager.SetFactory("BlockNumber", fact);
1547  }
1548 }
1549 
1550 // =====================================================================================================
1551 // =========================================== Restriction =============================================
1552 // =====================================================================================================
1553 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1556  int levelID, std::vector<keep_pair>& /* keeps */) const {
1557  MUELU_SET_VAR_2LIST(paramList, defaultList, "multigrid algorithm", std::string, multigridAlgo);
1558  bool have_userR = false;
1559  if (paramList.isParameter("R") && !paramList.get<RCP<Matrix> >("R").is_null())
1560  have_userR = true;
1561 
1562  // === Restriction ===
1563  RCP<Factory> R;
1564  if (!this->implicitTranspose_) {
1565  MUELU_SET_VAR_2LIST(paramList, defaultList, "problem: symmetric", bool, isSymmetric);
1566 
1567  if (isSymmetric == false && (multigridAlgo == "unsmoothed" || multigridAlgo == "emin")) {
1568  this->GetOStream(Warnings0) << "Switching \"problem: symmetric\" parameter to symmetric as multigrid algorithm. " << multigridAlgo << " is primarily supposed to be used for symmetric problems.\n\n"
1569  << "Please note: if you are using \"unsmoothed\" transfer operators the \"problem: symmetric\" parameter "
1570  << "has no real mathematical meaning, i.e. you can use it for non-symmetric\n"
1571  << "problems, too. With \"problem: symmetric\"=\"symmetric\" you can use implicit transpose for building "
1572  << "the restriction operators which may drastically reduce the amount of consumed memory." << std::endl;
1573  isSymmetric = true;
1574  }
1575  TEUCHOS_TEST_FOR_EXCEPTION(multigridAlgo == "pg" && isSymmetric == true, Exceptions::RuntimeError,
1576  "Petrov-Galerkin smoothed transfer operators are only allowed for non-symmetric problems: Set \"problem: symmetric\" to false!\n"
1577  "While PG smoothed transfer operators generally would also work for symmetric problems this is an unusual use case. "
1578  "You can use the factory-based xml interface though if you need PG-AMG for symmetric problems.");
1579 
1580  if (have_userR) {
1581  manager.SetFactory("R", NoFactory::getRCP());
1582  } else {
1583  if (isSymmetric)
1584  R = rcp(new TransPFactory());
1585  else
1586  R = rcp(new GenericRFactory());
1587 
1588  R->SetFactory("P", manager.GetFactory("P"));
1589  manager.SetFactory("R", R);
1590  }
1591 
1592  } else {
1593  manager.SetFactory("R", Teuchos::null);
1594  }
1595 
1596  // === Restriction: Nullspace Scaling ===
1597  if (paramList.isParameter("restriction: scale nullspace") && paramList.get<bool>("restriction: scale nullspace")) {
1598  RCP<TentativePFactory> tentPFactory = rcp(new TentativePFactory());
1599  Teuchos::ParameterList tentPlist;
1600  tentPlist.set("Nullspace name", "Scaled Nullspace");
1601  tentPFactory->SetParameterList(tentPlist);
1602  tentPFactory->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1603  tentPFactory->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1604 
1605  if (R.is_null()) R = rcp(new TransPFactory());
1606  R->SetFactory("P", tentPFactory);
1607  }
1608 }
1609 
1610 // =====================================================================================================
1611 // ========================================= Repartition ===============================================
1612 // =====================================================================================================
1613 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1616  int levelID, std::vector<keep_pair>& keeps, RCP<Factory>& nullSpaceFactory) const {
1617  // === Repartitioning ===
1618  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
1619  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: enable", bool, enableRepart);
1620  if (enableRepart) {
1621 #if defined(HAVE_MPI) && (defined(HAVE_MUELU_ZOLTAN) || defined(HAVE_MUELU_ZOLTAN2)) // skip to the end, print warning, and turn off repartitioning if we don't have MPI and Zoltan/Zoltan2
1622  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: use subcommunicators in place", bool, enableInPlace);
1623  // Short summary of the issue: RebalanceTransferFactory shares ownership
1624  // of "P" with SaPFactory, and therefore, changes the stored version.
1625  // That means that if SaPFactory generated P, and stored it on the level,
1626  // then after rebalancing the value in that storage changed. It goes
1627  // against the concept of factories (I think), that every factory is
1628  // responsible for its own objects, and they are immutable outside.
1629  //
1630  // In reuse, this is what happens: as we reuse Importer across setups,
1631  // the order of factories changes, and coupled with shared ownership
1632  // leads to problems.
1633  // *First setup*
1634  // SaP builds P [and stores it]
1635  // TransP builds R [and stores it]
1636  // RAP builds A [and stores it]
1637  // RebalanceTransfer rebalances P [and changes the P stored by SaP] (*)
1638  // RebalanceTransfer rebalances R
1639  // RebalanceAc rebalances A
1640  // *Second setup* ("RP" reuse)
1641  // RebalanceTransfer rebalances P [which is incorrect due to (*)]
1642  // RebalanceTransfer rebalances R
1643  // RAP builds A [which is incorrect due to (*)]
1644  // RebalanceAc rebalances A [which throws due to map inconsistency]
1645  // ...
1646  // *Second setup* ("tP" reuse)
1647  // SaP builds P [and stores it]
1648  // RebalanceTransfer rebalances P [and changes the P stored by SaP] (**)
1649  // TransP builds R [which is incorrect due to (**)]
1650  // RebalanceTransfer rebalances R
1651  // ...
1652  //
1653  // Couple solutions to this:
1654  // 1. [implemented] Requre "tP" and "PR" reuse to only be used with
1655  // implicit rebalancing.
1656  // 2. Do deep copy of P, and changed domain map and importer there.
1657  // Need to investigate how expensive this is.
1658  TEUCHOS_TEST_FOR_EXCEPTION(this->doPRrebalance_ && (reuseType == "tP" || reuseType == "RP"), Exceptions::InvalidArgument,
1659  "Reuse types \"tP\" and \"PR\" require \"repartition: rebalance P and R\" set to \"false\"");
1660 
1661  // TEUCHOS_TEST_FOR_EXCEPTION(aggType == "brick", Exceptions::InvalidArgument,
1662  // "Aggregation type \"brick\" requires \"repartition: enable\" set to \"false\"");
1663 
1664  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: partitioner", std::string, partName);
1665  TEUCHOS_TEST_FOR_EXCEPTION(partName != "zoltan" && partName != "zoltan2", Exceptions::InvalidArgument,
1666  "Invalid partitioner name: \"" << partName << "\". Valid options: \"zoltan\", \"zoltan2\"");
1667 
1668 #ifndef HAVE_MUELU_ZOLTAN
1669  bool switched = false;
1670  if (partName == "zoltan") {
1671  this->GetOStream(Warnings0) << "Zoltan interface is not available, trying to switch to Zoltan2" << std::endl;
1672  partName = "zoltan2";
1673  switched = true;
1674  }
1675 #else
1676 #ifndef HAVE_MUELU_ZOLTAN2
1677  bool switched = false;
1678 #endif // HAVE_MUELU_ZOLTAN2
1679 #endif // HAVE_MUELU_ZOLTAN
1680 
1681 #ifndef HAVE_MUELU_ZOLTAN2
1682  if (partName == "zoltan2" && !switched) {
1683  this->GetOStream(Warnings0) << "Zoltan2 interface is not available, trying to switch to Zoltan" << std::endl;
1684  partName = "zoltan";
1685  }
1686 #endif // HAVE_MUELU_ZOLTAN2
1687 
1688  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: node repartition level", int, nodeRepartitionLevel);
1689 
1690  // RepartitionHeuristic
1691  auto repartheurFactory = rcp(new RepartitionHeuristicFactory());
1692  ParameterList repartheurParams;
1693  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: node repartition level", int, repartheurParams);
1694  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: start level", int, repartheurParams);
1695  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: min rows per proc", int, repartheurParams);
1696  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: target rows per proc", int, repartheurParams);
1697  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: min rows per thread", int, repartheurParams);
1698  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: target rows per thread", int, repartheurParams);
1699  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: max imbalance", double, repartheurParams);
1700  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: put on single proc", int, repartheurParams);
1701  repartheurFactory->SetParameterList(repartheurParams);
1702  repartheurFactory->SetFactory("A", manager.GetFactory("A"));
1703  manager.SetFactory("number of partitions", repartheurFactory);
1704  manager.SetFactory("repartition: heuristic target rows per process", repartheurFactory);
1705 
1706  // Partitioner
1707  RCP<Factory> partitioner;
1708  if (levelID == nodeRepartitionLevel) {
1709  // partitioner = rcp(new NodePartitionInterface());
1711  ParameterList partParams;
1712  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: node id", int, repartheurParams);
1713  partitioner->SetParameterList(partParams);
1714  partitioner->SetFactory("Node Comm", manager.GetFactory("Node Comm"));
1715  } else if (partName == "zoltan") {
1716 #ifdef HAVE_MUELU_ZOLTAN
1717  partitioner = rcp(new ZoltanInterface());
1718  // NOTE: ZoltanInterface ("zoltan") does not support external parameters through ParameterList
1719 #else
1720  throw Exceptions::RuntimeError("Zoltan interface is not available");
1721 #endif // HAVE_MUELU_ZOLTAN
1722  } else if (partName == "zoltan2") {
1723 #ifdef HAVE_MUELU_ZOLTAN2
1724  partitioner = rcp(new Zoltan2Interface());
1725  ParameterList partParams;
1726  RCP<const ParameterList> partpartParams = rcp(new ParameterList(paramList.sublist("repartition: params", false)));
1727  partParams.set("ParameterList", partpartParams);
1728  partitioner->SetParameterList(partParams);
1729  partitioner->SetFactory("repartition: heuristic target rows per process",
1730  manager.GetFactory("repartition: heuristic target rows per process"));
1731 #else
1732  throw Exceptions::RuntimeError("Zoltan2 interface is not available");
1733 #endif // HAVE_MUELU_ZOLTAN2
1734  }
1735 
1736  partitioner->SetFactory("A", manager.GetFactory("A"));
1737  partitioner->SetFactory("number of partitions", manager.GetFactory("number of partitions"));
1738  if (useCoordinates_)
1739  partitioner->SetFactory("Coordinates", manager.GetFactory("Coordinates"));
1740  manager.SetFactory("Partition", partitioner);
1741 
1742  // Repartitioner
1743  auto repartFactory = rcp(new RepartitionFactory());
1744  ParameterList repartParams;
1745  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: print partition distribution", bool, repartParams);
1746  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap parts", bool, repartParams);
1747  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap num values", int, repartParams);
1748  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: save importer", bool, repartParams);
1749  repartFactory->SetParameterList(repartParams);
1750  repartFactory->SetFactory("A", manager.GetFactory("A"));
1751  repartFactory->SetFactory("number of partitions", manager.GetFactory("number of partitions"));
1752  repartFactory->SetFactory("Partition", manager.GetFactory("Partition"));
1753  manager.SetFactory("Importer", repartFactory);
1754  if (reuseType != "none" && reuseType != "S" && levelID)
1755  keeps.push_back(keep_pair("Importer", manager.GetFactory("Importer").get()));
1756 
1757  if (enableInPlace) {
1758  // Rebalanced A (in place)
1759  // NOTE: This is for when we want to constrain repartitioning to match some other idea of what's going on.
1760  // The major application is the (1,1) hierarchy in the Maxwell1 preconditioner.
1761  auto newA = rcp(new RebalanceAcFactory());
1762  ParameterList rebAcParams;
1763  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, rebAcParams);
1764  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators in place", bool, rebAcParams);
1765  newA->SetParameterList(rebAcParams);
1766  newA->SetFactory("A", manager.GetFactory("A"));
1767  newA->SetFactory("InPlaceMap", manager.GetFactory("InPlaceMap"));
1768  manager.SetFactory("A", newA);
1769  } else {
1770  // Rebalanced A
1771  auto newA = rcp(new RebalanceAcFactory());
1772  ParameterList rebAcParams;
1773  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, rebAcParams);
1774  newA->SetParameterList(rebAcParams);
1775  newA->SetFactory("A", manager.GetFactory("A"));
1776  newA->SetFactory("Importer", manager.GetFactory("Importer"));
1777  manager.SetFactory("A", newA);
1778 
1779  // Rebalanced P
1780  auto newP = rcp(new RebalanceTransferFactory());
1781  ParameterList newPparams;
1782  newPparams.set("type", "Interpolation");
1783  if (changedPRrebalance_)
1784  newPparams.set("repartition: rebalance P and R", this->doPRrebalance_);
1785  if (changedPRViaCopyrebalance_)
1786  newPparams.set("repartition: explicit via new copy rebalance P and R", true);
1787  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, newPparams);
1788  newP->SetParameterList(newPparams);
1789  newP->SetFactory("Importer", manager.GetFactory("Importer"));
1790  newP->SetFactory("P", manager.GetFactory("P"));
1791  manager.SetFactory("P", newP);
1792  if (!paramList.isParameter("semicoarsen: number of levels"))
1793  newP->SetFactory("Nullspace", manager.GetFactory("Ptent"));
1794  else
1795  newP->SetFactory("Nullspace", manager.GetFactory("P")); // TogglePFactory
1796  if (useCoordinates_) {
1797  newP->SetFactory("Coordinates", manager.GetFactory("Coordinates"));
1798  manager.SetFactory("Coordinates", newP);
1799  }
1800  if (useMaterial_) {
1801  newP->SetFactory("Material", manager.GetFactory("Material"));
1802  manager.SetFactory("Material", newP);
1803  }
1804  if (useBlockNumber_ && (levelID > 0)) {
1805  newP->SetFactory("BlockNumber", manager.GetFactory("BlockNumber"));
1806  manager.SetFactory("BlockNumber", newP);
1807  }
1808 
1809  // Rebalanced R
1810  auto newR = rcp(new RebalanceTransferFactory());
1811  ParameterList newRparams;
1812  newRparams.set("type", "Restriction");
1813  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, newRparams);
1814  if (changedPRrebalance_)
1815  newRparams.set("repartition: rebalance P and R", this->doPRrebalance_);
1816  if (changedPRViaCopyrebalance_)
1817  newPparams.set("repartition: explicit via new copy rebalance P and R", true);
1818  if (changedImplicitTranspose_)
1819  newRparams.set("transpose: use implicit", this->implicitTranspose_);
1820  newR->SetParameterList(newRparams);
1821  newR->SetFactory("Importer", manager.GetFactory("Importer"));
1822  if (!this->implicitTranspose_) {
1823  newR->SetFactory("R", manager.GetFactory("R"));
1824  manager.SetFactory("R", newR);
1825  }
1826 
1827  // NOTE: the role of NullspaceFactory is to provide nullspace on the finest
1828  // level if a user does not do that. For all other levels it simply passes
1829  // nullspace from a real factory to whoever needs it. If we don't use
1830  // repartitioning, that factory is "TentativePFactory"; if we do, it is
1831  // "RebalanceTransferFactory". But we still have to have NullspaceFactory as
1832  // the "Nullspace" of the manager
1833  // NOTE: This really needs to be set on the *NullSpaceFactory*, not manager.get("Nullspace").
1834  ParameterList newNullparams;
1835  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "nullspace: calculate rotations", bool, newNullparams);
1836  nullSpaceFactory->SetFactory("Nullspace", newP);
1837  nullSpaceFactory->SetParameterList(newNullparams);
1838  }
1839 #else
1840  paramList.set("repartition: enable", false);
1841 #ifndef HAVE_MPI
1842  this->GetOStream(Warnings0) << "No repartitioning available for a serial run\n";
1843 #else
1844  this->GetOStream(Warnings0) << "Zoltan/Zoltan2 are unavailable for repartitioning\n";
1845 #endif // HAVE_MPI
1846 #endif // defined(HAVE_MPI) && (defined(HAVE_MUELU_ZOLTAN) || defined(HAVE_MUELU_ZOLTAN2))
1847  }
1848 }
1849 
1850 // =====================================================================================================
1851 // ========================================= Low precision transfers ===================================
1852 // =====================================================================================================
1853 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1856  int levelID, std::vector<keep_pair>& keeps) const {
1857  MUELU_SET_VAR_2LIST(paramList, defaultList, "transfers: half precision", bool, enableLowPrecision);
1858 
1859  if (enableLowPrecision) {
1860  // Low precision P
1861  auto newP = rcp(new LowPrecisionFactory());
1862  ParameterList newPparams;
1863  newPparams.set("matrix key", "P");
1864  newP->SetParameterList(newPparams);
1865  newP->SetFactory("P", manager.GetFactory("P"));
1866  manager.SetFactory("P", newP);
1867 
1868  if (!this->implicitTranspose_) {
1869  // Low precision R
1870  auto newR = rcp(new LowPrecisionFactory());
1871  ParameterList newRparams;
1872  newRparams.set("matrix key", "R");
1873  newR->SetParameterList(newRparams);
1874  newR->SetFactory("R", manager.GetFactory("R"));
1875  manager.SetFactory("R", newR);
1876  }
1877  }
1878 }
1879 
1880 // =====================================================================================================
1881 // =========================================== Nullspace ===============================================
1882 // =====================================================================================================
1883 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1886  int /* levelID */, std::vector<keep_pair>& /* keeps */, RCP<Factory>& nullSpaceFactory) const {
1887  // Nullspace
1888  RCP<Factory> nullSpace = rcp(new NullspaceFactory());
1889 
1890  bool have_userNS = false;
1891  if (paramList.isParameter("Nullspace") && !paramList.get<RCP<MultiVector> >("Nullspace").is_null())
1892  have_userNS = true;
1893 
1894  if (!have_userNS) {
1895  ParameterList newNullparams;
1896  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "nullspace: calculate rotations", bool, newNullparams);
1897  nullSpace->SetParameterList(newNullparams);
1898  nullSpace->SetFactory("Nullspace", manager.GetFactory("Ptent"));
1899  manager.SetFactory("Nullspace", nullSpace);
1900  }
1901  nullSpaceFactory = nullSpace;
1902 
1903  if (paramList.isParameter("restriction: scale nullspace") && paramList.get<bool>("restriction: scale nullspace")) {
1904  RCP<ScaledNullspaceFactory> scaledNSfactory = rcp(new ScaledNullspaceFactory());
1905  scaledNSfactory->SetFactory("Nullspace", nullSpaceFactory);
1906  manager.SetFactory("Scaled Nullspace", scaledNSfactory);
1907  }
1908 }
1909 
1910 // =====================================================================================================
1911 // ================================= Algorithm: SemiCoarsening =========================================
1912 // =====================================================================================================
1913 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1916  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
1917  // === Semi-coarsening ===
1918  RCP<Factory> semicoarsenFactory = Teuchos::null;
1919  if (paramList.isParameter("semicoarsen: number of levels") &&
1920  paramList.get<int>("semicoarsen: number of levels") > 0) {
1921  ParameterList togglePParams;
1922  ParameterList semicoarsenPParams;
1923  ParameterList linedetectionParams;
1924  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: number of levels", int, togglePParams);
1925  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: coarsen rate", int, semicoarsenPParams);
1926  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: piecewise constant", bool, semicoarsenPParams);
1927  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: piecewise linear", bool, semicoarsenPParams);
1928  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: calculate nonsym restriction", bool, semicoarsenPParams);
1929  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "linedetection: orientation", std::string, linedetectionParams);
1930  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "linedetection: num layers", int, linedetectionParams);
1931 
1933  RCP<LineDetectionFactory> linedetectionFactory = rcp(new LineDetectionFactory());
1934  RCP<TogglePFactory> togglePFactory = rcp(new TogglePFactory());
1935 
1936  linedetectionFactory->SetParameterList(linedetectionParams);
1937  semicoarsenFactory->SetParameterList(semicoarsenPParams);
1938  togglePFactory->SetParameterList(togglePParams);
1939 
1940  togglePFactory->AddCoarseNullspaceFactory(semicoarsenFactory);
1941  togglePFactory->AddProlongatorFactory(semicoarsenFactory);
1942  togglePFactory->AddPtentFactory(semicoarsenFactory);
1943  togglePFactory->AddCoarseNullspaceFactory(manager.GetFactory("Ptent"));
1944  togglePFactory->AddProlongatorFactory(manager.GetFactory("P"));
1945  togglePFactory->AddPtentFactory(manager.GetFactory("Ptent"));
1946 
1947  manager.SetFactory("CoarseNumZLayers", linedetectionFactory);
1948  manager.SetFactory("LineDetection_Layers", linedetectionFactory);
1949  manager.SetFactory("LineDetection_VertLineIds", linedetectionFactory);
1950 
1951  manager.SetFactory("P", togglePFactory);
1952  manager.SetFactory("Ptent", togglePFactory);
1953  manager.SetFactory("Nullspace", togglePFactory);
1954  }
1955 
1956  if (paramList.isParameter("semicoarsen: number of levels")) {
1957  auto tf = rcp(new ToggleCoordinatesTransferFactory());
1958  tf->SetFactory("Chosen P", manager.GetFactory("P"));
1959  tf->AddCoordTransferFactory(semicoarsenFactory);
1960 
1962  coords->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1963  coords->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1964  tf->AddCoordTransferFactory(coords);
1965  manager.SetFactory("Coordinates", tf);
1966  }
1967 }
1968 
1969 // =====================================================================================================
1970 // ================================== Algorithm: P-Coarsening ==========================================
1971 // =====================================================================================================
1972 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1975  int levelID, std::vector<keep_pair>& keeps) const {
1976 #ifdef HAVE_MUELU_INTREPID2
1977  // This only makes sense to invoke from the default list.
1978  if (defaultList.isParameter("pcoarsen: schedule") && defaultList.isParameter("pcoarsen: element")) {
1979  // P-Coarsening by schedule (new interface)
1980  // NOTE: levelID represents the *coarse* level in this case
1981  auto pcoarsen_schedule = Teuchos::getArrayFromStringParameter<int>(defaultList, "pcoarsen: schedule");
1982  auto pcoarsen_element = defaultList.get<std::string>("pcoarsen: element");
1983 
1984  if (levelID >= (int)pcoarsen_schedule.size()) {
1985  // Past the p-coarsening levels, we do Smoothed Aggregation
1986  // NOTE: We should probably consider allowing other options past p-coarsening
1987  UpdateFactoryManager_SA(paramList, defaultList, manager, levelID, keeps);
1988 
1989  } else {
1990  // P-Coarsening
1991  ParameterList Pparams;
1992  auto P = rcp(new IntrepidPCoarsenFactory());
1993  std::string lo = pcoarsen_element + std::to_string(pcoarsen_schedule[levelID]);
1994  std::string hi = (levelID ? pcoarsen_element + std::to_string(pcoarsen_schedule[levelID - 1]) : lo);
1995  Pparams.set("pcoarsen: hi basis", hi);
1996  Pparams.set("pcoarsen: lo basis", lo);
1997  P->SetParameterList(Pparams);
1998  manager.SetFactory("P", P);
1999 
2000  // Add special nullspace handling
2001  rcp_dynamic_cast<Factory>(manager.GetFactoryNonConst("Nullspace"))->SetFactory("Nullspace", manager.GetFactory("P"));
2002  }
2003 
2004  } else {
2005  // P-Coarsening by manual specification (old interface)
2006  ParameterList Pparams;
2007  auto P = rcp(new IntrepidPCoarsenFactory());
2008  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "pcoarsen: hi basis", std::string, Pparams);
2009  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "pcoarsen: lo basis", std::string, Pparams);
2010  P->SetParameterList(Pparams);
2011  manager.SetFactory("P", P);
2012 
2013  // Add special nullspace handling
2014  rcp_dynamic_cast<Factory>(manager.GetFactoryNonConst("Nullspace"))->SetFactory("Nullspace", manager.GetFactory("P"));
2015  }
2016 
2017 #endif
2018 }
2019 
2020 // =====================================================================================================
2021 // ============================== Algorithm: Smoothed Aggregation ======================================
2022 // =====================================================================================================
2023 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2025  UpdateFactoryManager_SA(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& keeps) const {
2026  // Smoothed aggregation
2027  RCP<Factory> P = rcp(new SaPFactory());
2028  ParameterList Pparams;
2029  if (paramList.isSublist("matrixmatrix: kernel params"))
2030  Pparams.sublist("matrixmatrix: kernel params", false) = paramList.sublist("matrixmatrix: kernel params");
2031  if (defaultList.isSublist("matrixmatrix: kernel params"))
2032  Pparams.sublist("matrixmatrix: kernel params", false) = defaultList.sublist("matrixmatrix: kernel params");
2033  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: damping factor", double, Pparams);
2034  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: calculate eigenvalue estimate", bool, Pparams);
2035  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: max eigenvalue", double, Pparams);
2036  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: eigenvalue estimate num iterations", int, Pparams);
2037  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: use rowsumabs diagonal scaling", bool, Pparams);
2038  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs diagonal replacement tolerance", double, Pparams);
2039  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs diagonal replacement value", double, Pparams);
2040  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs use automatic diagonal tolerance", bool, Pparams);
2041  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: enforce constraints", bool, Pparams);
2042  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: eigen-analysis type", std::string, Pparams);
2043  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: calculate qr", bool, Pparams);
2044 
2045  P->SetParameterList(Pparams);
2046 
2047  // Filtering
2048  MUELU_SET_VAR_2LIST(paramList, defaultList, "sa: use filtered matrix", bool, useFiltering);
2049  if (useFiltering) {
2050  // NOTE: Here, non-Kokkos and Kokkos versions diverge in the way the
2051  // dependency tree is setup. The Kokkos version has merged the the
2052  // FilteredAFactory into the CoalesceDropFactory.
2053  if (!useKokkos_) {
2054  RCP<Factory> filterFactory = rcp(new FilteredAFactory());
2055 
2056  ParameterList fParams;
2057  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use lumping", bool, fParams);
2058  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse graph", bool, fParams);
2059  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse eigenvalue", bool, fParams);
2060  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use root stencil", bool, fParams);
2061  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: Dirichlet threshold", double, fParams);
2062  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use spread lumping", bool, fParams);
2063  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom growth factor", double, fParams);
2064  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom cap", double, fParams);
2065  filterFactory->SetParameterList(fParams);
2066  filterFactory->SetFactory("Graph", manager.GetFactory("Graph"));
2067  filterFactory->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
2068  filterFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
2069  // I'm not sure why we need this line. See comments for DofsPerNode for UncoupledAggregation above
2070  filterFactory->SetFactory("Filtering", manager.GetFactory("Graph"));
2071 
2072  P->SetFactory("A", filterFactory);
2073 
2074  } else {
2075  P->SetFactory("A", manager.GetFactory("Graph"));
2076  }
2077  }
2078 
2079  P->SetFactory("P", manager.GetFactory("Ptent"));
2080  manager.SetFactory("P", P);
2081 
2082  bool filteringChangesMatrix = useFiltering && !MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol", double, 0);
2083  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
2084  if (reuseType == "tP" && !filteringChangesMatrix)
2085  keeps.push_back(keep_pair("AP reuse data", P.get()));
2086 }
2087 
2088 // =====================================================================================================
2089 // =============================== Algorithm: Energy Minimization ======================================
2090 // =====================================================================================================
2091 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2094  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
2095  MUELU_SET_VAR_2LIST(paramList, defaultList, "emin: pattern", std::string, patternType);
2096  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
2097  TEUCHOS_TEST_FOR_EXCEPTION(patternType != "AkPtent", Exceptions::InvalidArgument,
2098  "Invalid pattern name: \"" << patternType << "\". Valid options: \"AkPtent\"");
2099  // Pattern
2100  auto patternFactory = rcp(new PatternFactory());
2101  ParameterList patternParams;
2102  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: pattern order", int, patternParams);
2103  patternFactory->SetParameterList(patternParams);
2104  patternFactory->SetFactory("P", manager.GetFactory("Ptent"));
2105 
2106  // Filtering
2107  MUELU_SET_VAR_2LIST(paramList, defaultList, "emin: use filtered matrix", bool, useFiltering);
2108  if (useFiltering) {
2109  // NOTE: Here, non-Kokkos and Kokkos versions diverge in the way the
2110  // dependency tree is setup. The Kokkos version has merged the the
2111  // FilteredAFactory into the CoalesceDropFactory.
2112  if (!useKokkos_) {
2113  RCP<Factory> filterFactory = rcp(new FilteredAFactory());
2114 
2115  ParameterList fParams;
2116  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use lumping", bool, fParams);
2117  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse graph", bool, fParams);
2118  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse eigenvalue", bool, fParams);
2119  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use root stencil", bool, fParams);
2120  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: Dirichlet threshold", double, fParams);
2121  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use spread lumping", bool, fParams);
2122  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom growth factor", double, fParams);
2123  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom cap", double, fParams);
2124  filterFactory->SetParameterList(fParams);
2125  filterFactory->SetFactory("Graph", manager.GetFactory("Graph"));
2126  filterFactory->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
2127  filterFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
2128  // I'm not sure why we need this line. See comments for DofsPerNode for UncoupledAggregation above
2129  filterFactory->SetFactory("Filtering", manager.GetFactory("Graph"));
2130 
2131  patternFactory->SetFactory("A", filterFactory);
2132 
2133  } else {
2134  patternFactory->SetFactory("A", manager.GetFactory("Graph"));
2135  }
2136  }
2137 
2138  manager.SetFactory("Ppattern", patternFactory);
2139 
2140  // Constraint
2141  auto constraintFactory = rcp(new ConstraintFactory());
2142  constraintFactory->SetFactory("Ppattern", manager.GetFactory("Ppattern"));
2143  constraintFactory->SetFactory("CoarseNullspace", manager.GetFactory("Ptent"));
2144  manager.SetFactory("Constraint", constraintFactory);
2145 
2146  // Energy minimization
2147  ParameterList Pparams;
2148  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: num iterations", int, Pparams);
2149  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: iterative method", std::string, Pparams);
2150  if (reuseType == "emin") {
2151  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: num reuse iterations", int, Pparams);
2152  Pparams.set("Keep P0", true);
2153  Pparams.set("Keep Constraint0", true);
2154  }
2155 
2156  // Emin Factory
2157  auto P = rcp(new EminPFactory());
2158  P->SetParameterList(Pparams);
2159  P->SetFactory("P", manager.GetFactory("Ptent"));
2160  P->SetFactory("Constraint", manager.GetFactory("Constraint"));
2161  manager.SetFactory("P", P);
2162 }
2163 
2164 // =====================================================================================================
2165 // ================================= Algorithm: Petrov-Galerkin ========================================
2166 // =====================================================================================================
2167 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2169  UpdateFactoryManager_PG(ParameterList& /* paramList */, const ParameterList& /* defaultList */, FactoryManager& manager,
2170  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
2171  TEUCHOS_TEST_FOR_EXCEPTION(this->implicitTranspose_, Exceptions::RuntimeError,
2172  "Implicit transpose not supported with Petrov-Galerkin smoothed transfer operators: Set \"transpose: use implicit\" to false!\n"
2173  "Petrov-Galerkin transfer operator smoothing for non-symmetric problems requires a separate handling of the restriction operator which "
2174  "does not allow the usage of implicit transpose easily.");
2175 
2176  // Petrov-Galerkin
2177  auto P = rcp(new PgPFactory());
2178  P->SetFactory("P", manager.GetFactory("Ptent"));
2179  manager.SetFactory("P", P);
2180 }
2181 
2182 // =====================================================================================================
2183 // ================================= Algorithm: Replicate ========================================
2184 // =====================================================================================================
2185 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2187  UpdateFactoryManager_Replicate(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& keeps) const {
2189 
2190  ParameterList Pparams;
2191  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "replicate: npdes", int, Pparams);
2192 
2193  P->SetParameterList(Pparams);
2194  manager.SetFactory("P", P);
2195 }
2196 
2197 // =====================================================================================================
2198 // ====================================== Algorithm: Combine ============================================
2199 // =====================================================================================================
2200 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2202  UpdateFactoryManager_Combine(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& keeps) const {
2204 
2205  ParameterList Pparams;
2206  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "combine: numBlks", int, Pparams);
2207 
2208  P->SetParameterList(Pparams);
2209  manager.SetFactory("P", P);
2210 }
2211 
2212 // =====================================================================================================
2213 // ====================================== Algorithm: Matlab ============================================
2214 // =====================================================================================================
2215 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2217  UpdateFactoryManager_Matlab(ParameterList& paramList, const ParameterList& /* defaultList */, FactoryManager& manager,
2218  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
2219 #ifdef HAVE_MUELU_MATLAB
2220  ParameterList Pparams = paramList.sublist("transfer: params");
2221  auto P = rcp(new TwoLevelMatlabFactory());
2222  P->SetParameterList(Pparams);
2223  P->SetFactory("P", manager.GetFactory("Ptent"));
2224  manager.SetFactory("P", P);
2225 #else
2226  (void)paramList;
2227  (void)manager;
2228 #endif
2229 }
2230 
2231 #undef MUELU_SET_VAR_2LIST
2232 #undef MUELU_TEST_AND_SET_VAR
2233 #undef MUELU_TEST_AND_SET_PARAM_2LIST
2234 #undef MUELU_TEST_PARAM_2LIST
2235 #undef MUELU_KOKKOS_FACTORY
2236 
2237 size_t LevenshteinDistance(const char* s, size_t len_s, const char* t, size_t len_t);
2238 
2239 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2241  ParameterList paramList = constParamList;
2242  const ParameterList& validList = *MasterList::List();
2243  // Validate up to maxLevels level specific parameter sublists
2244  const int maxLevels = 100;
2245 
2246  // Extract level specific list
2247  std::vector<ParameterList> paramLists;
2248  for (int levelID = 0; levelID < maxLevels; levelID++) {
2249  std::string sublistName = "level " + toString(levelID);
2250  if (paramList.isSublist(sublistName)) {
2251  paramLists.push_back(paramList.sublist(sublistName));
2252  // paramLists.back().setName(sublistName);
2253  paramList.remove(sublistName);
2254  }
2255  }
2256  paramLists.push_back(paramList);
2257  // paramLists.back().setName("main");
2258 #ifdef HAVE_MUELU_MATLAB
2259  // If Muemex is supported, hide custom level variables from validator by removing them from paramList's sublists
2260  for (size_t i = 0; i < paramLists.size(); i++) {
2261  std::vector<std::string> customVars; // list of names (keys) to be removed from list
2262 
2263  for (Teuchos::ParameterList::ConstIterator it = paramLists[i].begin(); it != paramLists[i].end(); it++) {
2264  std::string paramName = paramLists[i].name(it);
2265 
2266  if (IsParamMuemexVariable(paramName))
2267  customVars.push_back(paramName);
2268  }
2269 
2270  // Remove the keys
2271  for (size_t j = 0; j < customVars.size(); j++)
2272  paramLists[i].remove(customVars[j], false);
2273  }
2274 #endif
2275 
2276  const int maxDepth = 0;
2277  for (size_t i = 0; i < paramLists.size(); i++) {
2278  // validate every sublist
2279  try {
2280  paramLists[i].validateParameters(validList, maxDepth);
2281 
2282  } catch (const Teuchos::Exceptions::InvalidParameterName& e) {
2283  std::string eString = e.what();
2284 
2285  // Parse name from: <Error, the parameter {name="smoothe: type",...>
2286  size_t nameStart = eString.find_first_of('"') + 1;
2287  size_t nameEnd = eString.find_first_of('"', nameStart);
2288  std::string name = eString.substr(nameStart, nameEnd - nameStart);
2289 
2290  size_t bestScore = 100;
2291  std::string bestName = "";
2292  for (ParameterList::ConstIterator it = validList.begin(); it != validList.end(); it++) {
2293  const std::string& pName = validList.name(it);
2294  this->GetOStream(Runtime1) << "| " << pName;
2295  size_t score = LevenshteinDistance(name.c_str(), name.length(), pName.c_str(), pName.length());
2296  this->GetOStream(Runtime1) << " -> " << score << std::endl;
2297  if (score < bestScore) {
2298  bestScore = score;
2299  bestName = pName;
2300  }
2301  }
2302  if (bestScore < 10 && bestName != "") {
2304  eString << "The parameter name \"" + name + "\" is not valid. Did you mean \"" + bestName << "\"?\n");
2305 
2306  } else {
2308  eString << "The parameter name \"" + name + "\" is not valid.\n");
2309  }
2310  }
2311  }
2312 }
2313 
2314 // =====================================================================================================
2315 // ==================================== FACTORY interpreter ============================================
2316 // =====================================================================================================
2317 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2319  SetFactoryParameterList(const ParameterList& constParamList) {
2320  // Create a non const copy of the parameter list
2321  // Working with a modifiable list is much much easier than with original one
2322  ParameterList paramList = constParamList;
2323 
2324  // Parameter List Parsing:
2325  // ---------
2326  // <ParameterList name="MueLu">
2327  // <ParameterList name="Matrix">
2328  // </ParameterList>
2329  if (paramList.isSublist("Matrix")) {
2330  blockSize_ = paramList.sublist("Matrix").get<int>("PDE equations", MasterList::getDefault<int>("number of equations"));
2331  dofOffset_ = paramList.sublist("Matrix").get<GlobalOrdinal>("DOF offset", 0); // undocumented parameter allowing to define a DOF offset of the global dofs of an operator (defaul = 0)
2332  }
2333 
2334  // create new FactoryFactory object if necessary
2335  if (factFact_ == Teuchos::null)
2336  factFact_ = Teuchos::rcp(new FactoryFactory());
2337 
2338  // Parameter List Parsing:
2339  // ---------
2340  // <ParameterList name="MueLu">
2341  // <ParameterList name="Factories"> <== call BuildFactoryMap() on this parameter list
2342  // ...
2343  // </ParameterList>
2344  // </ParameterList>
2345  FactoryMap factoryMap;
2346  FactoryManagerMap factoryManagers;
2347  if (paramList.isSublist("Factories"))
2348  this->BuildFactoryMap(paramList.sublist("Factories"), factoryMap, factoryMap, factoryManagers);
2349 
2350  // Parameter List Parsing:
2351  // ---------
2352  // <ParameterList name="MueLu">
2353  // <ParameterList name="Hierarchy">
2354  // <Parameter name="verbose" type="string" value="Warnings"/> <== get
2355  // <Parameter name="numDesiredLevel" type="int" value="10"/> <== get
2356  //
2357  // <ParameterList name="firstLevel"> <== parse first args and call BuildFactoryMap() on the rest of this parameter list
2358  // ...
2359  // </ParameterList>
2360  // </ParameterList>
2361  // </ParameterList>
2362  if (paramList.isSublist("Hierarchy")) {
2363  ParameterList hieraList = paramList.sublist("Hierarchy"); // copy because list temporally modified (remove 'id')
2364 
2365  // Get hierarchy options
2366  if (hieraList.isParameter("max levels")) {
2367  this->numDesiredLevel_ = hieraList.get<int>("max levels");
2368  hieraList.remove("max levels");
2369  }
2370 
2371  if (hieraList.isParameter("coarse: max size")) {
2372  this->maxCoarseSize_ = hieraList.get<int>("coarse: max size");
2373  hieraList.remove("coarse: max size");
2374  }
2375 
2376  if (hieraList.isParameter("repartition: rebalance P and R")) {
2377  this->doPRrebalance_ = hieraList.get<bool>("repartition: rebalance P and R");
2378  hieraList.remove("repartition: rebalance P and R");
2379  }
2380 
2381  if (hieraList.isParameter("transpose: use implicit")) {
2382  this->implicitTranspose_ = hieraList.get<bool>("transpose: use implicit");
2383  hieraList.remove("transpose: use implicit");
2384  }
2385 
2386  if (hieraList.isParameter("fuse prolongation and update")) {
2387  this->fuseProlongationAndUpdate_ = hieraList.get<bool>("fuse prolongation and update");
2388  hieraList.remove("fuse prolongation and update");
2389  }
2390 
2391  if (hieraList.isParameter("nullspace: suppress dimension check")) {
2392  this->suppressNullspaceDimensionCheck_ = hieraList.get<bool>("nullspace: suppress dimension check");
2393  hieraList.remove("nullspace: suppress dimension check");
2394  }
2395 
2396  if (hieraList.isParameter("number of vectors")) {
2397  this->sizeOfMultiVectors_ = hieraList.get<int>("number of vectors");
2398  hieraList.remove("number of vectors");
2399  }
2400 
2401  if (hieraList.isSublist("matvec params"))
2402  this->matvecParams_ = Teuchos::parameterList(hieraList.sublist("matvec params"));
2403 
2404  if (hieraList.isParameter("coarse grid correction scaling factor")) {
2405  this->scalingFactor_ = hieraList.get<double>("coarse grid correction scaling factor");
2406  hieraList.remove("coarse grid correction scaling factor");
2407  }
2408 
2409  // Translate cycle type parameter
2410  if (hieraList.isParameter("cycle type")) {
2411  std::map<std::string, CycleType> cycleMap;
2412  cycleMap["V"] = VCYCLE;
2413  cycleMap["W"] = WCYCLE;
2414 
2415  std::string cycleType = hieraList.get<std::string>("cycle type");
2416  TEUCHOS_TEST_FOR_EXCEPTION(cycleMap.count(cycleType) == 0, Exceptions::RuntimeError, "Invalid cycle type: \"" << cycleType << "\"");
2417  this->Cycle_ = cycleMap[cycleType];
2418  }
2419 
2420  if (hieraList.isParameter("W cycle start level")) {
2421  this->WCycleStartLevel_ = hieraList.get<int>("W cycle start level");
2422  }
2423 
2424  if (hieraList.isParameter("verbosity")) {
2425  std::string vl = hieraList.get<std::string>("verbosity");
2426  hieraList.remove("verbosity");
2427  this->verbosity_ = toVerbLevel(vl);
2428  }
2429 
2430  if (hieraList.isParameter("output filename"))
2431  VerboseObject::SetMueLuOFileStream(hieraList.get<std::string>("output filename"));
2432 
2433  if (hieraList.isParameter("dependencyOutputLevel"))
2434  this->graphOutputLevel_ = hieraList.get<int>("dependencyOutputLevel");
2435 
2436  // Check for the reuse case
2437  if (hieraList.isParameter("reuse"))
2439 
2440  if (hieraList.isSublist("DataToWrite")) {
2441  // TODO We should be able to specify any data. If it exists, write it.
2442  // TODO This would requires something like std::set<dataName, Array<int> >
2443  ParameterList foo = hieraList.sublist("DataToWrite");
2444  std::string dataName = "Matrices";
2445  if (foo.isParameter(dataName))
2446  this->matricesToPrint_["A"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2447  dataName = "Prolongators";
2448  if (foo.isParameter(dataName))
2449  this->matricesToPrint_["P"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2450  dataName = "Restrictors";
2451  if (foo.isParameter(dataName))
2452  this->matricesToPrint_["R"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2453  dataName = "D0";
2454  if (foo.isParameter(dataName))
2455  this->matricesToPrint_["D0"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2456  }
2457 
2458  // Get level configuration
2459  for (ParameterList::ConstIterator param = hieraList.begin(); param != hieraList.end(); ++param) {
2460  const std::string& paramName = hieraList.name(param);
2461 
2462  if (paramName != "DataToWrite" && hieraList.isSublist(paramName)) {
2463  ParameterList levelList = hieraList.sublist(paramName); // copy because list temporally modified (remove 'id')
2464 
2465  int startLevel = 0;
2466  if (levelList.isParameter("startLevel")) {
2467  startLevel = levelList.get<int>("startLevel");
2468  levelList.remove("startLevel");
2469  }
2470  int numDesiredLevel = 1;
2471  if (levelList.isParameter("numDesiredLevel")) {
2472  numDesiredLevel = levelList.get<int>("numDesiredLevel");
2473  levelList.remove("numDesiredLevel");
2474  }
2475 
2476  // Parameter List Parsing:
2477  // ---------
2478  // <ParameterList name="firstLevel">
2479  // <Parameter name="startLevel" type="int" value="0"/>
2480  // <Parameter name="numDesiredLevel" type="int" value="1"/>
2481  // <Parameter name="verbose" type="string" value="Warnings"/>
2482  //
2483  // [] <== call BuildFactoryMap() on the rest of the parameter list
2484  //
2485  // </ParameterList>
2486  FactoryMap levelFactoryMap;
2487  BuildFactoryMap(levelList, factoryMap, levelFactoryMap, factoryManagers);
2488 
2489  RCP<FactoryManager> m = rcp(new FactoryManager(levelFactoryMap));
2490  if (hieraList.isParameter("use kokkos refactor"))
2491  m->SetKokkosRefactor(hieraList.get<bool>("use kokkos refactor"));
2492 
2493  if (startLevel >= 0)
2494  this->AddFactoryManager(startLevel, numDesiredLevel, m);
2495  else
2496  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::ParameterListInterpreter():: invalid level id");
2497  } /* TODO: else { } */
2498  }
2499  }
2500 }
2501 
2502 // TODO: static?
2536 
2588 
2625 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2627  BuildFactoryMap(const ParameterList& paramList, const FactoryMap& factoryMapIn, FactoryMap& factoryMapOut, FactoryManagerMap& factoryManagers) const {
2628  for (ParameterList::ConstIterator param = paramList.begin(); param != paramList.end(); ++param) {
2629  const std::string& paramName = paramList.name(param); //< paramName contains the user chosen factory name (e.g., "smootherFact1")
2630  const Teuchos::ParameterEntry& paramValue = paramList.entry(param); //< for factories, paramValue should be either a list or just a MueLu Factory (e.g., TrilinosSmoother)
2631 
2632  // TODO: do not allow name of existing MueLu classes (can be tested using FactoryFactory)
2633 
2634  if (paramValue.isList()) {
2635  ParameterList paramList1 = Teuchos::getValue<ParameterList>(paramValue);
2636  if (paramList1.isParameter("factory")) { // default: just a factory definition
2637  // New Factory is a sublist with internal parameters and/or data dependencies
2638  TEUCHOS_TEST_FOR_EXCEPTION(paramList1.isParameter("dependency for") == true, Exceptions::RuntimeError,
2639  "MueLu::ParameterListInterpreter(): It seems that in the parameter lists for defining " << paramName << " there is both a 'factory' and 'dependency for' parameter. This is not allowed. Please remove the 'dependency for' parameter.");
2640 
2641  factoryMapOut[paramName] = factFact_->BuildFactory(paramValue, factoryMapIn, factoryManagers);
2642 
2643  } else if (paramList1.isParameter("dependency for")) { // add more data dependencies to existing factory
2644  TEUCHOS_TEST_FOR_EXCEPTION(paramList1.isParameter("factory") == true, Exceptions::RuntimeError,
2645  "MueLu::ParameterListInterpreter(): It seems that in the parameter lists for defining " << paramName << " there is both a 'factory' and 'dependency for' parameter. This is not allowed.");
2646 
2647  std::string factoryName = paramList1.get<std::string>("dependency for");
2648 
2649  RCP<const FactoryBase> factbase = factoryMapIn.find(factoryName /*paramName*/)->second; // access previously defined factory
2650  TEUCHOS_TEST_FOR_EXCEPTION(factbase.is_null() == true, Exceptions::RuntimeError,
2651  "MueLu::ParameterListInterpreter(): could not find factory " + factoryName + " in factory map. Did you define it before?");
2652 
2653  RCP<const Factory> factoryconst = Teuchos::rcp_dynamic_cast<const Factory>(factbase);
2654  RCP<Factory> factory = Teuchos::rcp_const_cast<Factory>(factoryconst);
2655 
2656  // Read the RCP<Factory> parameters of the class T
2657  RCP<const ParameterList> validParamList = factory->GetValidParameterList();
2658  for (ParameterList::ConstIterator vparam = validParamList->begin(); vparam != validParamList->end(); ++vparam) {
2659  const std::string& pName = validParamList->name(vparam);
2660 
2661  if (!paramList1.isParameter(pName)) {
2662  // Ignore unknown parameters
2663  continue;
2664  }
2665 
2666  if (validParamList->isType<RCP<const FactoryBase> >(pName)) {
2667  // Generate or get factory described by pName and set dependency
2668  RCP<const FactoryBase> generatingFact = factFact_->BuildFactory(paramList1.getEntry(pName), factoryMapIn, factoryManagers);
2669  factory->SetFactory(pName, generatingFact.create_weak());
2670 
2671  } else if (validParamList->isType<RCP<const ParameterList> >(pName)) {
2672  if (pName == "ParameterList") {
2673  // NOTE: we cannot use
2674  // subList = sublist(rcpFromRef(paramList), pName)
2675  // here as that would result in sublist also being a reference to a temporary object.
2676  // The resulting dereferencing in the corresponding factory would then segfault
2677  RCP<const ParameterList> subList = Teuchos::sublist(rcp(new ParameterList(paramList1)), pName);
2678  factory->SetParameter(pName, ParameterEntry(subList));
2679  }
2680  } else {
2681  factory->SetParameter(pName, paramList1.getEntry(pName));
2682  }
2683  }
2684 
2685  } else if (paramList1.isParameter("group")) { // definitiion of a factory group (for a factory manager)
2686  // Define a new (sub) FactoryManager
2687  std::string groupType = paramList1.get<std::string>("group");
2688  TEUCHOS_TEST_FOR_EXCEPTION(groupType != "FactoryManager", Exceptions::RuntimeError,
2689  "group must be of type \"FactoryManager\".");
2690 
2691  ParameterList groupList = paramList1; // copy because list temporally modified (remove 'id')
2692  groupList.remove("group");
2693 
2694  bool setKokkosRefactor = false;
2695  bool kokkosRefactor = useKokkos_;
2696  if (groupList.isParameter("use kokkos refactor")) {
2697  kokkosRefactor = groupList.get<bool>("use kokkos refactor");
2698  groupList.remove("use kokkos refactor");
2699  setKokkosRefactor = true;
2700  }
2701 
2702  FactoryMap groupFactoryMap;
2703  BuildFactoryMap(groupList, factoryMapIn, groupFactoryMap, factoryManagers);
2704 
2705  // do not store groupFactoryMap in factoryMapOut
2706  // Create a factory manager object from groupFactoryMap
2707  RCP<FactoryManager> m = rcp(new FactoryManager(groupFactoryMap));
2708  if (setKokkosRefactor)
2709  m->SetKokkosRefactor(kokkosRefactor);
2710  factoryManagers[paramName] = m;
2711 
2712  } else {
2713  this->GetOStream(Warnings0) << "Could not interpret parameter list " << paramList1 << std::endl;
2715  "XML Parameter list must either be of type \"factory\" or of type \"group\".");
2716  }
2717  } else {
2718  // default: just a factory (no parameter list)
2719  factoryMapOut[paramName] = factFact_->BuildFactory(paramValue, factoryMapIn, factoryManagers);
2720  }
2721  }
2722 }
2723 
2724 // =====================================================================================================
2725 // ======================================= MISC functions ==============================================
2726 // =====================================================================================================
2727 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2729  try {
2730  Matrix& A = dynamic_cast<Matrix&>(Op);
2731  if (A.IsFixedBlockSizeSet() && (A.GetFixedBlockSize() != blockSize_))
2732  this->GetOStream(Warnings0) << "Setting matrix block size to " << blockSize_ << " (value of the parameter in the list) "
2733  << "instead of " << A.GetFixedBlockSize() << " (provided matrix)." << std::endl
2734  << "You may want to check \"number of equations\" (or \"PDE equations\" for factory style list) parameter." << std::endl;
2735 
2736  A.SetFixedBlockSize(blockSize_, dofOffset_);
2737 
2738 #ifdef HAVE_MUELU_DEBUG
2739  MatrixUtils::checkLocalRowMapMatchesColMap(A);
2740 #endif // HAVE_MUELU_DEBUG
2741 
2742  } catch (std::bad_cast&) {
2743  this->GetOStream(Warnings0) << "Skipping setting block size as the operator is not a matrix" << std::endl;
2744  }
2745 }
2746 
2747 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2749  H.SetCycle(Cycle_);
2750  H.SetCycleStartLevel(WCycleStartLevel_);
2751  H.SetProlongatorScalingFactor(scalingFactor_);
2753 }
2754 
2755 static bool compare(const ParameterList& list1, const ParameterList& list2) {
2756  // First loop through and validate the parameters at this level.
2757  // In addition, we generate a list of sublists that we will search next
2758  for (ParameterList::ConstIterator it = list1.begin(); it != list1.end(); it++) {
2759  const std::string& name = it->first;
2760  const Teuchos::ParameterEntry& entry1 = it->second;
2761 
2762  const Teuchos::ParameterEntry* entry2 = list2.getEntryPtr(name);
2763  if (!entry2) // entry is not present in the second list
2764  return false;
2765  if (entry1.isList() && entry2->isList()) { // sublist check
2766  compare(Teuchos::getValue<ParameterList>(entry1), Teuchos::getValue<ParameterList>(*entry2));
2767  continue;
2768  }
2769  if (entry1.getAny(false) != entry2->getAny(false)) // entries have different types or different values
2770  return false;
2771  }
2772 
2773  return true;
2774 }
2775 
2776 static inline bool areSame(const ParameterList& list1, const ParameterList& list2) {
2777  return compare(list1, list2) && compare(list2, list1);
2778 }
2779 
2780 } // namespace MueLu
2781 
2782 #define MUELU_PARAMETERLISTINTERPRETER_SHORT
2783 #endif /* MUELU_PARAMETERLISTINTERPRETER_DEF_HPP */
Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
virtual ~ParameterListInterpreter()
Destructor.
const std::string & name() const
This class specifies the default factory that should generate some data on a Level if the data does n...
void SetupHierarchy(Hierarchy &H) const
Call the SetupHierarchy routine from the HiearchyManager object.
ParameterList & setEntry(const std::string &name, U &&entry)
static VerbLevel GetDefaultVerbLevel()
Get the default (global) verbosity level.
ConstIterator end() const
void UpdateFactoryManager_Replicate(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
void UpdateFactoryManager_Reitzinger(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
Factory for determing the number of partitions for rebalancing.
#define MUELU_KOKKOS_FACTORY_NO_DECL(varName, oldFactory, newFactory)
Factory for generating coarse level map. Used by TentativePFactory.
bool is_null(const boost::shared_ptr< T > &p)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Factory for building transfer operators based on coarsening in polynomial degree, following the Intre...
virtual void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Configuration.
Factory for building coarse grid matrices, when the matrix is of the form K+a*M. Useful when you want...
Class for generating an initial LocalOrdinal-type BlockNumber vector, based on an input paraemter for...
void UpdateFactoryManager(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
Factory that can generate other factories from.
RCP< T > create_weak() const
T & get(const std::string &name, T def_value)
Class that encapsulates external library smoothers.
size_t LevenshteinDistance(const char *s, size_t len_s, const char *t, size_t len_t)
void UpdateFactoryManager_Matlab(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
#define MUELU_TEST_AND_SET_VAR(paramList, paramName, paramType, varName)
void AddProlongatorFactory(const RCP< const FactoryBase > &factory)
Add a prolongator factory in the end of list of prolongator factories.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Factory for building permutation matrix that can be be used to shuffle data (matrices, vectors) among processes.
Class for transferring a vector of local ordinals from a finer level to a coarser one...
Factory for converting matrices to half precision operators.
void SetKokkosRefactor(const bool useKokkos)
Class for restricting a MultiVector from a finer to a coarser level.
One-liner description of what is happening.
void UpdateFactoryManager_LowPrecision(ParameterList &paramList, const ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
T * get() const
Ordinal numParams() const
void SetCycle(CycleType Cycle)
Supports VCYCLE and WCYCLE types.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define MUELU_KOKKOS_FACTORY(varName, oldFactory, newFactory)
Interface to Zoltan library.This interface provides access to partitioning methods in Zoltan...
bool IsParamMuemexVariable(const std::string &name)
Factory for creating a graph base on a given matrix.
void UpdateFactoryManager_SemiCoarsen(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Get factory associated with a particular data name.
static void DisableMultipleCheckGlobally()
std::map< std::string, RCP< FactoryManagerBase > > FactoryManagerMap
std::map< std::string, RCP< const FactoryBase > > FactoryMap
void UpdateFactoryManager_Nullspace(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps, RCP< Factory > &nullSpaceFactory) const
static void EnableTimerSync()
Factory for building tentative prolongator.
static void SetMueLuOFileStream(const std::string &filename)
#define TEST_MUTUALLY_EXCLUSIVE(arg1, arg2)
MsgType toVerbLevel(const std::string &verbLevelStr)
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
virtual void SetupOperator(Operator &A) const
Setup Operator object.
Prolongator factory performing semi-coarsening.
Factory for building restriction operators using a prolongator factory.
static RCP< Time > getNewTimer(const std::string &name)
ParameterEntry * getEntryPtr(const std::string &name)
void UpdateFactoryManager_PG(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
Teuchos::RCP< MueLu::FacadeClassFactory< Scalar, LocalOrdinal, GlobalOrdinal, Node > > facadeFact_
FacadeClass factory.
Additional warnings.
Prolongator factory that replicates &#39;Psubblock&#39; matrix to create new prolongator suitable for PDE sys...
bool isParameter(const std::string &name) const
void UpdateFactoryManager_Material(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
bool remove(std::string const &name, bool throwIfNotExists=true)
void UpdateFactoryManager_Emin(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
#define MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, paramName, paramType, listWrite)
static bool compare(const ParameterList &list1, const ParameterList &list2)
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
static CycleType GetDefaultCycle()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
void UpdateFactoryManager_Coordinates(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
MueLu::DefaultGlobalOrdinal GlobalOrdinal
bool isSublist(const std::string &name) const
void SetCycleStartLevel(int cycleStart)
Factory for interacting with Matlab.
Factory for interacting with Matlab.
VerbLevel GetVerbLevel() const
Get the verbosity level.
Factory for building line detection information.
params_t::ConstIterator ConstIterator
void SetFactoryParameterList(const Teuchos::ParameterList &paramList)
Factory interpreter stuff.
Factory to export aggregation info or visualize aggregates using VTK.
Prolongator factory which allows switching between two different prolongator strategies.
Interface to Zoltan2 library.This interface provides access to partitioning methods in Zoltan2...
AmalgamationFactory for subblocks of strided map based amalgamation data.
static const RCP< const NoFactory > getRCP()
Static Get() functions.
void UpdateFactoryManager_RAP(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
ConstIterator begin() const
Factory for building the constraint operator.
void UpdateFactoryManager_Restriction(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
void UpdateFactoryManager_BlockNumber(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
Applies permutation to grid transfer operators.
Prolongator factory that replicates &#39;Psubblock&#39; matrix to create new prolongator suitable for PDE sys...
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
ParameterList & setParameters(const ParameterList &source)
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Set Factory.
const RCP< FactoryBase > GetFactoryNonConst(const std::string &varName)
Get factory associated with a particular data name (NONCONST version)
Factory for generating a very special nullspace.
const ParameterEntry & entry(ConstIterator i) const
any & getAny(bool activeQry=true)
void UpdateFactoryManager_Aggregation_TentativeP(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
TransListIter iter
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
Factory for creating a graph based on a given matrix.
void BuildFactoryMap(const Teuchos::ParameterList &paramList, const FactoryMap &factoryMapIn, FactoryMap &factoryMapOut, FactoryManagerMap &factoryManagers) const
Interpret &quot;Factories&quot; sublist.
Class that encapsulates Matlab smoothers.
This class checks matrix properties of A on current level. This factory can be plugged in everywhere ...
Partitioning within a node onlyThis interface provides partitioning within a node.
Factory for generating F/C-splitting and a coarse level map. Used by ClassicalPFactory.
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameter list for Parameter list interpreter.
void UpdateFactoryManager_Combine(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
Class for transferring coordinates from a finer level to a coarser one.
void UpdateFactoryManager_SA(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
void SetVerbLevel(const VerbLevel verbLevel)
Set the verbosity level of this object.
Factory for building nonzero patterns for energy minimization.
static Teuchos::RCP< Teuchos::ParameterList > GetProblemSpecificList(std::string const &problemType)
Return default parameter settings for the specified problem type.
bool isType(const std::string &name) const
Factory for building tentative prolongator.
Class for transferring coordinates from a finer level to a coarser one.
Factory for building restriction operators.
Factory for building Energy Minimization prolongators.
static bool areSame(const ParameterList &list1, const ParameterList &list2)
Helper functions to compare two paramter lists.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Factory for creating a graph based on a given matrix.
#define TEST_MUTUALLY_EXCLUSIVE_S(arg1, arg2)
static int GetDefaultCycleStartLevel()
void Validate(const Teuchos::ParameterList &paramList) const
void UpdateFactoryManager_PCoarsen(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
void SetProlongatorScalingFactor(double scalingFactor)
Specify damping factor alpha such that x = x + alpha*P*c, where c is the coarse grid correction...
#define MUELU_TEST_PARAM_2LIST(paramList, defaultList, paramName, paramType, cmpValue)
Factory for building coarse matrices.
Exception throws to report errors in the internal logical of the program.
Factory for building filtered matrices using filtered graphs.
std::pair< std::string, const FactoryBase * > keep_pair
void UpdateFactoryManager_Repartition(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps, RCP< Factory > &nullSpaceFactory) const
Description of what is happening (more verbose)
Factory for building coarse matrices.
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
An factory which assigns each aggregate a quality estimate. Originally developed by Napov and Notay i...
void UpdateFactoryManager_Smoothers(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
ParameterEntry & getEntry(const std::string &name)
void AddPtentFactory(const RCP< const FactoryBase > &factory)
Add a tentative prolongator factory in the end of list of prolongator factories.
Factory for building Smoothed Aggregation prolongators.Input/output of SaPFactory
Factory for building uncoupled aggregates.
static Teuchos::RCP< const Teuchos::ParameterList > List()
Return a &quot;master&quot; list of all valid parameters and their default values.
Prolongator factory performing semi-coarsening.
void UpdateFactoryManager_CoarseSolvers(Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const
#define MUELU_SET_VAR_2LIST(paramList, defaultList, paramName, paramType, varName)
void SetEasyParameterList(const Teuchos::ParameterList &paramList)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void AddCoarseNullspaceFactory(const RCP< const FactoryBase > &factory)
Add a coarse nullspace factory in the end of list of coarse nullspace factories.
Factory for generating nullspace.
long ExtractNonSerializableData(const Teuchos::ParameterList &inList, Teuchos::ParameterList &serialList, Teuchos::ParameterList &nonSerialList)
Extract non-serializable data from level-specific sublists and move it to a separate parameter list...
Exception throws to report invalid user entry.
#define TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(throw_exception_test, Exception, msg)
bool is_null() const
void UpdateFactoryManager_LocalOrdinalTransfer(const std::string &VarName, const std::string &multigridAlgo, Teuchos::ParameterList &paramList, const Teuchos::ParameterList &defaultList, FactoryManager &manager, int levelID, std::vector< keep_pair > &keeps) const