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