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: lumping choice", std::string, dropParams);
1076  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom growth factor", double, dropParams);
1077  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom cap", double, dropParams);
1078  }
1079 
1080 #ifdef HAVE_MUELU_COALESCEDROP_ALLOW_OLD_PARAMETERS
1081  if (!dropParams.isParameter("aggregation: drop scheme") ||
1082  (dropParams.isParameter("aggregation: drop scheme") &&
1083  ((dropParams.get<std::string>("aggregation: drop scheme") != "point-wise") && (dropParams.get<std::string>("aggregation: drop scheme") != "cut-drop")))) {
1084  Teuchos::ParameterList dropParamsWithDefaults(dropParams);
1085 
1086 #define MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(paramList, paramName, paramType) \
1087  if (!paramList.isParameter(paramName)) { \
1088  paramList.set(paramName, MasterList::getDefault<paramType>(paramName)); \
1089  }
1090 
1091  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: drop scheme", std::string);
1092  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: strength-of-connection: matrix", std::string);
1093  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: strength-of-connection: measure", std::string);
1094  MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST(dropParamsWithDefaults, "aggregation: use blocking", bool);
1095 
1096 #undef MUELU_TEST_AND_SET_VAR_FROM_MASTERLIST
1097 
1098  // We are using the old style of dropping params
1099  TEUCHOS_TEST_FOR_EXCEPTION(dropParams.isParameter("aggregation: strength-of-connection: matrix") ||
1100  dropParams.isParameter("aggregation: strength-of-connection: measure") ||
1101  dropParams.isParameter("aggregation: use blocking"),
1103  "The inputs contain a mix of old and new dropping parameters:\n\n"
1104  << dropParams << "\n\nKeep in mind that defaults are set for old parameters, so this gets interpreted as\n\n"
1105  << dropParamsWithDefaults);
1106  }
1107 #endif
1108 
1109  if (!amalgFact.is_null())
1110  dropFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
1111 
1112  if (dropParams.isParameter("aggregation: drop scheme")) {
1113  std::string drop_scheme = dropParams.get<std::string>("aggregation: drop scheme");
1114  if (drop_scheme == "block diagonal colored signed classical")
1115  manager.SetFactory("Coloring Graph", dropFactory);
1116  if ((MUELU_TEST_PARAM_2LIST(dropParams, defaultList, "aggregation: use blocking", bool, true)) ||
1117  (drop_scheme.find("block diagonal") != std::string::npos || drop_scheme == "signed classical")) {
1118  if (levelID > 0)
1119  dropFactory->SetFactory("BlockNumber", this->GetFactoryManager(levelID - 1)->GetFactory("BlockNumber"));
1120  else
1121  dropFactory->SetFactory("BlockNumber", manager.GetFactory("BlockNumber"));
1122  }
1123  }
1124 
1125  dropFactory->SetParameterList(dropParams);
1126  }
1127  manager.SetFactory("Graph", dropFactory);
1128 
1129 // Aggregation scheme
1130 #ifndef HAVE_MUELU_MATLAB
1131  if (aggType == "matlab")
1132  throw std::runtime_error("Cannot use MATLAB aggregation - MueLu was not configured with MATLAB support.");
1133 #endif
1134  RCP<Factory> aggFactory;
1135  if (aggType == "uncoupled") {
1136  aggFactory = rcp(new UncoupledAggregationFactory());
1137  ParameterList aggParams;
1138  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: mode", std::string, aggParams);
1139  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: ordering", std::string, aggParams);
1140  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: min agg size", int, aggParams);
1141  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: max agg size", int, aggParams);
1142  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: max selected neighbors", int, aggParams);
1143  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: backend", std::string, aggParams);
1144  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: phase 1 algorithm", std::string, aggParams);
1145  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: deterministic", bool, aggParams);
1146  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: coloring algorithm", std::string, aggParams);
1147  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 1", bool, aggParams);
1148  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 2a", bool, aggParams);
1149  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 2b", bool, aggParams);
1150  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 3", bool, aggParams);
1151  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: match ML phase1", bool, aggParams);
1152  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: match ML phase2a", bool, aggParams);
1153  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: match ML phase2b", bool, aggParams);
1154  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: phase2a agg factor", double, aggParams);
1155  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: preserve Dirichlet points", bool, aggParams);
1156  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: error on nodes with no on-rank neighbors", bool, aggParams);
1157  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: phase3 avoid singletons", bool, aggParams);
1158  aggFactory->SetParameterList(aggParams);
1159  // make sure that the aggregation factory has all necessary data
1160  aggFactory->SetFactory("DofsPerNode", manager.GetFactory("Graph"));
1161  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
1162  // aggFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
1163 
1164  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: coloring algorithm", std::string, "mis2 aggregation") ||
1165  MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: coloring algorithm", std::string, "mis2 coarsening")) {
1166  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: symmetrize graph after dropping", bool, false))
1169  "MIS2 algorithms require the use of a symmetrized graph. Please set \"aggregation: symmetrize graph after dropping\" to \"true\".");
1170  }
1171  } else if (aggType == "brick") {
1172  aggFactory = rcp(new BrickAggregationFactory());
1173  ParameterList aggParams;
1174  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick x size", int, aggParams);
1175  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick y size", int, aggParams);
1176  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick z size", int, aggParams);
1177  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick x Dirichlet", bool, aggParams);
1178  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick y Dirichlet", bool, aggParams);
1179  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick z Dirichlet", bool, aggParams);
1180  aggFactory->SetParameterList(aggParams);
1181 
1182  // Unlike other factories, BrickAggregationFactory makes the Graph/DofsPerNode itself
1183  manager.SetFactory("Graph", aggFactory);
1184  manager.SetFactory("DofsPerNode", aggFactory);
1185  manager.SetFactory("Filtering", aggFactory);
1186  if (levelID > 1) {
1187  // We check for levelID > 0, as in the interpreter aggFactory for
1188  // levelID really corresponds to level 0. Managers are clunky, as they
1189  // contain factories for two different levels
1190  aggFactory->SetFactory("Coordinates", this->GetFactoryManager(levelID - 1)->GetFactory("Coordinates"));
1191  }
1192  } else if (aggType == "classical") {
1193  // Map and coloring
1194  RCP<Factory> mapFact = rcp(new ClassicalMapFactory());
1195  ParameterList mapParams;
1196  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: deterministic", bool, mapParams);
1197  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: coloring algorithm", std::string, mapParams);
1198 
1199  ParameterList tempParams;
1200  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop scheme", std::string, tempParams);
1201  std::string drop_algo = tempParams.get<std::string>("aggregation: drop scheme");
1202  if (drop_algo == "block diagonal colored signed classical") {
1203  mapParams.set("aggregation: coloring: use color graph", true);
1204  mapFact->SetFactory("Coloring Graph", manager.GetFactory("Coloring Graph"));
1205  }
1206  mapFact->SetParameterList(mapParams);
1207  mapFact->SetFactory("Graph", manager.GetFactory("Graph"));
1208  mapFact->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
1209 
1210  manager.SetFactory("FC Splitting", mapFact);
1211  manager.SetFactory("CoarseMap", mapFact);
1212 
1213  aggFactory = rcp(new ClassicalPFactory());
1214  ParameterList aggParams;
1215  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: classical scheme", std::string, aggParams);
1216  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop scheme", std::string, aggParams);
1217  aggFactory->SetParameterList(aggParams);
1218  aggFactory->SetFactory("FC Splitting", manager.GetFactory("FC Splitting"));
1219  aggFactory->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1220  aggFactory->SetFactory("DofsPerNode", manager.GetFactory("Graph"));
1221  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
1222 
1223  if (drop_algo.find("block diagonal") != std::string::npos || drop_algo == "signed classical") {
1224  if (levelID > 0)
1225  aggFactory->SetFactory("BlockNumber", this->GetFactoryManager(levelID - 1)->GetFactory("BlockNumber"));
1226  else
1227  aggFactory->SetFactory("BlockNumber", manager.GetFactory("BlockNumber"));
1228  }
1229 
1230  // Now we short-circuit, because we neither need nor want TentativePFactory here
1231  manager.SetFactory("Ptent", aggFactory);
1232  manager.SetFactory("P Graph", aggFactory);
1233 
1234  if (reuseType == "tP" && levelID) {
1235  // keeps.push_back(keep_pair("Nullspace", Ptent.get()));
1236  keeps.push_back(keep_pair("Ptent", aggFactory.get()));
1237  }
1238  return;
1239  } else if (aggType == "notay") {
1240  aggFactory = rcp(new NotayAggregationFactory());
1241  ParameterList aggParams;
1242  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: pairwise: size", int, aggParams);
1243  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: pairwise: tie threshold", double, aggParams);
1244  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: Dirichlet threshold", double, aggParams);
1245  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: ordering", std::string, aggParams);
1246  aggFactory->SetParameterList(aggParams);
1247  aggFactory->SetFactory("DofsPerNode", manager.GetFactory("Graph"));
1248  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
1249  }
1250 #ifdef HAVE_MUELU_MATLAB
1251  else if (aggType == "matlab") {
1252  ParameterList aggParams = paramList.sublist("aggregation: params");
1253  aggFactory = rcp(new SingleLevelMatlabFactory());
1254  aggFactory->SetParameterList(aggParams);
1255  }
1256 #endif
1257 
1258  manager.SetFactory("Aggregates", aggFactory);
1259 
1260  // Coarse map
1261  RCP<Factory> coarseMap = rcp(new CoarseMapFactory());
1262  coarseMap->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1263  manager.SetFactory("CoarseMap", coarseMap);
1264 
1265  // Tentative P
1267  ParameterList ptentParams;
1268  if (paramList.isSublist("matrixmatrix: kernel params"))
1269  ptentParams.sublist("matrixmatrix: kernel params", false) = paramList.sublist("matrixmatrix: kernel params");
1270  if (defaultList.isSublist("matrixmatrix: kernel params"))
1271  ptentParams.sublist("matrixmatrix: kernel params", false) = defaultList.sublist("matrixmatrix: kernel params");
1272  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: calculate qr", bool, ptentParams);
1273  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: build coarse coordinates", bool, ptentParams);
1274  Ptent->SetParameterList(ptentParams);
1275  Ptent->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1276  Ptent->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1277  manager.SetFactory("Ptent", Ptent);
1278 
1279  if (reuseType == "tP" && levelID) {
1280  keeps.push_back(keep_pair("Nullspace", Ptent.get()));
1281  keeps.push_back(keep_pair("P", Ptent.get()));
1282  }
1283 }
1284 
1285 // =====================================================================================================
1286 // ============================================ RAP ====================================================
1287 // =====================================================================================================
1288 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1290  UpdateFactoryManager_RAP(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager,
1291  int levelID, std::vector<keep_pair>& keeps) const {
1292  if (paramList.isParameter("A") && !paramList.get<RCP<Matrix> >("A").is_null()) {
1293  // We have user matrix A
1294  manager.SetFactory("A", NoFactory::getRCP());
1295  return;
1296  }
1297 
1298  ParameterList RAPparams;
1299 
1300  RCP<RAPFactory> RAP;
1301  RCP<RAPShiftFactory> RAPs;
1302  // Allow for Galerkin or shifted RAP
1303  // FIXME: Should this not be some form of MUELU_SET_VAR_2LIST?
1304  std::string alg = paramList.get("rap: algorithm", "galerkin");
1305  if (alg == "shift" || alg == "non-galerkin") {
1306  RAPs = rcp(new RAPShiftFactory());
1307  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift", double, RAPparams);
1308  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift diagonal M", bool, RAPparams);
1309  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift low storage", bool, RAPparams);
1310  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: shift array", Teuchos::Array<double>, RAPparams);
1311  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: cfl array", Teuchos::Array<double>, RAPparams);
1312 
1313  } else {
1314  RAP = rcp(new RAPFactory());
1315  }
1316 
1317  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: relative diagonal floor", Teuchos::Array<double>, RAPparams);
1318 
1319  if (paramList.isSublist("matrixmatrix: kernel params"))
1320  RAPparams.sublist("matrixmatrix: kernel params", false) = paramList.sublist("matrixmatrix: kernel params");
1321  if (defaultList.isSublist("matrixmatrix: kernel params"))
1322  RAPparams.sublist("matrixmatrix: kernel params", false) = defaultList.sublist("matrixmatrix: kernel params");
1323  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "transpose: use implicit", bool, RAPparams);
1324  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: fix zero diagonals", bool, RAPparams);
1325  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: fix zero diagonals threshold", double, RAPparams);
1326  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: fix zero diagonals replacement", Scalar, RAPparams);
1327 
1328  // if "rap: triple product" has not been set and algorithm is "unsmoothed" switch triple product on
1329  if (!paramList.isParameter("rap: triple product") &&
1330  paramList.isType<std::string>("multigrid algorithm") &&
1331  paramList.get<std::string>("multigrid algorithm") == "unsmoothed")
1332  paramList.set("rap: triple product", true);
1333  else
1334  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "rap: triple product", bool, RAPparams);
1335 
1336  try {
1337  if (paramList.isParameter("aggregation: allow empty prolongator columns")) {
1338  RAPparams.set("CheckMainDiagonal", paramList.get<bool>("aggregation: allow empty prolongator columns"));
1339  RAPparams.set("RepairMainDiagonal", paramList.get<bool>("aggregation: allow empty prolongator columns"));
1340  } else if (defaultList.isParameter("aggregation: allow empty prolongator columns")) {
1341  RAPparams.set("CheckMainDiagonal", defaultList.get<bool>("aggregation: allow empty prolongator columns"));
1342  RAPparams.set("RepairMainDiagonal", defaultList.get<bool>("aggregation: allow empty prolongator columns"));
1343  }
1344 
1347  "Error: parameter \"aggregation: allow empty prolongator columns\" must be of type " << Teuchos::TypeNameTraits<bool>::name());
1348  }
1349 
1350  if (!RAP.is_null()) {
1351  RAP->SetParameterList(RAPparams);
1352  RAP->SetFactory("P", manager.GetFactory("P"));
1353  } else {
1354  RAPs->SetParameterList(RAPparams);
1355  RAPs->SetFactory("P", manager.GetFactory("P"));
1356  }
1357 
1358  if (!this->implicitTranspose_) {
1359  if (!RAP.is_null())
1360  RAP->SetFactory("R", manager.GetFactory("R"));
1361  else
1362  RAPs->SetFactory("R", manager.GetFactory("R"));
1363  }
1364 
1365  // Matrix analysis
1366  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "matrix: compute analysis", bool, true)) {
1367  RCP<Factory> matrixAnalysisFact = rcp(new MatrixAnalysisFactory());
1368 
1369  if (!RAP.is_null())
1370  RAP->AddTransferFactory(matrixAnalysisFact);
1371  else
1372  RAPs->AddTransferFactory(matrixAnalysisFact);
1373  }
1374 
1375  // Aggregate qualities
1376  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: compute aggregate qualities", bool, true)) {
1377  RCP<Factory> aggQualityFact = rcp(new AggregateQualityEstimateFactory());
1378  ParameterList aggQualityParams;
1379  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: good aggregate threshold", double, aggQualityParams);
1380  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: file output", bool, aggQualityParams);
1381  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: file base", std::string, aggQualityParams);
1382  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: check symmetry", bool, aggQualityParams);
1383  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: algorithm", std::string, aggQualityParams);
1384  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: zero threshold", double, aggQualityParams);
1385  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: percentiles", Teuchos::Array<double>, aggQualityParams);
1386  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregate qualities: mode", std::string, aggQualityParams);
1387  aggQualityFact->SetParameterList(aggQualityParams);
1388  aggQualityFact->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1389  aggQualityFact->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1390  manager.SetFactory("AggregateQualities", aggQualityFact);
1391 
1392  if (!RAP.is_null())
1393  RAP->AddTransferFactory(aggQualityFact);
1394  else
1395  RAPs->AddTransferFactory(aggQualityFact);
1396  }
1397 
1398  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: export visualization data", bool, true)) {
1400  ParameterList aggExportParams;
1401  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output filename", std::string, aggExportParams);
1402  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: agg style", std::string, aggExportParams);
1403  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: iter", int, aggExportParams);
1404  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: time step", int, aggExportParams);
1405  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: fine graph edges", bool, aggExportParams);
1406  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: coarse graph edges", bool, aggExportParams);
1407  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: build colormap", bool, aggExportParams);
1408  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: aggregate qualities", bool, aggExportParams);
1409  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: material", bool, aggExportParams);
1410  aggExport->SetParameterList(aggExportParams);
1411  aggExport->SetFactory("AggregateQualities", manager.GetFactory("AggregateQualities"));
1412  aggExport->SetFactory("DofsPerNode", manager.GetFactory("DofsPerNode"));
1413  aggExport->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1414  aggExport->SetFactory("Graph", manager.GetFactory("Graph"));
1415 
1416  if (!RAP.is_null())
1417  RAP->AddTransferFactory(aggExport);
1418  else
1419  RAPs->AddTransferFactory(aggExport);
1420  }
1421  if (!RAP.is_null())
1422  manager.SetFactory("A", RAP);
1423  else
1424  manager.SetFactory("A", RAPs);
1425 
1426  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
1427  MUELU_SET_VAR_2LIST(paramList, defaultList, "sa: use filtered matrix", bool, useFiltering);
1428  bool filteringChangesMatrix = useFiltering && !MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol", double, 0);
1429 
1430  if (reuseType == "RP" || (reuseType == "tP" && !filteringChangesMatrix)) {
1431  if (!RAP.is_null()) {
1432  keeps.push_back(keep_pair("AP reuse data", RAP.get()));
1433  keeps.push_back(keep_pair("RAP reuse data", RAP.get()));
1434 
1435  } else {
1436  keeps.push_back(keep_pair("AP reuse data", RAPs.get()));
1437  keeps.push_back(keep_pair("RAP reuse data", RAPs.get()));
1438  }
1439  }
1440 }
1441 
1442 // =====================================================================================================
1443 // ======================================= Coordinates =================================================
1444 // =====================================================================================================
1445 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1448  FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
1449  bool have_userCO = false;
1450  if (paramList.isParameter("Coordinates") && !paramList.get<RCP<MultiVector> >("Coordinates").is_null())
1451  have_userCO = true;
1452 
1453  if (useCoordinates_) {
1454  if (have_userCO) {
1455  manager.SetFactory("Coordinates", NoFactory::getRCP());
1456 
1457  } else {
1459  coords->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1460  coords->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1461  manager.SetFactory("Coordinates", coords);
1462 
1463  auto RAP = rcp_const_cast<RAPFactory>(rcp_dynamic_cast<const RAPFactory>(manager.GetFactory("A")));
1464  if (!RAP.is_null()) {
1465  RAP->AddTransferFactory(manager.GetFactory("Coordinates"));
1466  } else {
1467  auto RAPs = rcp_const_cast<RAPShiftFactory>(rcp_dynamic_cast<const RAPShiftFactory>(manager.GetFactory("A")));
1468  RAPs->AddTransferFactory(manager.GetFactory("Coordinates"));
1469  }
1470  }
1471  }
1472 }
1473 
1474 // ======================================================================================================
1475 // ======================================== Material ==================================================
1476 // =====================================================================================================
1477 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1479  UpdateFactoryManager_Material(ParameterList& paramList, const ParameterList& /* defaultList */,
1480  FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
1481  bool have_userMaterial = false;
1482  if (paramList.isParameter("Material") && !paramList.get<RCP<MultiVector> >("Material").is_null())
1483  have_userMaterial = true;
1484 
1485  if (useMaterial_) {
1486  if (have_userMaterial) {
1487  manager.SetFactory("Material", NoFactory::getRCP());
1488  } else {
1489  RCP<Factory> materialTransfer = rcp(new MultiVectorTransferFactory());
1490  ParameterList materialTransferParameters;
1491  materialTransferParameters.set("Vector name", "Material");
1492  materialTransferParameters.set("Transfer name", "Aggregates");
1493  materialTransferParameters.set("Normalize", true);
1494  materialTransfer->SetParameterList(materialTransferParameters);
1495  materialTransfer->SetFactory("Transfer factory", manager.GetFactory("Aggregates"));
1496  materialTransfer->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1497  manager.SetFactory("Material", materialTransfer);
1498 
1499  auto RAP = rcp_const_cast<RAPFactory>(rcp_dynamic_cast<const RAPFactory>(manager.GetFactory("A")));
1500  if (!RAP.is_null()) {
1501  RAP->AddTransferFactory(manager.GetFactory("Material"));
1502  } else {
1503  auto RAPs = rcp_const_cast<RAPShiftFactory>(rcp_dynamic_cast<const RAPShiftFactory>(manager.GetFactory("A")));
1504  RAPs->AddTransferFactory(manager.GetFactory("Material"));
1505  }
1506  }
1507  }
1508 }
1509 
1510 // =====================================================================================================
1511 // ================================= LocalOrdinalTransfer =============================================
1512 // =====================================================================================================
1513 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1515  UpdateFactoryManager_LocalOrdinalTransfer(const std::string& VarName, const std::string& multigridAlgo, ParameterList& paramList, const ParameterList& /* defaultList */,
1516  FactoryManager& manager, int levelID, std::vector<keep_pair>& /* keeps */) const {
1517  // NOTE: You would think this would be levelID > 0, but you'd be wrong, since the FactoryManager is basically
1518  // offset by a level from the things which actually do the work.
1519  if (useBlockNumber_ && (levelID > 0)) {
1520  auto RAP = rcp_const_cast<RAPFactory>(rcp_dynamic_cast<const RAPFactory>(manager.GetFactory("A")));
1521  auto RAPs = rcp_const_cast<RAPShiftFactory>(rcp_dynamic_cast<const RAPShiftFactory>(manager.GetFactory("A")));
1522  if (!RAP.is_null() || !RAPs.is_null()) {
1523  RCP<Factory> fact = rcp(new LocalOrdinalTransferFactory(VarName, multigridAlgo));
1524  if (multigridAlgo == "classical")
1525  fact->SetFactory("P Graph", manager.GetFactory("P Graph"));
1526  else
1527  fact->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1528  fact->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1529 
1530  fact->SetFactory(VarName, this->GetFactoryManager(levelID - 1)->GetFactory(VarName));
1531 
1532  manager.SetFactory(VarName, fact);
1533 
1534  if (!RAP.is_null())
1535  RAP->AddTransferFactory(manager.GetFactory(VarName));
1536  else
1537  RAPs->AddTransferFactory(manager.GetFactory(VarName));
1538  }
1539  }
1540 }
1541 
1542 // ======================================================================================================
1543 // ====================================== BlockNumber =================================================
1544 // =====================================================================================================
1545 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1548  FactoryManager& manager, int levelID, std::vector<keep_pair>& keeps) const {
1549  if (useBlockNumber_) {
1550  ParameterList myParams;
1552  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: block diagonal: interleaved blocksize", int, myParams);
1553  fact->SetParameterList(myParams);
1554  manager.SetFactory("BlockNumber", fact);
1555  }
1556 }
1557 
1558 // =====================================================================================================
1559 // =========================================== Restriction =============================================
1560 // =====================================================================================================
1561 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1564  int levelID, std::vector<keep_pair>& /* keeps */) const {
1565  MUELU_SET_VAR_2LIST(paramList, defaultList, "multigrid algorithm", std::string, multigridAlgo);
1566  bool have_userR = false;
1567  if (paramList.isParameter("R") && !paramList.get<RCP<Matrix> >("R").is_null())
1568  have_userR = true;
1569 
1570  // === Restriction ===
1571  RCP<Factory> R;
1572  if (!this->implicitTranspose_) {
1573  MUELU_SET_VAR_2LIST(paramList, defaultList, "problem: symmetric", bool, isSymmetric);
1574 
1575  if (isSymmetric == false && (multigridAlgo == "unsmoothed" || multigridAlgo == "emin")) {
1576  this->GetOStream(Warnings0) << "Switching \"problem: symmetric\" parameter to symmetric as multigrid algorithm. " << multigridAlgo << " is primarily supposed to be used for symmetric problems.\n\n"
1577  << "Please note: if you are using \"unsmoothed\" transfer operators the \"problem: symmetric\" parameter "
1578  << "has no real mathematical meaning, i.e. you can use it for non-symmetric\n"
1579  << "problems, too. With \"problem: symmetric\"=\"symmetric\" you can use implicit transpose for building "
1580  << "the restriction operators which may drastically reduce the amount of consumed memory." << std::endl;
1581  isSymmetric = true;
1582  }
1583  TEUCHOS_TEST_FOR_EXCEPTION(multigridAlgo == "pg" && isSymmetric == true, Exceptions::RuntimeError,
1584  "Petrov-Galerkin smoothed transfer operators are only allowed for non-symmetric problems: Set \"problem: symmetric\" to false!\n"
1585  "While PG smoothed transfer operators generally would also work for symmetric problems this is an unusual use case. "
1586  "You can use the factory-based xml interface though if you need PG-AMG for symmetric problems.");
1587 
1588  if (have_userR) {
1589  manager.SetFactory("R", NoFactory::getRCP());
1590  } else {
1591  if (isSymmetric)
1592  R = rcp(new TransPFactory());
1593  else
1594  R = rcp(new GenericRFactory());
1595 
1596  R->SetFactory("P", manager.GetFactory("P"));
1597  manager.SetFactory("R", R);
1598  }
1599 
1600  } else {
1601  manager.SetFactory("R", Teuchos::null);
1602  }
1603 
1604  // === Restriction: Nullspace Scaling ===
1605  if (paramList.isParameter("restriction: scale nullspace") && paramList.get<bool>("restriction: scale nullspace")) {
1606  RCP<TentativePFactory> tentPFactory = rcp(new TentativePFactory());
1607  Teuchos::ParameterList tentPlist;
1608  tentPlist.set("Nullspace name", "Scaled Nullspace");
1609  tentPFactory->SetParameterList(tentPlist);
1610  tentPFactory->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1611  tentPFactory->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1612 
1613  if (R.is_null()) R = rcp(new TransPFactory());
1614  R->SetFactory("P", tentPFactory);
1615  }
1616 }
1617 
1618 // =====================================================================================================
1619 // ========================================= Repartition ===============================================
1620 // =====================================================================================================
1621 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1624  int levelID, std::vector<keep_pair>& keeps, RCP<Factory>& nullSpaceFactory) const {
1625  // === Repartitioning ===
1626  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
1627  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: enable", bool, enableRepart);
1628  if (enableRepart) {
1629 #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
1630  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: use subcommunicators in place", bool, enableInPlace);
1631  // Short summary of the issue: RebalanceTransferFactory shares ownership
1632  // of "P" with SaPFactory, and therefore, changes the stored version.
1633  // That means that if SaPFactory generated P, and stored it on the level,
1634  // then after rebalancing the value in that storage changed. It goes
1635  // against the concept of factories (I think), that every factory is
1636  // responsible for its own objects, and they are immutable outside.
1637  //
1638  // In reuse, this is what happens: as we reuse Importer across setups,
1639  // the order of factories changes, and coupled with shared ownership
1640  // leads to problems.
1641  // *First setup*
1642  // SaP builds P [and stores it]
1643  // TransP builds R [and stores it]
1644  // RAP builds A [and stores it]
1645  // RebalanceTransfer rebalances P [and changes the P stored by SaP] (*)
1646  // RebalanceTransfer rebalances R
1647  // RebalanceAc rebalances A
1648  // *Second setup* ("RP" reuse)
1649  // RebalanceTransfer rebalances P [which is incorrect due to (*)]
1650  // RebalanceTransfer rebalances R
1651  // RAP builds A [which is incorrect due to (*)]
1652  // RebalanceAc rebalances A [which throws due to map inconsistency]
1653  // ...
1654  // *Second setup* ("tP" reuse)
1655  // SaP builds P [and stores it]
1656  // RebalanceTransfer rebalances P [and changes the P stored by SaP] (**)
1657  // TransP builds R [which is incorrect due to (**)]
1658  // RebalanceTransfer rebalances R
1659  // ...
1660  //
1661  // Couple solutions to this:
1662  // 1. [implemented] Requre "tP" and "PR" reuse to only be used with
1663  // implicit rebalancing.
1664  // 2. Do deep copy of P, and changed domain map and importer there.
1665  // Need to investigate how expensive this is.
1666  TEUCHOS_TEST_FOR_EXCEPTION(this->doPRrebalance_ && (reuseType == "tP" || reuseType == "RP"), Exceptions::InvalidArgument,
1667  "Reuse types \"tP\" and \"PR\" require \"repartition: rebalance P and R\" set to \"false\"");
1668 
1669  // TEUCHOS_TEST_FOR_EXCEPTION(aggType == "brick", Exceptions::InvalidArgument,
1670  // "Aggregation type \"brick\" requires \"repartition: enable\" set to \"false\"");
1671 
1672  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: partitioner", std::string, partName);
1673  TEUCHOS_TEST_FOR_EXCEPTION(partName != "zoltan" && partName != "zoltan2", Exceptions::InvalidArgument,
1674  "Invalid partitioner name: \"" << partName << "\". Valid options: \"zoltan\", \"zoltan2\"");
1675 
1676 #ifndef HAVE_MUELU_ZOLTAN
1677  bool switched = false;
1678  if (partName == "zoltan") {
1679  this->GetOStream(Warnings0) << "Zoltan interface is not available, trying to switch to Zoltan2" << std::endl;
1680  partName = "zoltan2";
1681  switched = true;
1682  }
1683 #else
1684 #ifndef HAVE_MUELU_ZOLTAN2
1685  bool switched = false;
1686 #endif // HAVE_MUELU_ZOLTAN2
1687 #endif // HAVE_MUELU_ZOLTAN
1688 
1689 #ifndef HAVE_MUELU_ZOLTAN2
1690  if (partName == "zoltan2" && !switched) {
1691  this->GetOStream(Warnings0) << "Zoltan2 interface is not available, trying to switch to Zoltan" << std::endl;
1692  partName = "zoltan";
1693  }
1694 #endif // HAVE_MUELU_ZOLTAN2
1695 
1696  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: node repartition level", int, nodeRepartitionLevel);
1697 
1698  // RepartitionHeuristic
1699  auto repartheurFactory = rcp(new RepartitionHeuristicFactory());
1700  ParameterList repartheurParams;
1701  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: node repartition level", int, repartheurParams);
1702  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: start level", int, repartheurParams);
1703  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: min rows per proc", int, repartheurParams);
1704  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: target rows per proc", int, repartheurParams);
1705  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: min rows per thread", int, repartheurParams);
1706  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: target rows per thread", int, repartheurParams);
1707  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: max imbalance", double, repartheurParams);
1708  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: put on single proc", int, repartheurParams);
1709  repartheurFactory->SetParameterList(repartheurParams);
1710  repartheurFactory->SetFactory("A", manager.GetFactory("A"));
1711  manager.SetFactory("number of partitions", repartheurFactory);
1712  manager.SetFactory("repartition: heuristic target rows per process", repartheurFactory);
1713 
1714  // Partitioner
1715  RCP<Factory> partitioner;
1716  if (levelID == nodeRepartitionLevel) {
1717  // partitioner = rcp(new NodePartitionInterface());
1719  ParameterList partParams;
1720  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: node id", int, repartheurParams);
1721  partitioner->SetParameterList(partParams);
1722  partitioner->SetFactory("Node Comm", manager.GetFactory("Node Comm"));
1723  } else if (partName == "zoltan") {
1724 #ifdef HAVE_MUELU_ZOLTAN
1725  partitioner = rcp(new ZoltanInterface());
1726  // NOTE: ZoltanInterface ("zoltan") does not support external parameters through ParameterList
1727 #else
1728  throw Exceptions::RuntimeError("Zoltan interface is not available");
1729 #endif // HAVE_MUELU_ZOLTAN
1730  } else if (partName == "zoltan2") {
1731 #ifdef HAVE_MUELU_ZOLTAN2
1732  partitioner = rcp(new Zoltan2Interface());
1733  ParameterList partParams;
1734  RCP<const ParameterList> partpartParams = rcp(new ParameterList(paramList.sublist("repartition: params", false)));
1735  partParams.set("ParameterList", partpartParams);
1736  partitioner->SetParameterList(partParams);
1737  partitioner->SetFactory("repartition: heuristic target rows per process",
1738  manager.GetFactory("repartition: heuristic target rows per process"));
1739 #else
1740  throw Exceptions::RuntimeError("Zoltan2 interface is not available");
1741 #endif // HAVE_MUELU_ZOLTAN2
1742  }
1743 
1744  partitioner->SetFactory("A", manager.GetFactory("A"));
1745  partitioner->SetFactory("number of partitions", manager.GetFactory("number of partitions"));
1746  if (useCoordinates_)
1747  partitioner->SetFactory("Coordinates", manager.GetFactory("Coordinates"));
1748  manager.SetFactory("Partition", partitioner);
1749 
1750  // Repartitioner
1751  auto repartFactory = rcp(new RepartitionFactory());
1752  ParameterList repartParams;
1753  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: print partition distribution", bool, repartParams);
1754  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap parts", bool, repartParams);
1755  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap num values", int, repartParams);
1756  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: save importer", bool, repartParams);
1757  repartFactory->SetParameterList(repartParams);
1758  repartFactory->SetFactory("A", manager.GetFactory("A"));
1759  repartFactory->SetFactory("number of partitions", manager.GetFactory("number of partitions"));
1760  repartFactory->SetFactory("Partition", manager.GetFactory("Partition"));
1761  manager.SetFactory("Importer", repartFactory);
1762  if (reuseType != "none" && reuseType != "S" && levelID)
1763  keeps.push_back(keep_pair("Importer", manager.GetFactory("Importer").get()));
1764 
1765  if (enableInPlace) {
1766  // Rebalanced A (in place)
1767  // NOTE: This is for when we want to constrain repartitioning to match some other idea of what's going on.
1768  // The major application is the (1,1) hierarchy in the Maxwell1 preconditioner.
1769  auto newA = rcp(new RebalanceAcFactory());
1770  ParameterList rebAcParams;
1771  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, rebAcParams);
1772  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators in place", bool, rebAcParams);
1773  newA->SetParameterList(rebAcParams);
1774  newA->SetFactory("A", manager.GetFactory("A"));
1775  newA->SetFactory("InPlaceMap", manager.GetFactory("InPlaceMap"));
1776  manager.SetFactory("A", newA);
1777  } else {
1778  // Rebalanced A
1779  auto newA = rcp(new RebalanceAcFactory());
1780  ParameterList rebAcParams;
1781  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, rebAcParams);
1782  newA->SetParameterList(rebAcParams);
1783  newA->SetFactory("A", manager.GetFactory("A"));
1784  newA->SetFactory("Importer", manager.GetFactory("Importer"));
1785  manager.SetFactory("A", newA);
1786 
1787  // Rebalanced P
1788  auto newP = rcp(new RebalanceTransferFactory());
1789  ParameterList newPparams;
1790  newPparams.set("type", "Interpolation");
1791  if (changedPRrebalance_)
1792  newPparams.set("repartition: rebalance P and R", this->doPRrebalance_);
1793  if (changedPRViaCopyrebalance_)
1794  newPparams.set("repartition: explicit via new copy rebalance P and R", true);
1795  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, newPparams);
1796  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: send type", std::string, newPparams);
1797  newP->SetParameterList(newPparams);
1798  newP->SetFactory("Importer", manager.GetFactory("Importer"));
1799  newP->SetFactory("P", manager.GetFactory("P"));
1800  manager.SetFactory("P", newP);
1801  if (!paramList.isParameter("semicoarsen: number of levels"))
1802  newP->SetFactory("Nullspace", manager.GetFactory("Ptent"));
1803  else
1804  newP->SetFactory("Nullspace", manager.GetFactory("P")); // TogglePFactory
1805  if (useCoordinates_) {
1806  newP->SetFactory("Coordinates", manager.GetFactory("Coordinates"));
1807  manager.SetFactory("Coordinates", newP);
1808  }
1809  if (useMaterial_) {
1810  newP->SetFactory("Material", manager.GetFactory("Material"));
1811  manager.SetFactory("Material", newP);
1812  }
1813  if (useBlockNumber_ && (levelID > 0)) {
1814  newP->SetFactory("BlockNumber", manager.GetFactory("BlockNumber"));
1815  manager.SetFactory("BlockNumber", newP);
1816  }
1817 
1818  // Rebalanced R
1819  auto newR = rcp(new RebalanceTransferFactory());
1820  ParameterList newRparams;
1821  newRparams.set("type", "Restriction");
1822  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, newRparams);
1823  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: send type", std::string, newRparams);
1824  if (changedPRrebalance_)
1825  newRparams.set("repartition: rebalance P and R", this->doPRrebalance_);
1826  if (changedPRViaCopyrebalance_)
1827  newPparams.set("repartition: explicit via new copy rebalance P and R", true);
1828  if (changedImplicitTranspose_)
1829  newRparams.set("transpose: use implicit", this->implicitTranspose_);
1830  newR->SetParameterList(newRparams);
1831  newR->SetFactory("Importer", manager.GetFactory("Importer"));
1832  if (!this->implicitTranspose_) {
1833  newR->SetFactory("R", manager.GetFactory("R"));
1834  manager.SetFactory("R", newR);
1835  }
1836 
1837  // NOTE: the role of NullspaceFactory is to provide nullspace on the finest
1838  // level if a user does not do that. For all other levels it simply passes
1839  // nullspace from a real factory to whoever needs it. If we don't use
1840  // repartitioning, that factory is "TentativePFactory"; if we do, it is
1841  // "RebalanceTransferFactory". But we still have to have NullspaceFactory as
1842  // the "Nullspace" of the manager
1843  // NOTE: This really needs to be set on the *NullSpaceFactory*, not manager.get("Nullspace").
1844  ParameterList newNullparams;
1845  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "nullspace: calculate rotations", bool, newNullparams);
1846  nullSpaceFactory->SetFactory("Nullspace", newP);
1847  nullSpaceFactory->SetParameterList(newNullparams);
1848  }
1849 #else
1850  paramList.set("repartition: enable", false);
1851 #ifndef HAVE_MPI
1852  this->GetOStream(Warnings0) << "No repartitioning available for a serial run\n";
1853 #else
1854  this->GetOStream(Warnings0) << "Zoltan/Zoltan2 are unavailable for repartitioning\n";
1855 #endif // HAVE_MPI
1856 #endif // defined(HAVE_MPI) && (defined(HAVE_MUELU_ZOLTAN) || defined(HAVE_MUELU_ZOLTAN2))
1857  }
1858 }
1859 
1860 // =====================================================================================================
1861 // ========================================= Low precision transfers ===================================
1862 // =====================================================================================================
1863 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1866  int levelID, std::vector<keep_pair>& keeps) const {
1867  MUELU_SET_VAR_2LIST(paramList, defaultList, "transfers: half precision", bool, enableLowPrecision);
1868 
1869  if (enableLowPrecision) {
1870  // Low precision P
1871  auto newP = rcp(new LowPrecisionFactory());
1872  ParameterList newPparams;
1873  newPparams.set("matrix key", "P");
1874  newP->SetParameterList(newPparams);
1875  newP->SetFactory("P", manager.GetFactory("P"));
1876  manager.SetFactory("P", newP);
1877 
1878  if (!this->implicitTranspose_) {
1879  // Low precision R
1880  auto newR = rcp(new LowPrecisionFactory());
1881  ParameterList newRparams;
1882  newRparams.set("matrix key", "R");
1883  newR->SetParameterList(newRparams);
1884  newR->SetFactory("R", manager.GetFactory("R"));
1885  manager.SetFactory("R", newR);
1886  }
1887  }
1888 }
1889 
1890 // =====================================================================================================
1891 // =========================================== Nullspace ===============================================
1892 // =====================================================================================================
1893 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1896  int /* levelID */, std::vector<keep_pair>& /* keeps */, RCP<Factory>& nullSpaceFactory) const {
1897  // Nullspace
1898  RCP<Factory> nullSpace = rcp(new NullspaceFactory());
1899 
1900  bool have_userNS = false;
1901  if (paramList.isParameter("Nullspace") && !paramList.get<RCP<MultiVector> >("Nullspace").is_null())
1902  have_userNS = true;
1903 
1904  if (!have_userNS) {
1905  ParameterList newNullparams;
1906  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "nullspace: calculate rotations", bool, newNullparams);
1907  nullSpace->SetParameterList(newNullparams);
1908  nullSpace->SetFactory("Nullspace", manager.GetFactory("Ptent"));
1909  manager.SetFactory("Nullspace", nullSpace);
1910  }
1911  nullSpaceFactory = nullSpace;
1912 
1913  if (paramList.isParameter("restriction: scale nullspace") && paramList.get<bool>("restriction: scale nullspace")) {
1914  RCP<ScaledNullspaceFactory> scaledNSfactory = rcp(new ScaledNullspaceFactory());
1915  scaledNSfactory->SetFactory("Nullspace", nullSpaceFactory);
1916  manager.SetFactory("Scaled Nullspace", scaledNSfactory);
1917  }
1918 }
1919 
1920 // =====================================================================================================
1921 // ================================= Algorithm: SemiCoarsening =========================================
1922 // =====================================================================================================
1923 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1926  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
1927  // === Semi-coarsening ===
1928  RCP<Factory> semicoarsenFactory = Teuchos::null;
1929  if (paramList.isParameter("semicoarsen: number of levels") &&
1930  paramList.get<int>("semicoarsen: number of levels") > 0) {
1931  ParameterList togglePParams;
1932  ParameterList semicoarsenPParams;
1933  ParameterList linedetectionParams;
1934  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: number of levels", int, togglePParams);
1935  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: coarsen rate", int, semicoarsenPParams);
1936  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: piecewise constant", bool, semicoarsenPParams);
1937  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: piecewise linear", bool, semicoarsenPParams);
1938  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: calculate nonsym restriction", bool, semicoarsenPParams);
1939  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "linedetection: orientation", std::string, linedetectionParams);
1940  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "linedetection: num layers", int, linedetectionParams);
1941 
1943  RCP<LineDetectionFactory> linedetectionFactory = rcp(new LineDetectionFactory());
1944  RCP<TogglePFactory> togglePFactory = rcp(new TogglePFactory());
1945 
1946  linedetectionFactory->SetParameterList(linedetectionParams);
1947  semicoarsenFactory->SetParameterList(semicoarsenPParams);
1948  togglePFactory->SetParameterList(togglePParams);
1949 
1950  togglePFactory->AddCoarseNullspaceFactory(semicoarsenFactory);
1951  togglePFactory->AddProlongatorFactory(semicoarsenFactory);
1952  togglePFactory->AddPtentFactory(semicoarsenFactory);
1953  togglePFactory->AddCoarseNullspaceFactory(manager.GetFactory("Ptent"));
1954  togglePFactory->AddProlongatorFactory(manager.GetFactory("P"));
1955  togglePFactory->AddPtentFactory(manager.GetFactory("Ptent"));
1956 
1957  manager.SetFactory("CoarseNumZLayers", linedetectionFactory);
1958  manager.SetFactory("LineDetection_Layers", linedetectionFactory);
1959  manager.SetFactory("LineDetection_VertLineIds", linedetectionFactory);
1960 
1961  manager.SetFactory("P", togglePFactory);
1962  manager.SetFactory("Ptent", togglePFactory);
1963  manager.SetFactory("Nullspace", togglePFactory);
1964  }
1965 
1966  if (paramList.isParameter("semicoarsen: number of levels")) {
1967  auto tf = rcp(new ToggleCoordinatesTransferFactory());
1968  tf->SetFactory("Chosen P", manager.GetFactory("P"));
1969  tf->AddCoordTransferFactory(semicoarsenFactory);
1970 
1972  coords->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
1973  coords->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
1974  tf->AddCoordTransferFactory(coords);
1975  manager.SetFactory("Coordinates", tf);
1976  }
1977 }
1978 
1979 // =====================================================================================================
1980 // ================================== Algorithm: P-Coarsening ==========================================
1981 // =====================================================================================================
1982 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1985  int levelID, std::vector<keep_pair>& keeps) const {
1986 #ifdef HAVE_MUELU_INTREPID2
1987  // This only makes sense to invoke from the default list.
1988  if (defaultList.isParameter("pcoarsen: schedule") && defaultList.isParameter("pcoarsen: element")) {
1989  // P-Coarsening by schedule (new interface)
1990  // NOTE: levelID represents the *coarse* level in this case
1991  auto pcoarsen_schedule = Teuchos::getArrayFromStringParameter<int>(defaultList, "pcoarsen: schedule");
1992  auto pcoarsen_element = defaultList.get<std::string>("pcoarsen: element");
1993 
1994  if (levelID >= (int)pcoarsen_schedule.size()) {
1995  // Past the p-coarsening levels, we do Smoothed Aggregation
1996  // NOTE: We should probably consider allowing other options past p-coarsening
1997  UpdateFactoryManager_SA(paramList, defaultList, manager, levelID, keeps);
1998 
1999  } else {
2000  // P-Coarsening
2001  ParameterList Pparams;
2002  auto P = rcp(new IntrepidPCoarsenFactory());
2003  std::string lo = pcoarsen_element + std::to_string(pcoarsen_schedule[levelID]);
2004  std::string hi = (levelID ? pcoarsen_element + std::to_string(pcoarsen_schedule[levelID - 1]) : lo);
2005  Pparams.set("pcoarsen: hi basis", hi);
2006  Pparams.set("pcoarsen: lo basis", lo);
2007  P->SetParameterList(Pparams);
2008  manager.SetFactory("P", P);
2009 
2010  // Add special nullspace handling
2011  rcp_dynamic_cast<Factory>(manager.GetFactoryNonConst("Nullspace"))->SetFactory("Nullspace", manager.GetFactory("P"));
2012  }
2013 
2014  } else {
2015  // P-Coarsening by manual specification (old interface)
2016  ParameterList Pparams;
2017  auto P = rcp(new IntrepidPCoarsenFactory());
2018  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "pcoarsen: hi basis", std::string, Pparams);
2019  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "pcoarsen: lo basis", std::string, Pparams);
2020  P->SetParameterList(Pparams);
2021  manager.SetFactory("P", P);
2022 
2023  // Add special nullspace handling
2024  rcp_dynamic_cast<Factory>(manager.GetFactoryNonConst("Nullspace"))->SetFactory("Nullspace", manager.GetFactory("P"));
2025  }
2026 
2027 #endif
2028 }
2029 
2030 // =====================================================================================================
2031 // ============================== Algorithm: Smoothed Aggregation ======================================
2032 // =====================================================================================================
2033 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2035  UpdateFactoryManager_SA(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& keeps) const {
2036  // Smoothed aggregation
2037  RCP<Factory> P = rcp(new SaPFactory());
2038  ParameterList Pparams;
2039  if (paramList.isSublist("matrixmatrix: kernel params"))
2040  Pparams.sublist("matrixmatrix: kernel params", false) = paramList.sublist("matrixmatrix: kernel params");
2041  if (defaultList.isSublist("matrixmatrix: kernel params"))
2042  Pparams.sublist("matrixmatrix: kernel params", false) = defaultList.sublist("matrixmatrix: kernel params");
2043  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: damping factor", double, Pparams);
2044  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: calculate eigenvalue estimate", bool, Pparams);
2045  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: max eigenvalue", double, Pparams);
2046  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: eigenvalue estimate num iterations", int, Pparams);
2047  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: use rowsumabs diagonal scaling", bool, Pparams);
2048  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs diagonal replacement tolerance", double, Pparams);
2049  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs diagonal replacement value", double, Pparams);
2050  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs use automatic diagonal tolerance", bool, Pparams);
2051  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: enforce constraints", bool, Pparams);
2052  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: eigen-analysis type", std::string, Pparams);
2053  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: calculate qr", bool, Pparams);
2054 
2055  P->SetParameterList(Pparams);
2056 
2057  // Filtering
2058  MUELU_SET_VAR_2LIST(paramList, defaultList, "sa: use filtered matrix", bool, useFiltering);
2059  if (useFiltering) {
2060  // NOTE: Here, non-Kokkos and Kokkos versions diverge in the way the
2061  // dependency tree is setup. The Kokkos version has merged the the
2062  // FilteredAFactory into the CoalesceDropFactory.
2063  if (!useKokkos_) {
2064  RCP<Factory> filterFactory = rcp(new FilteredAFactory());
2065 
2066  ParameterList fParams;
2067  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use lumping", bool, fParams);
2068  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse graph", bool, fParams);
2069  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse eigenvalue", bool, fParams);
2070  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use root stencil", bool, fParams);
2071  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: Dirichlet threshold", double, fParams);
2072  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use spread lumping", bool, fParams);
2073  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: lumping choice", std::string, fParams);
2074  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom growth factor", double, fParams);
2075  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom cap", double, fParams);
2076  filterFactory->SetParameterList(fParams);
2077  filterFactory->SetFactory("Graph", manager.GetFactory("Graph"));
2078  filterFactory->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
2079  filterFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
2080  // I'm not sure why we need this line. See comments for DofsPerNode for UncoupledAggregation above
2081  filterFactory->SetFactory("Filtering", manager.GetFactory("Graph"));
2082 
2083  P->SetFactory("A", filterFactory);
2084 
2085  } else {
2086  P->SetFactory("A", manager.GetFactory("Graph"));
2087  }
2088  }
2089 
2090  P->SetFactory("P", manager.GetFactory("Ptent"));
2091  manager.SetFactory("P", P);
2092 
2093  bool filteringChangesMatrix = useFiltering && !MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol", double, 0);
2094  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
2095  if (reuseType == "tP" && !filteringChangesMatrix)
2096  keeps.push_back(keep_pair("AP reuse data", P.get()));
2097 }
2098 
2099 // =====================================================================================================
2100 // =============================== Algorithm: Energy Minimization ======================================
2101 // =====================================================================================================
2102 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2105  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
2106  MUELU_SET_VAR_2LIST(paramList, defaultList, "emin: pattern", std::string, patternType);
2107  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
2108  TEUCHOS_TEST_FOR_EXCEPTION(patternType != "AkPtent", Exceptions::InvalidArgument,
2109  "Invalid pattern name: \"" << patternType << "\". Valid options: \"AkPtent\"");
2110  // Pattern
2111  auto patternFactory = rcp(new PatternFactory());
2112  ParameterList patternParams;
2113  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: pattern order", int, patternParams);
2114  patternFactory->SetParameterList(patternParams);
2115  patternFactory->SetFactory("P", manager.GetFactory("Ptent"));
2116 
2117  // Filtering
2118  MUELU_SET_VAR_2LIST(paramList, defaultList, "emin: use filtered matrix", bool, useFiltering);
2119  if (useFiltering) {
2120  // NOTE: Here, non-Kokkos and Kokkos versions diverge in the way the
2121  // dependency tree is setup. The Kokkos version has merged the the
2122  // FilteredAFactory into the CoalesceDropFactory.
2123  if (!useKokkos_) {
2124  RCP<Factory> filterFactory = rcp(new FilteredAFactory());
2125 
2126  ParameterList fParams;
2127  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use lumping", bool, fParams);
2128  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse graph", bool, fParams);
2129  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse eigenvalue", bool, fParams);
2130  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use root stencil", bool, fParams);
2131  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: Dirichlet threshold", double, fParams);
2132  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use spread lumping", bool, fParams);
2133  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: lumping choice", std::string, fParams);
2134  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom growth factor", double, fParams);
2135  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: spread lumping diag dom cap", double, fParams);
2136  filterFactory->SetParameterList(fParams);
2137  filterFactory->SetFactory("Graph", manager.GetFactory("Graph"));
2138  filterFactory->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
2139  filterFactory->SetFactory("UnAmalgamationInfo", manager.GetFactory("UnAmalgamationInfo"));
2140  // I'm not sure why we need this line. See comments for DofsPerNode for UncoupledAggregation above
2141  filterFactory->SetFactory("Filtering", manager.GetFactory("Graph"));
2142 
2143  patternFactory->SetFactory("A", filterFactory);
2144 
2145  } else {
2146  patternFactory->SetFactory("A", manager.GetFactory("Graph"));
2147  }
2148  }
2149 
2150  manager.SetFactory("Ppattern", patternFactory);
2151 
2152  // Constraint
2153  auto constraintFactory = rcp(new ConstraintFactory());
2154  constraintFactory->SetFactory("Ppattern", manager.GetFactory("Ppattern"));
2155  constraintFactory->SetFactory("CoarseNullspace", manager.GetFactory("Ptent"));
2156  manager.SetFactory("Constraint", constraintFactory);
2157 
2158  // Energy minimization
2159  ParameterList Pparams;
2160  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: num iterations", int, Pparams);
2161  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: iterative method", std::string, Pparams);
2162  if (reuseType == "emin") {
2163  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: num reuse iterations", int, Pparams);
2164  Pparams.set("Keep P0", true);
2165  Pparams.set("Keep Constraint0", true);
2166  }
2167 
2168  // Emin Factory
2169  auto P = rcp(new EminPFactory());
2170  P->SetParameterList(Pparams);
2171  P->SetFactory("P", manager.GetFactory("Ptent"));
2172  P->SetFactory("Constraint", manager.GetFactory("Constraint"));
2173  manager.SetFactory("P", P);
2174 }
2175 
2176 // =====================================================================================================
2177 // ================================= Algorithm: Petrov-Galerkin ========================================
2178 // =====================================================================================================
2179 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2181  UpdateFactoryManager_PG(ParameterList& /* paramList */, const ParameterList& /* defaultList */, FactoryManager& manager,
2182  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
2183  TEUCHOS_TEST_FOR_EXCEPTION(this->implicitTranspose_, Exceptions::RuntimeError,
2184  "Implicit transpose not supported with Petrov-Galerkin smoothed transfer operators: Set \"transpose: use implicit\" to false!\n"
2185  "Petrov-Galerkin transfer operator smoothing for non-symmetric problems requires a separate handling of the restriction operator which "
2186  "does not allow the usage of implicit transpose easily.");
2187 
2188  // Petrov-Galerkin
2189  auto P = rcp(new PgPFactory());
2190  P->SetFactory("P", manager.GetFactory("Ptent"));
2191  manager.SetFactory("P", P);
2192 }
2193 
2194 // =====================================================================================================
2195 // ================================= Algorithm: Replicate ========================================
2196 // =====================================================================================================
2197 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2199  UpdateFactoryManager_Replicate(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& keeps) const {
2201 
2202  ParameterList Pparams;
2203  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "replicate: npdes", int, Pparams);
2204 
2205  P->SetParameterList(Pparams);
2206  manager.SetFactory("P", P);
2207 }
2208 
2209 // =====================================================================================================
2210 // ====================================== Algorithm: Combine ============================================
2211 // =====================================================================================================
2212 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2214  UpdateFactoryManager_Combine(ParameterList& paramList, const ParameterList& defaultList, FactoryManager& manager, int /* levelID */, std::vector<keep_pair>& keeps) const {
2216 
2217  ParameterList Pparams;
2218  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "combine: numBlks", int, Pparams);
2219 
2220  P->SetParameterList(Pparams);
2221  manager.SetFactory("P", P);
2222 }
2223 
2224 // =====================================================================================================
2225 // ====================================== Algorithm: Matlab ============================================
2226 // =====================================================================================================
2227 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2229  UpdateFactoryManager_Matlab(ParameterList& paramList, const ParameterList& /* defaultList */, FactoryManager& manager,
2230  int /* levelID */, std::vector<keep_pair>& /* keeps */) const {
2231 #ifdef HAVE_MUELU_MATLAB
2232  ParameterList Pparams = paramList.sublist("transfer: params");
2233  auto P = rcp(new TwoLevelMatlabFactory());
2234  P->SetParameterList(Pparams);
2235  P->SetFactory("P", manager.GetFactory("Ptent"));
2236  manager.SetFactory("P", P);
2237 #else
2238  (void)paramList;
2239  (void)manager;
2240 #endif
2241 }
2242 
2243 #undef MUELU_SET_VAR_2LIST
2244 #undef MUELU_TEST_AND_SET_VAR
2245 #undef MUELU_TEST_AND_SET_PARAM_2LIST
2246 #undef MUELU_TEST_PARAM_2LIST
2247 #undef MUELU_KOKKOS_FACTORY
2248 
2249 size_t LevenshteinDistance(const char* s, size_t len_s, const char* t, size_t len_t);
2250 
2251 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2253  ParameterList paramList = constParamList;
2254  const ParameterList& validList = *MasterList::List();
2255  // Validate up to maxLevels level specific parameter sublists
2256  const int maxLevels = 100;
2257 
2258  // Extract level specific list
2259  std::vector<ParameterList> paramLists;
2260  for (int levelID = 0; levelID < maxLevels; levelID++) {
2261  std::string sublistName = "level " + toString(levelID);
2262  if (paramList.isSublist(sublistName)) {
2263  paramLists.push_back(paramList.sublist(sublistName));
2264  // paramLists.back().setName(sublistName);
2265  paramList.remove(sublistName);
2266  }
2267  }
2268  paramLists.push_back(paramList);
2269  // paramLists.back().setName("main");
2270 #ifdef HAVE_MUELU_MATLAB
2271  // If Muemex is supported, hide custom level variables from validator by removing them from paramList's sublists
2272  for (size_t i = 0; i < paramLists.size(); i++) {
2273  std::vector<std::string> customVars; // list of names (keys) to be removed from list
2274 
2275  for (Teuchos::ParameterList::ConstIterator it = paramLists[i].begin(); it != paramLists[i].end(); it++) {
2276  std::string paramName = paramLists[i].name(it);
2277 
2278  if (IsParamMuemexVariable(paramName))
2279  customVars.push_back(paramName);
2280  }
2281 
2282  // Remove the keys
2283  for (size_t j = 0; j < customVars.size(); j++)
2284  paramLists[i].remove(customVars[j], false);
2285  }
2286 #endif
2287 
2288  const int maxDepth = 0;
2289  for (size_t i = 0; i < paramLists.size(); i++) {
2290  // validate every sublist
2291  try {
2292  paramLists[i].validateParameters(validList, maxDepth);
2293 
2294  } catch (const Teuchos::Exceptions::InvalidParameterName& e) {
2295  std::string eString = e.what();
2296 
2297  // Parse name from: <Error, the parameter {name="smoothe: type",...>
2298  size_t nameStart = eString.find_first_of('"') + 1;
2299  size_t nameEnd = eString.find_first_of('"', nameStart);
2300  std::string name = eString.substr(nameStart, nameEnd - nameStart);
2301 
2302  size_t bestScore = 100;
2303  std::string bestName = "";
2304  for (ParameterList::ConstIterator it = validList.begin(); it != validList.end(); it++) {
2305  const std::string& pName = validList.name(it);
2306  this->GetOStream(Runtime1) << "| " << pName;
2307  size_t score = LevenshteinDistance(name.c_str(), name.length(), pName.c_str(), pName.length());
2308  this->GetOStream(Runtime1) << " -> " << score << std::endl;
2309  if (score < bestScore) {
2310  bestScore = score;
2311  bestName = pName;
2312  }
2313  }
2314  if (bestScore < 10 && bestName != "") {
2316  eString << "The parameter name \"" + name + "\" is not valid. Did you mean \"" + bestName << "\"?\n");
2317 
2318  } else {
2320  eString << "The parameter name \"" + name + "\" is not valid.\n");
2321  }
2322  }
2323  }
2324 }
2325 
2326 // =====================================================================================================
2327 // ==================================== FACTORY interpreter ============================================
2328 // =====================================================================================================
2329 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2331  SetFactoryParameterList(const ParameterList& constParamList) {
2332  // Create a non const copy of the parameter list
2333  // Working with a modifiable list is much much easier than with original one
2334  ParameterList paramList = constParamList;
2335 
2336  // Parameter List Parsing:
2337  // ---------
2338  // <ParameterList name="MueLu">
2339  // <ParameterList name="Matrix">
2340  // </ParameterList>
2341  if (paramList.isSublist("Matrix")) {
2342  blockSize_ = paramList.sublist("Matrix").get<int>("PDE equations", MasterList::getDefault<int>("number of equations"));
2343  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)
2344  }
2345 
2346  // create new FactoryFactory object if necessary
2347  if (factFact_ == Teuchos::null)
2348  factFact_ = Teuchos::rcp(new FactoryFactory());
2349 
2350  // Parameter List Parsing:
2351  // ---------
2352  // <ParameterList name="MueLu">
2353  // <ParameterList name="Factories"> <== call BuildFactoryMap() on this parameter list
2354  // ...
2355  // </ParameterList>
2356  // </ParameterList>
2357  FactoryMap factoryMap;
2358  FactoryManagerMap factoryManagers;
2359  if (paramList.isSublist("Factories"))
2360  this->BuildFactoryMap(paramList.sublist("Factories"), factoryMap, factoryMap, factoryManagers);
2361 
2362  // Parameter List Parsing:
2363  // ---------
2364  // <ParameterList name="MueLu">
2365  // <ParameterList name="Hierarchy">
2366  // <Parameter name="verbose" type="string" value="Warnings"/> <== get
2367  // <Parameter name="numDesiredLevel" type="int" value="10"/> <== get
2368  //
2369  // <ParameterList name="firstLevel"> <== parse first args and call BuildFactoryMap() on the rest of this parameter list
2370  // ...
2371  // </ParameterList>
2372  // </ParameterList>
2373  // </ParameterList>
2374  if (paramList.isSublist("Hierarchy")) {
2375  ParameterList hieraList = paramList.sublist("Hierarchy"); // copy because list temporally modified (remove 'id')
2376 
2377  // Get hierarchy options
2378  if (hieraList.isParameter("max levels")) {
2379  this->numDesiredLevel_ = hieraList.get<int>("max levels");
2380  hieraList.remove("max levels");
2381  }
2382 
2383  if (hieraList.isParameter("coarse: max size")) {
2384  this->maxCoarseSize_ = hieraList.get<int>("coarse: max size");
2385  hieraList.remove("coarse: max size");
2386  }
2387 
2388  if (hieraList.isParameter("repartition: rebalance P and R")) {
2389  this->doPRrebalance_ = hieraList.get<bool>("repartition: rebalance P and R");
2390  hieraList.remove("repartition: rebalance P and R");
2391  }
2392 
2393  if (hieraList.isParameter("transpose: use implicit")) {
2394  this->implicitTranspose_ = hieraList.get<bool>("transpose: use implicit");
2395  hieraList.remove("transpose: use implicit");
2396  }
2397 
2398  if (hieraList.isParameter("fuse prolongation and update")) {
2399  this->fuseProlongationAndUpdate_ = hieraList.get<bool>("fuse prolongation and update");
2400  hieraList.remove("fuse prolongation and update");
2401  }
2402 
2403  if (hieraList.isParameter("nullspace: suppress dimension check")) {
2404  this->suppressNullspaceDimensionCheck_ = hieraList.get<bool>("nullspace: suppress dimension check");
2405  hieraList.remove("nullspace: suppress dimension check");
2406  }
2407 
2408  if (hieraList.isParameter("number of vectors")) {
2409  this->sizeOfMultiVectors_ = hieraList.get<int>("number of vectors");
2410  hieraList.remove("number of vectors");
2411  }
2412 
2413  if (hieraList.isSublist("matvec params"))
2414  this->matvecParams_ = Teuchos::parameterList(hieraList.sublist("matvec params"));
2415 
2416  if (hieraList.isParameter("coarse grid correction scaling factor")) {
2417  this->scalingFactor_ = hieraList.get<double>("coarse grid correction scaling factor");
2418  hieraList.remove("coarse grid correction scaling factor");
2419  }
2420 
2421  // Translate cycle type parameter
2422  if (hieraList.isParameter("cycle type")) {
2423  std::map<std::string, CycleType> cycleMap;
2424  cycleMap["V"] = VCYCLE;
2425  cycleMap["W"] = WCYCLE;
2426 
2427  std::string cycleType = hieraList.get<std::string>("cycle type");
2428  TEUCHOS_TEST_FOR_EXCEPTION(cycleMap.count(cycleType) == 0, Exceptions::RuntimeError, "Invalid cycle type: \"" << cycleType << "\"");
2429  this->Cycle_ = cycleMap[cycleType];
2430  }
2431 
2432  if (hieraList.isParameter("W cycle start level")) {
2433  this->WCycleStartLevel_ = hieraList.get<int>("W cycle start level");
2434  }
2435 
2436  if (hieraList.isParameter("verbosity")) {
2437  std::string vl = hieraList.get<std::string>("verbosity");
2438  hieraList.remove("verbosity");
2439  this->verbosity_ = toVerbLevel(vl);
2440  }
2441 
2442  if (hieraList.isParameter("output filename"))
2443  VerboseObject::SetMueLuOFileStream(hieraList.get<std::string>("output filename"));
2444 
2445  if (hieraList.isParameter("dependencyOutputLevel"))
2446  this->graphOutputLevel_ = hieraList.get<int>("dependencyOutputLevel");
2447 
2448  // Check for the reuse case
2449  if (hieraList.isParameter("reuse"))
2451 
2452  if (hieraList.isSublist("DataToWrite")) {
2453  // TODO We should be able to specify any data. If it exists, write it.
2454  // TODO This would requires something like std::set<dataName, Array<int> >
2455  ParameterList foo = hieraList.sublist("DataToWrite");
2456  std::string dataName = "Matrices";
2457  if (foo.isParameter(dataName))
2458  this->matricesToPrint_["A"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2459  dataName = "Prolongators";
2460  if (foo.isParameter(dataName))
2461  this->matricesToPrint_["P"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2462  dataName = "Restrictors";
2463  if (foo.isParameter(dataName))
2464  this->matricesToPrint_["R"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2465  dataName = "D0";
2466  if (foo.isParameter(dataName))
2467  this->matricesToPrint_["D0"] = Teuchos::getArrayFromStringParameter<int>(foo, dataName);
2468  }
2469 
2470  // Get level configuration
2471  for (ParameterList::ConstIterator param = hieraList.begin(); param != hieraList.end(); ++param) {
2472  const std::string& paramName = hieraList.name(param);
2473 
2474  if (paramName != "DataToWrite" && hieraList.isSublist(paramName)) {
2475  ParameterList levelList = hieraList.sublist(paramName); // copy because list temporally modified (remove 'id')
2476 
2477  int startLevel = 0;
2478  if (levelList.isParameter("startLevel")) {
2479  startLevel = levelList.get<int>("startLevel");
2480  levelList.remove("startLevel");
2481  }
2482  int numDesiredLevel = 1;
2483  if (levelList.isParameter("numDesiredLevel")) {
2484  numDesiredLevel = levelList.get<int>("numDesiredLevel");
2485  levelList.remove("numDesiredLevel");
2486  }
2487 
2488  // Parameter List Parsing:
2489  // ---------
2490  // <ParameterList name="firstLevel">
2491  // <Parameter name="startLevel" type="int" value="0"/>
2492  // <Parameter name="numDesiredLevel" type="int" value="1"/>
2493  // <Parameter name="verbose" type="string" value="Warnings"/>
2494  //
2495  // [] <== call BuildFactoryMap() on the rest of the parameter list
2496  //
2497  // </ParameterList>
2498  FactoryMap levelFactoryMap;
2499  BuildFactoryMap(levelList, factoryMap, levelFactoryMap, factoryManagers);
2500 
2501  RCP<FactoryManager> m = rcp(new FactoryManager(levelFactoryMap));
2502  if (hieraList.isParameter("use kokkos refactor"))
2503  m->SetKokkosRefactor(hieraList.get<bool>("use kokkos refactor"));
2504 
2505  if (startLevel >= 0)
2506  this->AddFactoryManager(startLevel, numDesiredLevel, m);
2507  else
2508  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::ParameterListInterpreter():: invalid level id");
2509  } /* TODO: else { } */
2510  }
2511  }
2512 }
2513 
2514 // TODO: static?
2548 
2600 
2637 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2639  BuildFactoryMap(const ParameterList& paramList, const FactoryMap& factoryMapIn, FactoryMap& factoryMapOut, FactoryManagerMap& factoryManagers) const {
2640  for (ParameterList::ConstIterator param = paramList.begin(); param != paramList.end(); ++param) {
2641  const std::string& paramName = paramList.name(param); //< paramName contains the user chosen factory name (e.g., "smootherFact1")
2642  const Teuchos::ParameterEntry& paramValue = paramList.entry(param); //< for factories, paramValue should be either a list or just a MueLu Factory (e.g., TrilinosSmoother)
2643 
2644  // TODO: do not allow name of existing MueLu classes (can be tested using FactoryFactory)
2645 
2646  if (paramValue.isList()) {
2647  ParameterList paramList1 = Teuchos::getValue<ParameterList>(paramValue);
2648  if (paramList1.isParameter("factory")) { // default: just a factory definition
2649  // New Factory is a sublist with internal parameters and/or data dependencies
2650  TEUCHOS_TEST_FOR_EXCEPTION(paramList1.isParameter("dependency for") == true, Exceptions::RuntimeError,
2651  "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.");
2652 
2653  factoryMapOut[paramName] = factFact_->BuildFactory(paramValue, factoryMapIn, factoryManagers);
2654 
2655  } else if (paramList1.isParameter("dependency for")) { // add more data dependencies to existing factory
2656  TEUCHOS_TEST_FOR_EXCEPTION(paramList1.isParameter("factory") == true, Exceptions::RuntimeError,
2657  "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.");
2658 
2659  std::string factoryName = paramList1.get<std::string>("dependency for");
2660 
2661  RCP<const FactoryBase> factbase = factoryMapIn.find(factoryName /*paramName*/)->second; // access previously defined factory
2662  TEUCHOS_TEST_FOR_EXCEPTION(factbase.is_null() == true, Exceptions::RuntimeError,
2663  "MueLu::ParameterListInterpreter(): could not find factory " + factoryName + " in factory map. Did you define it before?");
2664 
2665  RCP<const Factory> factoryconst = Teuchos::rcp_dynamic_cast<const Factory>(factbase);
2666  RCP<Factory> factory = Teuchos::rcp_const_cast<Factory>(factoryconst);
2667 
2668  // Read the RCP<Factory> parameters of the class T
2669  RCP<const ParameterList> validParamList = factory->GetValidParameterList();
2670  for (ParameterList::ConstIterator vparam = validParamList->begin(); vparam != validParamList->end(); ++vparam) {
2671  const std::string& pName = validParamList->name(vparam);
2672 
2673  if (!paramList1.isParameter(pName)) {
2674  // Ignore unknown parameters
2675  continue;
2676  }
2677 
2678  if (validParamList->isType<RCP<const FactoryBase> >(pName)) {
2679  // Generate or get factory described by pName and set dependency
2680  RCP<const FactoryBase> generatingFact = factFact_->BuildFactory(paramList1.getEntry(pName), factoryMapIn, factoryManagers);
2681  factory->SetFactory(pName, generatingFact.create_weak());
2682 
2683  } else if (validParamList->isType<RCP<const ParameterList> >(pName)) {
2684  if (pName == "ParameterList") {
2685  // NOTE: we cannot use
2686  // subList = sublist(rcpFromRef(paramList), pName)
2687  // here as that would result in sublist also being a reference to a temporary object.
2688  // The resulting dereferencing in the corresponding factory would then segfault
2689  RCP<const ParameterList> subList = Teuchos::sublist(rcp(new ParameterList(paramList1)), pName);
2690  factory->SetParameter(pName, ParameterEntry(subList));
2691  }
2692  } else {
2693  factory->SetParameter(pName, paramList1.getEntry(pName));
2694  }
2695  }
2696 
2697  } else if (paramList1.isParameter("group")) { // definitiion of a factory group (for a factory manager)
2698  // Define a new (sub) FactoryManager
2699  std::string groupType = paramList1.get<std::string>("group");
2700  TEUCHOS_TEST_FOR_EXCEPTION(groupType != "FactoryManager", Exceptions::RuntimeError,
2701  "group must be of type \"FactoryManager\".");
2702 
2703  ParameterList groupList = paramList1; // copy because list temporally modified (remove 'id')
2704  groupList.remove("group");
2705 
2706  bool setKokkosRefactor = false;
2707  bool kokkosRefactor = useKokkos_;
2708  if (groupList.isParameter("use kokkos refactor")) {
2709  kokkosRefactor = groupList.get<bool>("use kokkos refactor");
2710  groupList.remove("use kokkos refactor");
2711  setKokkosRefactor = true;
2712  }
2713 
2714  FactoryMap groupFactoryMap;
2715  BuildFactoryMap(groupList, factoryMapIn, groupFactoryMap, factoryManagers);
2716 
2717  // do not store groupFactoryMap in factoryMapOut
2718  // Create a factory manager object from groupFactoryMap
2719  RCP<FactoryManager> m = rcp(new FactoryManager(groupFactoryMap));
2720  if (setKokkosRefactor)
2721  m->SetKokkosRefactor(kokkosRefactor);
2722  factoryManagers[paramName] = m;
2723 
2724  } else {
2725  this->GetOStream(Warnings0) << "Could not interpret parameter list " << paramList1 << std::endl;
2727  "XML Parameter list must either be of type \"factory\" or of type \"group\".");
2728  }
2729  } else {
2730  // default: just a factory (no parameter list)
2731  factoryMapOut[paramName] = factFact_->BuildFactory(paramValue, factoryMapIn, factoryManagers);
2732  }
2733  }
2734 }
2735 
2736 // =====================================================================================================
2737 // ======================================= MISC functions ==============================================
2738 // =====================================================================================================
2739 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2741  try {
2742  Matrix& A = dynamic_cast<Matrix&>(Op);
2743  if (A.IsFixedBlockSizeSet() && (A.GetFixedBlockSize() != blockSize_))
2744  this->GetOStream(Warnings0) << "Setting matrix block size to " << blockSize_ << " (value of the parameter in the list) "
2745  << "instead of " << A.GetFixedBlockSize() << " (provided matrix)." << std::endl
2746  << "You may want to check \"number of equations\" (or \"PDE equations\" for factory style list) parameter." << std::endl;
2747 
2748  A.SetFixedBlockSize(blockSize_, dofOffset_);
2749 
2750 #ifdef HAVE_MUELU_DEBUG
2751  MatrixUtils::checkLocalRowMapMatchesColMap(A);
2752 #endif // HAVE_MUELU_DEBUG
2753 
2754  } catch (std::bad_cast&) {
2755  this->GetOStream(Warnings0) << "Skipping setting block size as the operator is not a matrix" << std::endl;
2756  }
2757 }
2758 
2759 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2761  H.SetCycle(Cycle_);
2762  H.SetCycleStartLevel(WCycleStartLevel_);
2763  H.SetProlongatorScalingFactor(scalingFactor_);
2765 }
2766 
2767 static bool compare(const ParameterList& list1, const ParameterList& list2) {
2768  // First loop through and validate the parameters at this level.
2769  // In addition, we generate a list of sublists that we will search next
2770  for (ParameterList::ConstIterator it = list1.begin(); it != list1.end(); it++) {
2771  const std::string& name = it->first;
2772  const Teuchos::ParameterEntry& entry1 = it->second;
2773 
2774  const Teuchos::ParameterEntry* entry2 = list2.getEntryPtr(name);
2775  if (!entry2) // entry is not present in the second list
2776  return false;
2777  if (entry1.isList() && entry2->isList()) { // sublist check
2778  compare(Teuchos::getValue<ParameterList>(entry1), Teuchos::getValue<ParameterList>(*entry2));
2779  continue;
2780  }
2781  if (entry1.getAny(false) != entry2->getAny(false)) // entries have different types or different values
2782  return false;
2783  }
2784 
2785  return true;
2786 }
2787 
2788 static inline bool areSame(const ParameterList& list1, const ParameterList& list2) {
2789  return compare(list1, list2) && compare(list2, list1);
2790 }
2791 
2792 } // namespace MueLu
2793 
2794 #define MUELU_PARAMETERLISTINTERPRETER_SHORT
2795 #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