MueLu  Version of the Day
MueLu_ParameterListInterpreter_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_PARAMETERLISTINTERPRETER_DEF_HPP
47 #define MUELU_PARAMETERLISTINTERPRETER_DEF_HPP
48 
49 #include <Teuchos_XMLParameterListHelpers.hpp>
50 
51 #include <Xpetra_Matrix.hpp>
52 
53 #include "MueLu_ConfigDefs.hpp"
54 
56 
57 #include "MueLu_MasterList.hpp"
58 #include "MueLu_Level.hpp"
59 #include "MueLu_Hierarchy.hpp"
60 #include "MueLu_FactoryManager.hpp"
61 
62 #include "MueLu_AggregationExportFactory.hpp"
63 #include "MueLu_BrickAggregationFactory.hpp"
64 #include "MueLu_CoalesceDropFactory.hpp"
65 #include "MueLu_CoarseMapFactory.hpp"
66 #include "MueLu_ConstraintFactory.hpp"
67 #include "MueLu_CoordinatesTransferFactory.hpp"
68 #include "MueLu_CoupledAggregationFactory.hpp"
69 #include "MueLu_DirectSolver.hpp"
70 #include "MueLu_EminPFactory.hpp"
71 #include "MueLu_Exceptions.hpp"
72 #include "MueLu_FactoryFactory.hpp"
73 #include "MueLu_FilteredAFactory.hpp"
74 #include "MueLu_GenericRFactory.hpp"
75 #include "MueLu_LineDetectionFactory.hpp"
76 #include "MueLu_MasterList.hpp"
77 #include "MueLu_NullspaceFactory.hpp"
78 #include "MueLu_PatternFactory.hpp"
79 #include "MueLu_PgPFactory.hpp"
80 #include "MueLu_RAPFactory.hpp"
81 #include "MueLu_RebalanceAcFactory.hpp"
82 #include "MueLu_RebalanceTransferFactory.hpp"
83 #include "MueLu_RepartitionFactory.hpp"
84 #include "MueLu_SaPFactory.hpp"
85 #include "MueLu_SemiCoarsenPFactory.hpp"
86 #include "MueLu_SmootherFactory.hpp"
87 #include "MueLu_TentativePFactory.hpp"
88 #include "MueLu_TogglePFactory.hpp"
89 #include "MueLu_ToggleCoordinatesTransferFactory.hpp"
90 #include "MueLu_TransPFactory.hpp"
91 #include "MueLu_UncoupledAggregationFactory.hpp"
92 #include "MueLu_ZoltanInterface.hpp"
93 #include "MueLu_Zoltan2Interface.hpp"
94 
95 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
96 #include "MueLu_CoalesceDropFactory_kokkos.hpp"
97 #include "MueLu_CoarseMapFactory_kokkos.hpp"
98 #include "MueLu_CoordinatesTransferFactory_kokkos.hpp"
99 #include "MueLu_NullspaceFactory_kokkos.hpp"
100 #include "MueLu_SaPFactory_kokkos.hpp"
101 #include "MueLu_TentativePFactory_kokkos.hpp"
102 #include "MueLu_UncoupledAggregationFactory_kokkos.hpp"
103 #endif
104 
105 #ifdef HAVE_MUELU_MATLAB
106 #include "../matlab/src/MueLu_MatlabSmoother_decl.hpp"
107 #include "../matlab/src/MueLu_MatlabSmoother_def.hpp"
108 #include "../matlab/src/MueLu_TwoLevelMatlabFactory_decl.hpp"
109 #include "../matlab/src/MueLu_TwoLevelMatlabFactory_def.hpp"
110 #include "../matlab/src/MueLu_SingleLevelMatlabFactory_decl.hpp"
111 #include "../matlab/src/MueLu_SingleLevelMatlabFactory_def.hpp"
112 #endif
113 
114 // These code chunks should only be enabled once Tpetra supports proper graph
115 // reuse in MMM. At the moment, only Epetra does, while Tpetra throws
116 // #define REUSE_MATRIX_GRAPHS
117 
118 namespace MueLu {
119 
120  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121  ParameterListInterpreter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ParameterListInterpreter(ParameterList& paramList, Teuchos::RCP<const Teuchos::Comm<int> > comm, Teuchos::RCP<FactoryFactory> factFact) : factFact_(factFact) {
122 
123  if(paramList.isParameter("xml parameter file")){
124  std::string filename = paramList.get("xml parameter file","");
125  if(filename.length()!=0) {
126  if(comm.is_null()) throw Exceptions::RuntimeError("xml parameter file requires a valid comm");
127  Teuchos::ParameterList paramList2 = paramList;
128  Teuchos::updateParametersFromXmlFileAndBroadcast(filename, Teuchos::Ptr<Teuchos::ParameterList>(&paramList2),*comm);
129  SetParameterList(paramList2);
130  }
131  else
132  SetParameterList(paramList);
133  }
134  else
135  SetParameterList(paramList);
136  }
137 
138  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
139  ParameterListInterpreter<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ParameterListInterpreter(const std::string& xmlFileName, const Teuchos::Comm<int>& comm,Teuchos::RCP<FactoryFactory> factFact) : factFact_(factFact) {
140  ParameterList paramList;
141  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<ParameterList>(&paramList), comm);
142  SetParameterList(paramList);
143  }
144 
145  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
148  blockSize_ = 1;
149  dofOffset_ = 0;
150 
151  if (paramList.isSublist("Hierarchy")) {
152  SetFactoryParameterList(paramList);
153 
154  } else {
155  // The validator doesn't work correctly for non-serializable data (Hint: template parameters), so strip it out
156  ParameterList validList, nonSerialList;
157 
158  ExtractNonSerializableData(paramList, validList, nonSerialList);
159  Validate(validList);
160  SetEasyParameterList(paramList);
161  }
162  }
163 
164  // =====================================================================================================
165  // ====================================== EASY interpreter =============================================
166  // =====================================================================================================
168  static inline bool areSame(const ParameterList& list1, const ParameterList& list2);
169 
170  // Get value from one of the lists, or set it to default
171  // Use case: check for a parameter value in a level-specific sublist, then in a root level list;
172  // if it is absent from both, set it to default
173 #define MUELU_SET_VAR_2LIST(paramList, defaultList, paramName, paramType, varName) \
174  paramType varName; \
175  if (paramList.isParameter(paramName)) varName = paramList.get<paramType>(paramName); \
176  else if (defaultList.isParameter(paramName)) varName = defaultList.get<paramType>(paramName); \
177  else varName = MasterList::getDefault<paramType>(paramName);
178 
179 #define MUELU_TEST_AND_SET_VAR(paramList, paramName, paramType, varName) \
180  (paramList.isParameter(paramName) ? varName = paramList.get<paramType>(paramName), true : false)
181 
182  // Set parameter in a list if it is present in any of two lists
183  // User case: set factory specific parameter, first checking for a level-specific value, then cheking root level value
184 #define MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, paramName, paramType, listWrite) \
185  try { \
186  if (paramList .isParameter(paramName)) listWrite.set(paramName, paramList .get<paramType>(paramName)); \
187  else if (defaultList.isParameter(paramName)) listWrite.set(paramName, defaultList.get<paramType>(paramName)); \
188  } \
189  catch(Teuchos::Exceptions::InvalidParameterType) { \
190  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true, Teuchos::Exceptions::InvalidParameterType, \
191  "Error: parameter \"" << paramName << "\" must be of type " << Teuchos::TypeNameTraits<paramType>::name()); \
192  } \
193 
194 #define MUELU_TEST_PARAM_2LIST(paramList, defaultList, paramName, paramType, cmpValue) \
195  (cmpValue == ( \
196  paramList.isParameter(paramName) ? paramList .get<paramType>(paramName) : ( \
197  defaultList.isParameter(paramName) ? defaultList.get<paramType>(paramName) : \
198  MasterList::getDefault<paramType>(paramName) ) ) )
199 
200 #ifndef HAVE_MUELU_KOKKOS_REFACTOR
201 #define MUELU_KOKKOS_FACTORY(varName, oldFactory, newFactory) \
202  RCP<Factory> varName = rcp(new oldFactory());
203 #define MUELU_KOKKOS_FACTORY_NO_DECL(varName, oldFactory, newFactory) \
204  varName = rcp(new oldFactory());
205 #else
206 #define MUELU_KOKKOS_FACTORY(varName, oldFactory, newFactory) \
207  RCP<Factory> varName; \
208  if (!useKokkos) varName = rcp(new oldFactory()); \
209  else varName = rcp(new newFactory());
210 #define MUELU_KOKKOS_FACTORY_NO_DECL(varName, oldFactory, newFactory) \
211  if (!useKokkos) varName = rcp(new oldFactory()); \
212  else varName = rcp(new newFactory());
213 #endif
214 
215  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
217  ParameterList paramList;
218 
219  MUELU_SET_VAR_2LIST(constParamList, constParamList, "problem: type", std::string, problemType);
220  if (problemType != "unknown") {
221  paramList = *MasterList::GetProblemSpecificList(problemType);
222  paramList.setParameters(constParamList);
223  } else {
224  // Create a non const copy of the parameter list
225  // Working with a modifiable list is much much easier than with original one
226  paramList = constParamList;
227  }
228 
229  // Translate cycle type parameter
230  if (paramList.isParameter("cycle type")) {
231  std::map<std::string,CycleType> cycleMap;
232  cycleMap["V"] = VCYCLE;
233  cycleMap["W"] = WCYCLE;
234 
235  std::string cycleType = paramList.get<std::string>("cycle type");
236  TEUCHOS_TEST_FOR_EXCEPTION(cycleMap.count(cycleType) == 0, Exceptions::RuntimeError, "Invalid cycle type: \"" << cycleType << "\"");
237  Cycle_ = cycleMap[cycleType];
238  }
239 
240  this->maxCoarseSize_ = paramList.get<int> ("coarse: max size", MasterList::getDefault<int>("coarse: max size"));
241  this->numDesiredLevel_ = paramList.get<int> ("max levels", MasterList::getDefault<int>("max levels"));
242  blockSize_ = paramList.get<int> ("number of equations", MasterList::getDefault<int>("number of equations"));
243 
244  (void)MUELU_TEST_AND_SET_VAR(paramList, "debug: graph level", int, this->graphOutputLevel_);
245 
246  // Save level data
247  if (paramList.isSublist("export data")) {
248  ParameterList printList = paramList.sublist("export data");
249 
250  if (printList.isParameter("A"))
251  this->matricesToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "A");
252  if (printList.isParameter("P"))
253  this->prolongatorsToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "P");
254  if (printList.isParameter("R"))
255  this->restrictorsToPrint_ = Teuchos::getArrayFromStringParameter<int>(printList, "R");
256  }
257 
258  // Set verbosity parameter
259  {
260  std::map<std::string,MsgType> verbMap;
261  verbMap["none"] = None;
262  verbMap["low"] = Low;
263  verbMap["medium"] = Medium;
264  verbMap["high"] = High;
265  verbMap["extreme"] = Extreme;
266  verbMap["test"] = Test;
267 
268  MUELU_SET_VAR_2LIST(paramList, paramList, "verbosity", std::string, verbosityLevel);
269 
270  TEUCHOS_TEST_FOR_EXCEPTION(verbMap.count(verbosityLevel) == 0, Exceptions::RuntimeError,
271  "Invalid verbosity level: \"" << verbosityLevel << "\"");
272  this->verbosity_ = verbMap[verbosityLevel];
273  this->SetVerbLevel(this->verbosity_);
274  }
275 
276  // Detect if we need to transfer coordinates to coarse levels. We do that iff
277  // - we use "distance laplacian" dropping on some level, or
278  // - we use repartitioning on some level
279  // - we use brick aggregation
280  // This is not ideal, as we may have "repartition: enable" turned on by default
281  // and not present in the list, but it is better than nothing.
282  useCoordinates_ = false;
283  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "repartition: enable", bool, true) ||
284  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "distance laplacian") ||
285  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: type", std::string, "brick") ||
286  MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: export visualization data", bool, true)) {
287  useCoordinates_ = true;
288 
289  } else {
290  for (int levelID = 0; levelID < this->numDesiredLevel_; levelID++) {
291  std::string levelStr = "level " + toString(levelID);
292 
293  if (paramList.isSublist(levelStr)) {
294  const ParameterList& levelList = paramList.sublist(levelStr);
295 
296  if (MUELU_TEST_PARAM_2LIST(levelList, paramList, "repartition: enable", bool, true) ||
297  MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: drop scheme", std::string, "distance laplacian") ||
298  MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: type", std::string, "brick") ||
299  MUELU_TEST_PARAM_2LIST(levelList, paramList, "aggregation: export visualization data", bool, true)) {
300  useCoordinates_ = true;
301  break;
302  }
303  }
304  }
305  }
306 
307  // Detect if we do implicit P and R rebalance
308  changedPRrebalance_ = false;
309  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "repartition: enable", bool, true))
310  changedPRrebalance_ = MUELU_TEST_AND_SET_VAR(paramList, "repartition: rebalance P and R", bool, this->doPRrebalance_);
311 
312  // Detect if we use implicit transpose
313  changedImplicitTranspose_ = MUELU_TEST_AND_SET_VAR(paramList, "transpose: use implicit", bool, this->implicitTranspose_);
314 
315  // Create default manager
316  // FIXME: should it be here, or higher up
317  RCP<FactoryManager> defaultManager = rcp(new FactoryManager());
318  defaultManager->SetVerbLevel(this->verbosity_);
319 
320  // We will ignore keeps0
321  std::vector<keep_pair> keeps0;
322  UpdateFactoryManager(paramList, ParameterList(), *defaultManager, 0/*levelID*/, keeps0);
323 
324  // Create level specific factory managers
325  for (int levelID = 0; levelID < this->numDesiredLevel_; levelID++) {
326  // Note, that originally if there were no level specific parameters, we
327  // simply copied the defaultManager However, with the introduction of
328  // levelID to UpdateFactoryManager (required for reuse), we can no longer
329  // guarantee that the kept variables are the same for each level even if
330  // dependency structure does not change.
331  RCP<FactoryManager> levelManager = rcp(new FactoryManager(*defaultManager));
332  levelManager->SetVerbLevel(defaultManager->GetVerbLevel());
333 
334  std::vector<keep_pair> keeps;
335  if (paramList.isSublist("level " + toString(levelID))) {
336  // We do this so the parameters on the level get flagged correctly as "used"
337  ParameterList& levelList = paramList.sublist("level " + toString(levelID), true/*mustAlreadyExist*/);
338  UpdateFactoryManager(levelList, paramList, *levelManager, levelID, keeps);
339 
340  } else {
341  ParameterList levelList;
342  UpdateFactoryManager(levelList, paramList, *levelManager, levelID, keeps);
343  }
344 
345  this->keep_[levelID] = keeps;
346  this->AddFactoryManager(levelID, 1, levelManager);
347  }
348 
349  // FIXME: parameters passed to packages, like Ifpack2, are not touched by us, resulting in "[unused]" flag
350  // being displayed. On the other hand, we don't want to simply iterate through them touching. I don't know
351  // what a good solution looks like
352  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "print initial parameters", bool, true))
353  this->GetOStream(static_cast<MsgType>(Runtime1), 0) << paramList << std::endl;
354 
355  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "print unused parameters", bool, true)) {
356  // Check unused parameters
357  ParameterList unusedParamList;
358 
359  // Check for unused parameters that aren't lists
360  for (ParameterList::ConstIterator it = paramList.begin(); it != paramList.end(); it++) {
361  const ParameterEntry& entry = paramList.entry(it);
362 
363  if (!entry.isList() && !entry.isUsed())
364  unusedParamList.setEntry(paramList.name(it), entry);
365  }
366 
367  // Check for unused parameters in level-specific sublists
368  for (int levelID = 0; levelID < this->numDesiredLevel_; levelID++) {
369  std::string levelStr = "level " + toString(levelID);
370 
371  if (paramList.isSublist(levelStr)) {
372  const ParameterList& levelList = paramList.sublist(levelStr);
373 
374  for (ParameterList::ConstIterator itr = levelList.begin(); itr != levelList.end(); ++itr) {
375  const ParameterEntry& entry = levelList.entry(itr);
376 
377  if (!entry.isList() && !entry.isUsed())
378  unusedParamList.sublist(levelStr).setEntry(levelList.name(itr), entry);
379  }
380  }
381  }
382 
383  if (unusedParamList.numParams() > 0) {
384  std::ostringstream unusedParamsStream;
385  int indent = 4;
386  unusedParamList.print(unusedParamsStream, indent);
387 
388  this->GetOStream(Warnings1) << "The following parameters were not used:\n" << unusedParamsStream.str() << std::endl;
389  }
390  }
391 
392  }
393 
394  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
396  const ParameterList& defaultList, FactoryManager& manager, int levelID, std::vector<keep_pair>& keeps) const {
397  // NOTE: Factory::SetParameterList must be called prior to Factory::SetFactory, as
398  // SetParameterList sets default values for non mentioned parameters, including factories
399 
400  // shortcut
401  if (paramList.numParams() == 0 && defaultList.numParams() > 0)
402  paramList = ParameterList(defaultList);
403 
404  MUELU_SET_VAR_2LIST(paramList, defaultList, "reuse: type", std::string, reuseType);
405  TEUCHOS_TEST_FOR_EXCEPTION(reuseType != "none" && reuseType != "tP" && reuseType != "RP" && reuseType != "emin" && reuseType != "RAP" && reuseType != "full",
406  Exceptions::RuntimeError, "Unknown \"reuse: type\" value: \"" << reuseType << "\". Please consult User's Guide.");
407 
408  MUELU_SET_VAR_2LIST(paramList, defaultList, "multigrid algorithm", std::string, multigridAlgo);
409  TEUCHOS_TEST_FOR_EXCEPTION(multigridAlgo != "unsmoothed" && multigridAlgo != "sa" && multigridAlgo != "pg" && multigridAlgo != "emin" && multigridAlgo != "matlab", Exceptions::RuntimeError, "Unknown \"multigrid algorithm\" value: \"" << multigridAlgo << "\". Please consult User's Guide.");
410  #ifndef HAVE_MUELU_MATLAB
411  if(multigridAlgo == "matlab")
412  throw std::runtime_error("Cannot use matlab for multigrid algorithm - MueLu was not configured with MATLAB support.");
413  #endif
414  // Only some combinations of reuse and multigrid algorithms are tested, all
415  // other are considered invalid at the moment
416  if (reuseType == "none" || reuseType == "RP" || reuseType == "RAP") {
417  // This works for all kinds of multigrid algorithms
418 
419  } else if (reuseType == "tP" && (multigridAlgo != "sa" && multigridAlgo != "unsmoothed")) {
420  reuseType = "none";
421  this->GetOStream(Warnings0) << "Ignoring \"tP\" reuse option as it is only compatible with \"sa\", or \"unsmoothed\" multigrid algorithms" << std::endl;
422 
423  } else if (reuseType == "emin" && multigridAlgo != "emin") {
424  reuseType = "none";
425  this->GetOStream(Warnings0) << "Ignoring \"emin\" reuse option it is only compatible with \"emin\" multigrid algorithm" << std::endl;
426  }
427 
428  MUELU_SET_VAR_2LIST(paramList, defaultList, "use kokkos refactor", bool, useKokkos);
429  (void) useKokkos;
430 
431  // == Non-serializable data ===
432  // Check both the parameter and the type
433  bool have_userA = false, have_userP = false, have_userR = false, have_userNS = false, have_userCO = false;
434  if (paramList.isParameter("A") && !paramList.get<RCP<Matrix> > ("A") .is_null()) have_userA = true;
435  if (paramList.isParameter("P") && !paramList.get<RCP<Matrix> > ("P") .is_null()) have_userP = true;
436  if (paramList.isParameter("R") && !paramList.get<RCP<Matrix> > ("R") .is_null()) have_userR = true;
437  if (paramList.isParameter("Nullspace") && !paramList.get<RCP<MultiVector> >("Nullspace") .is_null()) have_userNS = true;
438  if (paramList.isParameter("Coordinates") && !paramList.get<RCP<MultiVector> >("Coordinates").is_null()) have_userCO = true;
439 
440  // === Smoothing ===
441  // FIXME: should custom smoother check default list too?
442  bool isCustomSmoother =
443  paramList.isParameter("smoother: pre or post") ||
444  paramList.isParameter("smoother: type") || paramList.isParameter("smoother: pre type") || paramList.isParameter("smoother: post type") ||
445  paramList.isSublist ("smoother: params") || paramList.isSublist ("smoother: pre params") || paramList.isSublist ("smoother: post params") ||
446  paramList.isParameter("smoother: sweeps") || paramList.isParameter("smoother: pre sweeps") || paramList.isParameter("smoother: post sweeps") ||
447  paramList.isParameter("smoother: overlap") || paramList.isParameter("smoother: pre overlap") || paramList.isParameter("smoother: post overlap");
448  MUELU_SET_VAR_2LIST(paramList, defaultList, "smoother: pre or post", std::string, PreOrPost);
449  if (PreOrPost == "none") {
450  manager.SetFactory("Smoother", Teuchos::null);
451 
452  } else if (isCustomSmoother) {
453  // FIXME: get default values from the factory
454  // NOTE: none of the smoothers at the moment use parameter validation framework, so we
455  // cannot get the default values from it.
456 #define TEST_MUTUALLY_EXCLUSIVE(arg1,arg2) \
457  TEUCHOS_TEST_FOR_EXCEPTION(paramList.isParameter(#arg1) && paramList.isParameter(#arg2), \
458  Exceptions::InvalidArgument, "You cannot specify both \""#arg1"\" and \""#arg2"\"");
459 #define TEST_MUTUALLY_EXCLUSIVE_S(arg1,arg2) \
460  TEUCHOS_TEST_FOR_EXCEPTION(paramList.isSublist(#arg1) && paramList.isSublist(#arg2), \
461  Exceptions::InvalidArgument, "You cannot specify both \""#arg1"\" and \""#arg2"\"");
462 
463  TEST_MUTUALLY_EXCLUSIVE ("smoother: type", "smoother: pre type");
464  TEST_MUTUALLY_EXCLUSIVE ("smoother: type", "smoother: post type");
465  TEST_MUTUALLY_EXCLUSIVE ("smoother: sweeps", "smoother: pre sweeps");
466  TEST_MUTUALLY_EXCLUSIVE ("smoother: sweeps", "smoother: post sweeps");
467  TEST_MUTUALLY_EXCLUSIVE ("smoother: overlap", "smoother: pre overlap");
468  TEST_MUTUALLY_EXCLUSIVE ("smoother: overlap", "smoother: post overlap");
469  TEST_MUTUALLY_EXCLUSIVE_S("smoother: params", "smoother: pre params");
470  TEST_MUTUALLY_EXCLUSIVE_S("smoother: params", "smoother: post params");
471  TEUCHOS_TEST_FOR_EXCEPTION(PreOrPost == "both" && (paramList.isParameter("smoother: pre type") != paramList.isParameter("smoother: post type")),
472  Exceptions::InvalidArgument, "You must specify both \"smoother: pre type\" and \"smoother: post type\"");
473 
474  // Default values
475  int overlap = 0;
476  ParameterList defaultSmootherParams;
477  defaultSmootherParams.set("relaxation: type", "Symmetric Gauss-Seidel");
478  defaultSmootherParams.set("relaxation: sweeps", Teuchos::OrdinalTraits<LO>::one());
479  defaultSmootherParams.set("relaxation: damping factor", Teuchos::ScalarTraits<Scalar>::one());
480 
481  RCP<SmootherFactory> preSmoother = Teuchos::null, postSmoother = Teuchos::null;
482  std::string preSmootherType, postSmootherType;
483  ParameterList preSmootherParams, postSmootherParams;
484 
485  if (paramList.isParameter("smoother: overlap"))
486  overlap = paramList.get<int>("smoother: overlap");
487 
488  if (PreOrPost == "pre" || PreOrPost == "both") {
489  if (paramList.isParameter("smoother: pre type")) {
490  preSmootherType = paramList.get<std::string>("smoother: pre type");
491  } else {
492  MUELU_SET_VAR_2LIST(paramList, defaultList, "smoother: type", std::string, preSmootherTypeTmp);
493  preSmootherType = preSmootherTypeTmp;
494  }
495  if (paramList.isParameter("smoother: pre overlap"))
496  overlap = paramList.get<int>("smoother: pre overlap");
497 
498  if (paramList.isSublist("smoother: pre params"))
499  preSmootherParams = paramList.sublist("smoother: pre params");
500  else if (paramList.isSublist("smoother: params"))
501  preSmootherParams = paramList.sublist("smoother: params");
502  else if (defaultList.isSublist("smoother: params"))
503  preSmootherParams = defaultList.sublist("smoother: params");
504  else if (preSmootherType == "RELAXATION")
505  preSmootherParams = defaultSmootherParams;
506 #ifdef HAVE_MUELU_MATLAB
507  if(preSmootherType == "matlab")
508  preSmoother = rcp(new SmootherFactory(rcp(new MatlabSmoother<Scalar,LocalOrdinal, GlobalOrdinal, Node>(preSmootherParams))));
509  else
510 #endif
511  preSmoother = rcp(new SmootherFactory(rcp(new TrilinosSmoother(preSmootherType, preSmootherParams, overlap))));
512  }
513 
514  if (PreOrPost == "post" || PreOrPost == "both") {
515  if (paramList.isParameter("smoother: post type"))
516  postSmootherType = paramList.get<std::string>("smoother: post type");
517  else {
518  MUELU_SET_VAR_2LIST(paramList, defaultList, "smoother: type", std::string, postSmootherTypeTmp);
519  postSmootherType = postSmootherTypeTmp;
520  }
521 
522  if (paramList.isSublist("smoother: post params"))
523  postSmootherParams = paramList.sublist("smoother: post params");
524  else if (paramList.isSublist("smoother: params"))
525  postSmootherParams = paramList.sublist("smoother: params");
526  else if (defaultList.isSublist("smoother: params"))
527  postSmootherParams = defaultList.sublist("smoother: params");
528  else if (postSmootherType == "RELAXATION")
529  postSmootherParams = defaultSmootherParams;
530  if (paramList.isParameter("smoother: post overlap"))
531  overlap = paramList.get<int>("smoother: post overlap");
532 
533  if (postSmootherType == preSmootherType && areSame(preSmootherParams, postSmootherParams))
534  postSmoother = preSmoother;
535  else
536 #ifdef HAVE_MUELU_MATLAB
537  if(postSmootherType == "matlab")
538  postSmoother = rcp(new SmootherFactory(rcp(new MatlabSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>(postSmootherParams))));
539  else
540 #endif
541  postSmoother = rcp(new SmootherFactory(rcp(new TrilinosSmoother(postSmootherType, postSmootherParams, overlap))));
542  }
543 
544  if (preSmoother == postSmoother)
545  manager.SetFactory("Smoother", preSmoother);
546  else {
547  manager.SetFactory("PreSmoother", preSmoother);
548  manager.SetFactory("PostSmoother", postSmoother);
549  }
550  }
551  if ((reuseType == "RAP" && levelID) || (reuseType == "full")) {
552  // The difference between "RAP" and "full" is keeping smoothers. However,
553  // as in both cases we keep coarse matrices, we do not need to update
554  // coarse smoothers. On the other hand, if a user changes fine level
555  // matrix, "RAP" would update the fine level smoother, while "full" would
556  // not
557  keeps.push_back(keep_pair("PreSmoother", manager.GetFactory("PreSmoother") .get()));
558  keeps.push_back(keep_pair("PostSmoother", manager.GetFactory("PostSmoother").get()));
559  }
560 
561  // === Coarse solver ===
562  // FIXME: should custom coarse solver check default list too?
563  bool isCustomCoarseSolver =
564  paramList.isParameter("coarse: type") ||
565  paramList.isParameter("coarse: params");
566  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "coarse: type", std::string, "none")) {
567  manager.SetFactory("CoarseSolver", Teuchos::null);
568 
569  } else if (isCustomCoarseSolver) {
570  // FIXME: get default values from the factory
571  // NOTE: none of the smoothers at the moment use parameter validation framework, so we
572  // cannot get the default values from it.
573  MUELU_SET_VAR_2LIST(paramList, defaultList, "coarse: type", std::string, coarseType);
574 
575  int overlap = 0;
576  if (paramList.isParameter("coarse: overlap"))
577  overlap = paramList.get<int>("coarse: overlap");
578 
579  ParameterList coarseParams;
580  if (paramList.isSublist("coarse: params"))
581  coarseParams = paramList.sublist("coarse: params");
582  else if (defaultList.isSublist("coarse: params"))
583  coarseParams = defaultList.sublist("coarse: params");
584 
585  RCP<SmootherPrototype> coarseSmoother;
586  // TODO: this is not a proper place to check. If we consider direct solver to be a special
587  // case of smoother, we would like to unify Amesos and Ifpack2 smoothers in src/Smoothers, and
588  // have a single factory responsible for those. Then, this check would belong there.
589  if (coarseType == "RELAXATION" || coarseType == "CHEBYSHEV" ||
590  coarseType == "ILUT" || coarseType == "ILU" || coarseType == "RILUK" || coarseType == "SCHWARZ" ||
591  coarseType == "Amesos" ||
592  coarseType == "LINESMOOTHING_BANDEDRELAXATION" ||
593  coarseType == "LINESMOOTHING_BANDED_RELAXATION" ||
594  coarseType == "LINESMOOTHING_BANDED RELAXATION")
595  coarseSmoother = rcp(new TrilinosSmoother(coarseType, coarseParams, overlap));
596  else
597 #ifdef HAVE_MUELU_MATLAB
598  if(coarseType == "matlab")
599  coarseSmoother = rcp(new MatlabSmoother<Scalar,LocalOrdinal, GlobalOrdinal, Node>(coarseParams));
600  else
601 #endif
602  coarseSmoother = rcp(new DirectSolver(coarseType, coarseParams));
603 
604  manager.SetFactory("CoarseSolver", rcp(new SmootherFactory(coarseSmoother)));
605  }
606  if ((reuseType == "RAP" && levelID) || (reuseType == "full")) {
607  // We do keep_pair("PreSmoother", // manager.GetFactory("CoarseSolver").get())
608  // as the coarse solver factory is in fact a smoothing factory, so the
609  // only pieces of data it generates are PreSmoother and PostSmoother
610  keeps.push_back(keep_pair("PreSmoother", manager.GetFactory("CoarseSolver").get()));
611  }
612 
613  // === Aggregation ===
614  // Aggregation graph
615  RCP<Factory> dropFactory;
616 
617  if (MUELU_TEST_PARAM_2LIST(paramList, paramList, "aggregation: drop scheme", std::string, "matlab")) {
618 #ifdef HAVE_MUELU_MATLAB
620  ParameterList socParams = paramList.sublist("strength-of-connection: params");
621  dropFactory->SetParameterList(socParams);
622 #else
623  throw std::runtime_error("Cannot use MATLAB evolutionary strength-of-connection - MueLu was not configured with MATLAB support.");
624 #endif
625  } else {
626  MUELU_KOKKOS_FACTORY_NO_DECL(dropFactory, CoalesceDropFactory, CoalesceDropFactory_kokkos);
627  ParameterList dropParams;
628  dropParams.set("lightweight wrap", true);
629  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop scheme", std::string, dropParams);
630  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol", double, dropParams);
631  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: Dirichlet threshold", double, dropParams);
632  dropFactory->SetParameterList(dropParams);
633  }
634  manager.SetFactory("Graph", dropFactory);
635 
636  // Aggregation scheme
637  MUELU_SET_VAR_2LIST(paramList, defaultList, "aggregation: type", std::string, aggType);
638  TEUCHOS_TEST_FOR_EXCEPTION(aggType != "uncoupled" && aggType != "coupled" && aggType != "brick" && aggType != "matlab",
639  Exceptions::RuntimeError, "Unknown aggregation algorithm: \"" << aggType << "\". Please consult User's Guide.");
640  #ifndef HAVE_MUELU_MATLAB
641  if(aggType == "matlab")
642  throw std::runtime_error("Cannot use MATLAB aggregation - MueLu was not configured with MATLAB support.");
643  #endif
644  RCP<Factory> aggFactory;
645  if (aggType == "uncoupled") {
646  MUELU_KOKKOS_FACTORY_NO_DECL(aggFactory, UncoupledAggregationFactory, UncoupledAggregationFactory_kokkos);
647  ParameterList aggParams;
648  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: mode", std::string, aggParams);
649  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: ordering", std::string, aggParams);
650  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: min agg size", int, aggParams);
651  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: max agg size", int, aggParams);
652  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: max selected neighbors", int, aggParams);
653  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 1", bool, aggParams);
654  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 2a", bool, aggParams);
655  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 2b", bool, aggParams);
656  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: enable phase 3", bool, aggParams);
657  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: preserve Dirichlet points", bool, aggParams);
658  aggFactory->SetParameterList(aggParams);
659  // make sure that the aggregation factory has all necessary data
660  aggFactory->SetFactory("DofsPerNode", manager.GetFactory("Graph"));
661  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
662  } else if (aggType == "coupled") {
663  aggFactory = rcp(new CoupledAggregationFactory());
664  aggFactory->SetFactory("Graph", manager.GetFactory("Graph"));
665  } else if (aggType == "brick") {
666  aggFactory = rcp(new BrickAggregationFactory());
667  ParameterList aggParams;
668  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick x size", int, aggParams);
669  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick y size", int, aggParams);
670  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick z size", int, aggParams);
671  aggFactory->SetParameterList(aggParams);
672  if (levelID > 1) {
673  // We check for levelID > 0, as in the interpreter aggFactory for
674  // levelID really corresponds to level 0. Managers are clunky, as they
675  // contain factories for two different levels
676  aggFactory->SetFactory("Coordinates", this->GetFactoryManager(levelID-1)->GetFactory("Coordinates"));
677  }
678  }
679 #ifdef HAVE_MUELU_MATLAB
680  else if(aggType == "matlab") {
681  ParameterList aggParams = paramList.sublist("aggregation: params");
683  aggFactory->SetParameterList(aggParams);
684  }
685 #endif
686  manager.SetFactory("Aggregates", aggFactory);
687 
688  // Coarse map
689  MUELU_KOKKOS_FACTORY(coarseMap, CoarseMapFactory, CoarseMapFactory_kokkos);
690  coarseMap->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
691  manager.SetFactory("CoarseMap", coarseMap);
692 
693  // Tentative P
694  MUELU_KOKKOS_FACTORY(Ptent, TentativePFactory, TentativePFactory_kokkos);
695  Ptent->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
696  Ptent->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
697 
698  manager.SetFactory("Ptent", Ptent);
699 
700  if (reuseType == "tP") {
701  keeps.push_back(keep_pair("Nullspace", manager.GetFactory("Ptent").get()));
702  keeps.push_back(keep_pair("P", manager.GetFactory("Ptent").get()));
703  }
704 
705  // Nullspace
706  MUELU_KOKKOS_FACTORY(nullSpace, NullspaceFactory, NullspaceFactory_kokkos);
707  if (!have_userNS) {
708  nullSpace->SetFactory("Nullspace", manager.GetFactory("Ptent"));
709  manager.SetFactory("Nullspace", nullSpace);
710  }
711 
712  // === Prolongation ===
713  TEUCHOS_TEST_FOR_EXCEPTION(multigridAlgo != "unsmoothed" && multigridAlgo != "sa" && multigridAlgo != "pg" && multigridAlgo != "emin" && multigridAlgo != "matlab", Exceptions::RuntimeError, "Unknown multigrid algorithm: \"" << multigridAlgo << "\". Please consult User's Guide.");
714  #ifndef HAVE_MUELU_MATLAB
715  if(multigridAlgo == "matlab")
716  throw std::runtime_error("Cannot use MATLAB prolongator factory - MueLu was not configured with MATLAB support.");
717  #endif
718  if (have_userP) {
719  // User prolongator
720  manager.SetFactory("P", NoFactory::getRCP());
721  } else if (multigridAlgo == "unsmoothed") {
722  // Unsmoothed aggregation
723  manager.SetFactory("P", Ptent);
724  } else if (multigridAlgo == "sa") {
725  // Smoothed aggregation
726  MUELU_KOKKOS_FACTORY(P, SaPFactory, SaPFactory_kokkos);
727  ParameterList Pparams;
728  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: damping factor", double, Pparams);
729 #if REUSE_MATRIX_GRAPHS
730  if (reuseType == "tP" && MUELU_TEST_PARAM_2LIST(paramList, defaultList, "sa: use filtered matrix", bool, false)) {
731  // Pattern can only be reuse when we don't use filtered matrix, as
732  // otherwise graph can potentially change
733  Pparams.set("Keep AP Pattern", true);
734  }
735 #endif
736  P->SetParameterList(Pparams);
737 
738  // Filtering
739  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "sa: use filtered matrix", bool, true)) {
740  RCP<FilteredAFactory> filterFactory = rcp(new FilteredAFactory());
741  ParameterList fParams;
742  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: use lumping", bool, fParams);
743  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse graph", bool, fParams);
744  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "filtered matrix: reuse eigenvalue", bool, fParams);
745  filterFactory->SetParameterList(fParams);
746  filterFactory->SetFactory("Graph", manager.GetFactory("Graph"));
747  // I'm not sure why we need this line. See comments for DofsPerNode for UncoupledAggregation above
748  filterFactory->SetFactory("Filtering", manager.GetFactory("Graph"));
749  P->SetFactory("A", filterFactory);
750  }
751 
752  P->SetFactory("P", manager.GetFactory("Ptent"));
753  manager.SetFactory("P", P);
754 
755  } else if (multigridAlgo == "emin") {
756  MUELU_SET_VAR_2LIST(paramList, defaultList, "emin: pattern", std::string, patternType);
757  TEUCHOS_TEST_FOR_EXCEPTION(patternType != "AkPtent", Exceptions::InvalidArgument,
758  "Invalid pattern name: \"" << patternType << "\". Valid options: \"AkPtent\"");
759  // Pattern
760  RCP<PatternFactory> patternFactory = rcp(new PatternFactory());
761  ParameterList patternParams;
762  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: pattern order", int, patternParams);
763  patternFactory->SetParameterList(patternParams);
764  patternFactory->SetFactory("P", manager.GetFactory("Ptent"));
765  manager.SetFactory("Ppattern", patternFactory);
766 
767  // Constraint
768  RCP<ConstraintFactory> constraintFactory = rcp(new ConstraintFactory());
769  constraintFactory->SetFactory("Ppattern", manager.GetFactory("Ppattern"));
770  constraintFactory->SetFactory("CoarseNullspace", manager.GetFactory("Ptent"));
771  manager.SetFactory("Constraint", constraintFactory);
772 
773  // Energy minimization
774  RCP<EminPFactory> P = rcp(new EminPFactory());
775  ParameterList Pparams;
776  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: num iterations", int, Pparams);
777  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: iterative method", std::string, Pparams);
778  if (reuseType == "emin") {
779  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "emin: num reuse iterations", int, Pparams);
780  Pparams.set("Keep P0", true);
781  Pparams.set("Keep Constraint0", true);
782  }
783  P->SetParameterList(Pparams);
784  P->SetFactory("P", manager.GetFactory("Ptent"));
785  P->SetFactory("Constraint", manager.GetFactory("Constraint"));
786  manager.SetFactory("P", P);
787 
788  } else if (multigridAlgo == "pg") {
789  TEUCHOS_TEST_FOR_EXCEPTION(this->implicitTranspose_, Exceptions::RuntimeError,
790  "Implicit transpose not supported with Petrov-Galerkin smoothed transfer operators: Set \"transpose: use implicit\" to false!\n" \
791  "Petrov-Galerkin transfer operator smoothing for non-symmetric problems requires a separate handling of the restriction operator which " \
792  "does not allow the usage of implicit transpose easily.");
793 
794  // Petrov-Galerkin
795  RCP<PgPFactory> P = rcp(new PgPFactory());
796  P->SetFactory("P", manager.GetFactory("Ptent"));
797  manager.SetFactory("P", P);
798  }
799 #ifdef HAVE_MUELU_MATLAB
800  else if(multigridAlgo == "matlab") {
801  ParameterList Pparams = paramList.sublist("transfer: params");
802  RCP<TwoLevelMatlabFactory<Scalar,LocalOrdinal, GlobalOrdinal, Node> > P = rcp(new TwoLevelMatlabFactory<Scalar,LocalOrdinal, GlobalOrdinal, Node>());
803  P->SetParameterList(Pparams);
804  P->SetFactory("P",manager.GetFactory("Ptent"));
805  manager.SetFactory("P", P);
806  }
807 #endif
808 
809  // === Semi-coarsening ===
810  RCP<SemiCoarsenPFactory> semicoarsenFactory = Teuchos::null;
811  if (paramList.isParameter("semicoarsen: number of levels") &&
812  paramList.get<int>("semicoarsen: number of levels") > 0) {
813 
814  ParameterList togglePParams;
815  ParameterList semicoarsenPParams;
816  ParameterList linedetectionParams;
817  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: number of levels", int, togglePParams);
818  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "semicoarsen: coarsen rate", int, semicoarsenPParams);
819  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "linedetection: orientation", std::string, linedetectionParams);
820  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "linedetection: num layers", int, linedetectionParams);
821 
822  semicoarsenFactory = rcp(new SemiCoarsenPFactory());
823  RCP<LineDetectionFactory> linedetectionFactory = rcp(new LineDetectionFactory());
824  RCP<TogglePFactory> togglePFactory = rcp(new TogglePFactory());
825 
826  linedetectionFactory->SetParameterList(linedetectionParams);
827  semicoarsenFactory->SetParameterList(semicoarsenPParams);
828  togglePFactory->SetParameterList(togglePParams);
829  togglePFactory->AddCoarseNullspaceFactory(semicoarsenFactory);
830  togglePFactory->AddProlongatorFactory(semicoarsenFactory);
831  togglePFactory->AddPtentFactory(semicoarsenFactory);
832  togglePFactory->AddCoarseNullspaceFactory(manager.GetFactory("Ptent"));
833  togglePFactory->AddProlongatorFactory(manager.GetFactory("P"));
834  togglePFactory->AddPtentFactory(manager.GetFactory("Ptent"));
835 
836  manager.SetFactory("CoarseNumZLayers", linedetectionFactory);
837  manager.SetFactory("LineDetection_Layers", linedetectionFactory);
838  manager.SetFactory("LineDetection_VertLineIds", linedetectionFactory);
839 
840  manager.SetFactory("P", togglePFactory);
841  manager.SetFactory("Ptent", togglePFactory);
842  manager.SetFactory("Nullspace", togglePFactory);
843  }
844 
845  // === Restriction ===
846  if (!this->implicitTranspose_) {
847  MUELU_SET_VAR_2LIST(paramList, defaultList, "problem: symmetric", bool, isSymmetric);
848  if (isSymmetric == false && (multigridAlgo == "unsmoothed" || multigridAlgo == "emin")) {
849  this->GetOStream(Warnings0) << "Switching \"problem: symmetric\" parameter to symmetric as multigrid algorithm. " << multigridAlgo << " is primarily supposed to be used for symmetric problems." << std::endl << std::endl;
850  this->GetOStream(Warnings0) << "Please note: if you are using \"unsmoothed\" transfer operators the \"problem: symmetric\" parameter has no real mathematical meaning, i.e. you can use it for non-symmetric" << std::endl;
851  this->GetOStream(Warnings0) << "problems, too. With \"problem: symmetric\"=\"symmetric\" you can use implicit transpose for building the restriction operators which may drastically reduce the amount of consumed memory." << std::endl;
852  isSymmetric = true;
853  }
854  TEUCHOS_TEST_FOR_EXCEPTION(multigridAlgo == "pg" && isSymmetric == true, Exceptions::RuntimeError,
855  "Petrov-Galerkin smoothed transfer operators are only allowed for non-symmetric problems: Set \"problem: symmetric\" to false!\n" \
856  "While PG smoothed transfer operators generally would also work for symmetric problems this is an unusual use case. " \
857  "You can use the factory-based xml interface though if you need PG-AMG for symmetric problems.");
858 
859  if (have_userR) {
860  manager.SetFactory("R", NoFactory::getRCP());
861  } else {
862  RCP<Factory> R;
863  if (isSymmetric) R = rcp(new TransPFactory());
864  else R = rcp(new GenericRFactory());
865 
866  R->SetFactory("P", manager.GetFactory("P"));
867  manager.SetFactory("R", R);
868  }
869 
870  } else {
871  manager.SetFactory("R", Teuchos::null);
872  }
873 
874  // === RAP ===
875  RCP<RAPFactory> RAP;
876  if (have_userA) {
877  manager.SetFactory("A", NoFactory::getRCP());
878 
879  } else {
880  RAP = rcp(new RAPFactory());
881  ParameterList RAPparams;
882  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "transpose: use implicit", bool, RAPparams);
883 #if REUSE_MATRIX_GRAPHS
884  // This should be enable once Tpetra supports proper graph reuse in MMM
885  // At the moment, only Epetra does, and Tpetra throws
886  if (reuseType == "RP") {
887  RAPparams.set("Keep AP Pattern", true);
888  RAPparams.set("Keep RAP Pattern", true);
889  }
890 #endif
891  RAP->SetParameterList(RAPparams);
892  RAP->SetFactory("P", manager.GetFactory("P"));
893  if (!this->implicitTranspose_)
894  RAP->SetFactory("R", manager.GetFactory("R"));
895  if (MUELU_TEST_PARAM_2LIST(paramList, defaultList, "aggregation: export visualization data", bool, true)) {
896  RCP<AggregationExportFactory> aggExport = rcp(new AggregationExportFactory());
897  ParameterList aggExportParams;
898  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output filename", std::string, aggExportParams);
899  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: agg style", std::string, aggExportParams);
900  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: iter", int, aggExportParams);
901  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: time step", int, aggExportParams);
902  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: fine graph edges", bool, aggExportParams);
903  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: coarse graph edges", bool, aggExportParams);
904  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: output file: build colormap", bool, aggExportParams);
905  aggExport->SetParameterList(aggExportParams);
906  aggExport->SetFactory("DofsPerNode", manager.GetFactory("DofsPerNode"));
907  RAP->AddTransferFactory(aggExport);
908  }
909  manager.SetFactory("A", RAP);
910  }
911 
912  // === Coordinates ===
913  if (useCoordinates_) {
914  if (have_userCO) {
915  manager.SetFactory("Coordinates", NoFactory::getRCP());
916 
917  } else {
918  MUELU_KOKKOS_FACTORY(coords, CoordinatesTransferFactory, CoordinatesTransferFactory_kokkos);
919  coords->SetFactory("Aggregates", manager.GetFactory("Aggregates"));
920  coords->SetFactory("CoarseMap", manager.GetFactory("CoarseMap"));
921  manager.SetFactory("Coordinates", coords);
922 
923  if (paramList.isParameter("semicoarsen: number of levels")) {
924  RCP<ToggleCoordinatesTransferFactory> tf = rcp(new ToggleCoordinatesTransferFactory());
925  tf->SetFactory("Chosen P", manager.GetFactory("P"));
926  tf->AddCoordTransferFactory(semicoarsenFactory);
927  tf->AddCoordTransferFactory(coords);
928  manager.SetFactory("Coordinates", tf);
929 
930  }
931  RAP->AddTransferFactory(manager.GetFactory("Coordinates"));
932  }
933  }
934 
935  if (reuseType == "RP" || reuseType == "RAP" || reuseType == "full")
936  keeps.push_back(keep_pair("Nullspace", manager.GetFactory("Nullspace").get()));
937 
938  if (reuseType == "RP") {
939  keeps.push_back(keep_pair("P", manager.GetFactory("P").get()));
940  if (!this->implicitTranspose_)
941  keeps.push_back(keep_pair("R", manager.GetFactory("R").get()));
942  }
943  if ((reuseType == "tP" || reuseType == "RP" || reuseType == "emin") && useCoordinates_)
944  keeps.push_back(keep_pair("Coordinates", manager.GetFactory("Coordinates").get()));
945 
946  // === Repartitioning ===
947  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: enable", bool, enableRepart);
948  if (enableRepart) {
949 #ifdef HAVE_MPI
950  // Short summary of the issue: RebalanceTransferFactory shares ownership
951  // of "P" with SaPFactory, and therefore, changes the stored version.
952  // That means that if SaPFactory generated P, and stored it on the level,
953  // then after rebalancing the value in that storage changed. It goes
954  // against the concept of factories (I think), that every factory is
955  // responsible for its own objects, and they are immutable outside.
956  //
957  // In reuse, this is what happens: as we reuse Importer across setups,
958  // the order of factories changes, and coupled with shared ownership
959  // leads to problems.
960  // *First setup*
961  // SaP builds P [and stores it]
962  // TransP builds R [and stores it]
963  // RAP builds A [and stores it]
964  // RebalanceTransfer rebalances P [and changes the P stored by SaP] (*)
965  // RebalanceTransfer rebalances R
966  // RebalanceAc rebalances A
967  // *Second setup* ("RP" reuse)
968  // RebalanceTransfer rebalances P [which is incorrect due to (*)]
969  // RebalanceTransfer rebalances R
970  // RAP builds A [which is incorrect due to (*)]
971  // RebalanceAc rebalances A [which throws due to map inconsistency]
972  // ...
973  // *Second setup* ("tP" reuse)
974  // SaP builds P [and stores it]
975  // RebalanceTransfer rebalances P [and changes the P stored by SaP] (**)
976  // TransP builds R [which is incorrect due to (**)]
977  // RebalanceTransfer rebalances R
978  // ...
979  //
980  // Couple solutions to this:
981  // 1. [implemented] Requre "tP" and "PR" reuse to only be used with
982  // implicit rebalancing.
983  // 2. Do deep copy of P, and changed domain map and importer there.
984  // Need to investigate how expensive this is.
985  TEUCHOS_TEST_FOR_EXCEPTION(this->doPRrebalance_ && (reuseType == "tP" || reuseType == "RP"), Exceptions::InvalidArgument,
986  "Reuse types \"tP\" and \"PR\" require \"repartition: rebalance P and R\" set to \"false\"");
987 
988  TEUCHOS_TEST_FOR_EXCEPTION(aggType == "brick", Exceptions::InvalidArgument,
989  "Aggregation type \"brick\" requires \"repartition: enable\" set to \"false\"");
990 
991  MUELU_SET_VAR_2LIST(paramList, defaultList, "repartition: partitioner", std::string, partName);
992  TEUCHOS_TEST_FOR_EXCEPTION(partName != "zoltan" && partName != "zoltan2", Exceptions::InvalidArgument,
993  "Invalid partitioner name: \"" << partName << "\". Valid options: \"zoltan\", \"zoltan2\"");
994  // Partitioner
995  RCP<Factory> partitioner;
996  if (partName == "zoltan") {
997 #ifdef HAVE_MUELU_ZOLTAN
998  partitioner = rcp(new ZoltanInterface());
999  // NOTE: ZoltanInteface ("zoltan") does not support external parameters through ParameterList
1000 #else
1001  throw Exceptions::RuntimeError("Zoltan interface is not available");
1002 #endif
1003  } else if (partName == "zoltan2") {
1004 #ifdef HAVE_MUELU_ZOLTAN2
1005  partitioner = rcp(new Zoltan2Interface());
1006  ParameterList partParams;
1007  RCP<const ParameterList> partpartParams = rcp(new ParameterList(paramList.sublist("repartition: params", false)));
1008  partParams.set("ParameterList", partpartParams);
1009  partitioner->SetParameterList(partParams);
1010 #else
1011  throw Exceptions::RuntimeError("Zoltan2 interface is not available");
1012 #endif
1013  }
1014  partitioner->SetFactory("A", manager.GetFactory("A"));
1015  partitioner->SetFactory("Coordinates", manager.GetFactory("Coordinates"));
1016  manager.SetFactory("Partition", partitioner);
1017 
1018  // Repartitioner
1019  RCP<RepartitionFactory> repartFactory = rcp(new RepartitionFactory());
1020  ParameterList repartParams;
1021  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: start level", int, repartParams);
1022  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: min rows per proc", int, repartParams);
1023  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: max imbalance", double, repartParams);
1024  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: keep proc 0", bool, repartParams);
1025  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: print partition distribution", bool, repartParams);
1026  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap parts", bool, repartParams);
1027  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap num values", int, repartParams);
1028  repartFactory->SetParameterList(repartParams);
1029  repartFactory->SetFactory("A", manager.GetFactory("A"));
1030  repartFactory->SetFactory("Partition", manager.GetFactory("Partition"));
1031  manager.SetFactory("Importer", repartFactory);
1032  if (reuseType != "none")
1033  keeps.push_back(keep_pair("Importer", manager.GetFactory("Importer").get()));
1034 
1035  // Rebalanced A
1036  RCP<RebalanceAcFactory> newA = rcp(new RebalanceAcFactory());
1037  ParameterList rebAcParams;
1038  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, rebAcParams);
1039  newA->SetParameterList(rebAcParams);
1040  newA->SetFactory("A", manager.GetFactory("A"));
1041  newA->SetFactory("Importer", manager.GetFactory("Importer"));
1042  manager.SetFactory("A", newA);
1043 
1044  // Rebalanced P
1045  RCP<RebalanceTransferFactory> newP = rcp(new RebalanceTransferFactory());
1046  ParameterList newPparams;
1047  newPparams.set("type", "Interpolation");
1048  if (changedPRrebalance_)
1049  newPparams.set("repartition: rebalance P and R", this->doPRrebalance_);
1050  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, newPparams);
1051  newP-> SetParameterList(newPparams);
1052  newP-> SetFactory("Importer", manager.GetFactory("Importer"));
1053  newP-> SetFactory("P", manager.GetFactory("P"));
1054  if (!paramList.isParameter("semicoarsen: number of levels"))
1055  newP-> SetFactory("Nullspace", manager.GetFactory("Ptent"));
1056  else
1057  newP-> SetFactory("Nullspace", manager.GetFactory("P")); // TogglePFactory
1058  newP-> SetFactory("Coordinates", manager.GetFactory("Coordinates"));
1059  manager.SetFactory("P", newP);
1060  manager.SetFactory("Coordinates", newP);
1061 
1062  // Rebalanced R
1063  RCP<RebalanceTransferFactory> newR = rcp(new RebalanceTransferFactory());
1064  ParameterList newRparams;
1065  newRparams.set("type", "Restriction");
1066  MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: use subcommunicators", bool, newRparams);
1067  if (changedPRrebalance_)
1068  newRparams.set("repartition: rebalance P and R", this->doPRrebalance_);
1070  newRparams.set("transpose: use implicit", this->implicitTranspose_);
1071  newR-> SetParameterList(newRparams);
1072  newR-> SetFactory("Importer", manager.GetFactory("Importer"));
1073  if (!this->implicitTranspose_) {
1074  newR->SetFactory("R", manager.GetFactory("R"));
1075  manager.SetFactory("R", newR);
1076  }
1077 
1078  // NOTE: the role of NullspaceFactory is to provide nullspace on the finest
1079  // level if a user does not do that. For all other levels it simply passes
1080  // nullspace from a real factory to whoever needs it. If we don't use
1081  // repartitioning, that factory is "TentativePFactory"; if we do, it is
1082  // "RebalanceTransferFactory". But we still have to have NullspaceFactory as
1083  // the "Nullspace" of the manager
1084  nullSpace->SetFactory("Nullspace", newP);
1085 #else
1086  throw Exceptions::RuntimeError("No repartitioning available for a serial run");
1087 #endif
1088  }
1089  if (reuseType == "RAP" || reuseType == "full") {
1090  keeps.push_back(keep_pair("P", manager.GetFactory("P").get()));
1091  if (!this->implicitTranspose_)
1092  keeps.push_back(keep_pair("R", manager.GetFactory("R").get()));
1093  keeps.push_back(keep_pair("A", manager.GetFactory("A").get()));
1094  }
1095  }
1096 #undef MUELU_SET_VAR_2LIST
1097 #undef MUELU_TEST_AND_SET_VAR
1098 #undef MUELU_TEST_AND_SET_PARAM_2LIST
1099 #undef MUELU_TEST_PARAM_2LIST
1100 
1101  int LevenshteinDistance(const char* s, size_t len_s, const char* t, size_t len_t);
1102 
1103  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1105  ParameterList paramList = constParamList;
1106  const ParameterList& validList = *MasterList::List();
1107  // Validate up to maxLevels level specific parameter sublists
1108  const int maxLevels = 100;
1109 
1110  // Extract level specific list
1111  std::vector<ParameterList> paramLists;
1112  for (int levelID = 0; levelID < maxLevels; levelID++) {
1113  std::string sublistName = "level " + toString(levelID);
1114  if (paramList.isSublist(sublistName)) {
1115  paramLists.push_back(paramList.sublist(sublistName));
1116  // paramLists.back().setName(sublistName);
1117  paramList.remove(sublistName);
1118  }
1119  }
1120  paramLists.push_back(paramList);
1121  // paramLists.back().setName("main");
1122  //If Muemex is supported, hide custom level variables from validator by removing them from paramList's sublists
1123 #ifdef HAVE_MUELU_MATLAB
1124  for (size_t i = 0; i < paramLists.size(); i++) {
1125  std::vector<std::string> customVars; // list of names (keys) to be removed from list
1126 
1127  for(Teuchos::ParameterList::ConstIterator it = paramLists[i].begin(); it != paramLists[i].end(); it++) {
1128  std::string paramName = paramLists[i].name(it);
1129 
1130  if (IsParamMuemexVariable(paramName))
1131  customVars.push_back(paramName);
1132  }
1133 
1134  // Remove the keys
1135  for (size_t j = 0; j < customVars.size(); j++)
1136  paramLists[i].remove(customVars[j], false);
1137  }
1138 #endif
1139 
1140  const int maxDepth = 0;
1141  for (size_t i = 0; i < paramLists.size(); i++) {
1142  // validate every sublist
1143  try {
1144  paramLists[i].validateParameters(validList, maxDepth);
1145 
1146  } catch (const Teuchos::Exceptions::InvalidParameterName& e) {
1147  std::string eString = e.what();
1148 
1149  // Parse name from: <Error, the parameter {name="smoothe: type",...>
1150  size_t nameStart = eString.find_first_of('"') + 1;
1151  size_t nameEnd = eString.find_first_of('"', nameStart);
1152  std::string name = eString.substr(nameStart, nameEnd - nameStart);
1153 
1154  int bestScore = 100;
1155  std::string bestName = "";
1156  for (ParameterList::ConstIterator it = validList.begin(); it != validList.end(); it++) {
1157  const std::string& pName = validList.name(it);
1158  this->GetOStream(Runtime1) << "| " << pName;
1159  int score = LevenshteinDistance(name.c_str(), name.length(), pName.c_str(), pName.length());
1160  this->GetOStream(Runtime1) << " -> " << score << std::endl;
1161  if (score < bestScore) {
1162  bestScore = score;
1163  bestName = pName;
1164  }
1165  }
1166  if (bestScore < 10 && bestName != "") {
1167  TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameterName,
1168  eString << "The parameter name \"" + name + "\" is not valid. Did you mean \"" + bestName << "\"?\n");
1169 
1170  } else {
1171  TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameterName,
1172  eString << "The parameter name \"" + name + "\" is not valid.\n");
1173  }
1174  }
1175  }
1176  }
1177 
1178  // =====================================================================================================
1179  // ==================================== FACTORY interpreter ============================================
1180  // =====================================================================================================
1181  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1183  // Create a non const copy of the parameter list
1184  // Working with a modifiable list is much much easier than with original one
1185  ParameterList paramList = constParamList;
1186 
1187  // Parameter List Parsing:
1188  // ---------
1189  // <ParameterList name="MueLu">
1190  // <ParameterList name="Matrix">
1191  // </ParameterList>
1192  if (paramList.isSublist("Matrix")) {
1193  blockSize_ = paramList.sublist("Matrix").get<int>("number of equations", MasterList::getDefault<int>("number of equations"));
1194  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)
1195  }
1196 
1197  // create new FactoryFactory object if necessary
1198  if (factFact_ == Teuchos::null)
1199  factFact_ = Teuchos::rcp(new FactoryFactory());
1200 
1201  // Parameter List Parsing:
1202  // ---------
1203  // <ParameterList name="MueLu">
1204  // <ParameterList name="Factories"> <== call BuildFactoryMap() on this parameter list
1205  // ...
1206  // </ParameterList>
1207  // </ParameterList>
1208  FactoryMap factoryMap;
1209  FactoryManagerMap factoryManagers;
1210  if (paramList.isSublist("Factories"))
1211  this->BuildFactoryMap(paramList.sublist("Factories"), factoryMap, factoryMap, factoryManagers);
1212 
1213  // Parameter List Parsing:
1214  // ---------
1215  // <ParameterList name="MueLu">
1216  // <ParameterList name="Hierarchy">
1217  // <Parameter name="verbose" type="string" value="Warnings"/> <== get
1218  // <Parameter name="numDesiredLevel" type="int" value="10"/> <== get
1219  //
1220  // <ParameterList name="firstLevel"> <== parse first args and call BuildFactoryMap() on the rest of this parameter list
1221  // ...
1222  // </ParameterList>
1223  // </ParameterList>
1224  // </ParameterList>
1225  if (paramList.isSublist("Hierarchy")) {
1226  ParameterList hieraList = paramList.sublist("Hierarchy"); // copy because list temporally modified (remove 'id')
1227 
1228  // Get hierarchy options
1229  if (hieraList.isParameter("max levels")) {
1230  this->numDesiredLevel_ = hieraList.get<int>("max levels");
1231  hieraList.remove("max levels");
1232  }
1233 
1234  if (hieraList.isParameter("coarse: max size")) {
1235  this->maxCoarseSize_ = hieraList.get<int>("coarse: max size");
1236  hieraList.remove("coarse: max size");
1237  }
1238 
1239  if (hieraList.isParameter("repartition: rebalance P and R")) {
1240  this->doPRrebalance_ = hieraList.get<bool>("repartition: rebalance P and R");
1241  hieraList.remove("repartition: rebalance P and R");
1242  }
1243 
1244  if (hieraList.isParameter("transpose: use implicit")) {
1245  this->implicitTranspose_ = hieraList.get<bool>("transpose: use implicit");
1246  hieraList.remove("transpose: use implicit");
1247  }
1248 
1249  //TODO Move this its own class or MueLu::Utils?
1250  std::map<std::string,MsgType> verbMap;
1251  //for developers
1252  verbMap["Errors"] = Errors;
1253  verbMap["Warnings0"] = Warnings0;
1254  verbMap["Warnings00"] = Warnings00;
1255  verbMap["Warnings1"] = Warnings1;
1256  verbMap["PerfWarnings"] = PerfWarnings;
1257  verbMap["Runtime0"] = Runtime0;
1258  verbMap["Runtime1"] = Runtime1;
1259  verbMap["RuntimeTimings"] = RuntimeTimings;
1260  verbMap["NoTimeReport"] = NoTimeReport;
1261  verbMap["Parameters0"] = Parameters0;
1262  verbMap["Parameters1"] = Parameters1;
1263  verbMap["Statistics0"] = Statistics0;
1264  verbMap["Statistics1"] = Statistics1;
1265  verbMap["Timings0"] = Timings0;
1266  verbMap["Timings1"] = Timings1;
1267  verbMap["TimingsByLevel"] = TimingsByLevel;
1268  verbMap["External"] = External;
1269  verbMap["Debug"] = Debug;
1270  verbMap["Test"] = Test;
1271  //for users and developers
1272  verbMap["None"] = None;
1273  verbMap["Low"] = Low;
1274  verbMap["Medium"] = Medium;
1275  verbMap["High"] = High;
1276  verbMap["Extreme"] = Extreme;
1277  if (hieraList.isParameter("verbosity")) {
1278  std::string vl = hieraList.get<std::string>("verbosity");
1279  hieraList.remove("verbosity");
1280  //TODO Move this to its own class or MueLu::Utils?
1281  if (verbMap.find(vl) != verbMap.end())
1282  this->verbosity_ = verbMap[vl];
1283  else
1284  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::ParameterListInterpreter():: invalid verbosity level");
1285  }
1286 
1287  if (hieraList.isParameter("dependencyOutputLevel"))
1288  this->graphOutputLevel_ = hieraList.get<int>("dependencyOutputLevel");
1289 
1290  // Check for the reuse case
1291  if (hieraList.isParameter("reuse"))
1293 
1294  if (hieraList.isSublist("DataToWrite")) {
1295  //TODO We should be able to specify any data. If it exists, write it.
1296  //TODO This would requires something like std::set<dataName,Array<int> >
1297  ParameterList foo = hieraList.sublist("DataToWrite");
1298  std::string dataName = "Matrices";
1299  if (foo.isParameter(dataName))
1300  this->matricesToPrint_ = Teuchos::getArrayFromStringParameter<int>(foo,dataName);
1301  dataName = "Prolongators";
1302  if (foo.isParameter(dataName))
1303  this->prolongatorsToPrint_ = Teuchos::getArrayFromStringParameter<int>(foo,dataName);
1304  dataName = "Restrictors";
1305  if (foo.isParameter(dataName))
1306  this->restrictorsToPrint_ = Teuchos::getArrayFromStringParameter<int>(foo,dataName);
1307  }
1308 
1309  // Get level configuration
1310  for (ParameterList::ConstIterator param = hieraList.begin(); param != hieraList.end(); ++param) {
1311  const std::string & paramName = hieraList.name(param);
1312 
1313  if (paramName != "DataToWrite" && hieraList.isSublist(paramName)) {
1314  ParameterList levelList = hieraList.sublist(paramName); // copy because list temporally modified (remove 'id')
1315 
1316  int startLevel = 0; if(levelList.isParameter("startLevel")) { startLevel = levelList.get<int>("startLevel"); levelList.remove("startLevel"); }
1317  int numDesiredLevel = 1; if(levelList.isParameter("numDesiredLevel")) { numDesiredLevel = levelList.get<int>("numDesiredLevel"); levelList.remove("numDesiredLevel"); }
1318 
1319  // Parameter List Parsing:
1320  // ---------
1321  // <ParameterList name="firstLevel">
1322  // <Parameter name="startLevel" type="int" value="0"/>
1323  // <Parameter name="numDesiredLevel" type="int" value="1"/>
1324  // <Parameter name="verbose" type="string" value="Warnings"/>
1325  //
1326  // [] <== call BuildFactoryMap() on the rest of the parameter list
1327  //
1328  // </ParameterList>
1329  FactoryMap levelFactoryMap;
1330  BuildFactoryMap(levelList, factoryMap, levelFactoryMap, factoryManagers);
1331 
1332  RCP<FactoryManagerBase> m = rcp(new FactoryManager(levelFactoryMap));
1333 
1334  if (startLevel >= 0)
1335  this->AddFactoryManager(startLevel, numDesiredLevel, m);
1336  else
1337  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::ParameterListInterpreter():: invalid level id");
1338  } /* TODO: else { } */
1339  }
1340  }
1341  }
1342 
1343  // Parameter List Parsing:
1344  // Create an entry in factoryMap for each parameter of the list paramList
1345  // ---------
1346  // <ParameterList name="...">
1347  // <Parameter name="smootherFact0" type="string" value="TrilinosSmoother"/>
1348  //
1349  // <ParameterList name="smootherFact1">
1350  // <Parameter name="type" type="string" value="TrilinosSmoother"/>
1351  // ...
1352  // </ParameterList>
1353  // </ParameterList>
1354  //
1355  //TODO: static?
1356  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1358  BuildFactoryMap(const ParameterList& paramList, const FactoryMap& factoryMapIn, FactoryMap& factoryMapOut, FactoryManagerMap& factoryManagers) const {
1359  for (ParameterList::ConstIterator param = paramList.begin(); param != paramList.end(); ++param) {
1360  const std::string & paramName = paramList.name(param);
1361  const Teuchos::ParameterEntry & paramValue = paramList.entry(param);
1362 
1363  //TODO: do not allow name of existing MueLu classes (can be tested using FactoryFactory)
1364 
1365  // TODO: add support for "factory groups" which are stored in a map.
1366  // A factory group has a name and a list of factories
1367 
1368  if (paramValue.isList()) {
1369  ParameterList paramList1 = Teuchos::getValue<ParameterList>(paramValue);
1370  if (paramList1.isParameter("factory")) { // default: just a factory definition
1371  factoryMapOut[paramName] = factFact_->BuildFactory(paramValue, factoryMapIn, factoryManagers);
1372 
1373  } else if (paramList1.isParameter("group")) { // definitiion of a factory group (for a factory manager)
1374  std::string groupType = paramList1.get<std::string>("group");
1375  TEUCHOS_TEST_FOR_EXCEPTION(groupType!="FactoryManager", Exceptions::RuntimeError, "group must be of type \"FactoryManager\".");
1376 
1377  ParameterList groupList = paramList1; // copy because list temporally modified (remove 'id')
1378  groupList.remove("group");
1379 
1380  FactoryMap groupFactoryMap;
1381  BuildFactoryMap(groupList, factoryMapIn, groupFactoryMap, factoryManagers);
1382 
1383  // do not store groupFactoryMap in factoryMapOut
1384  // Create a factory manager object from groupFactoryMap
1385  RCP<FactoryManagerBase> m = rcp(new FactoryManager(groupFactoryMap));
1386 
1387  factoryManagers[paramName] = m;
1388 
1389  } else {
1390  this->GetOStream(Warnings0) << "Could not interpret parameter list " << paramList1 << std::endl;
1391  TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "XML Parameter list must either be of type \"factory\" or of type \"group\".");
1392  }
1393  } else {
1394  // default: just a factory (no parameter list)
1395  factoryMapOut[paramName] = factFact_->BuildFactory(paramValue, factoryMapIn, factoryManagers);
1396  }
1397  }
1398  }
1399 
1400  // =====================================================================================================
1401  // ======================================= MISC functions ==============================================
1402  // =====================================================================================================
1403  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1405  try {
1406  Matrix& A = dynamic_cast<Matrix&>(Op);
1407  if (A.GetFixedBlockSize() != blockSize_)
1408  this->GetOStream(Warnings0) << "Setting matrix block size to " << blockSize_ << " (value of the parameter in the list) "
1409  << "instead of " << A.GetFixedBlockSize() << " (provided matrix)." << std::endl;
1410 
1411  A.SetFixedBlockSize(blockSize_, dofOffset_);
1412 
1413  } catch (std::bad_cast& e) {
1414  this->GetOStream(Warnings0) << "Skipping setting block size as the operator is not a matrix" << std::endl;
1415  }
1416  }
1417 
1418  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1421  H.SetCycle(Cycle_);
1422  }
1423 
1424  static bool compare(const ParameterList& list1, const ParameterList& list2) {
1425  // First loop through and validate the parameters at this level.
1426  // In addition, we generate a list of sublists that we will search next
1427  for (ParameterList::ConstIterator it = list1.begin(); it != list1.end(); it++) {
1428  const std::string& name = it->first;
1429  const Teuchos::ParameterEntry& entry1 = it->second;
1430 
1431  const Teuchos::ParameterEntry *entry2 = list2.getEntryPtr(name);
1432  if (!entry2) // entry is not present in the second list
1433  return false;
1434  if (entry1.isList() && entry2->isList()) { // sublist check
1435  compare(Teuchos::getValue<ParameterList>(entry1), Teuchos::getValue<ParameterList>(*entry2));
1436  continue;
1437  }
1438  if (entry1.getAny(false) != entry2->getAny(false)) // entries have different types or different values
1439  return false;
1440  }
1441 
1442  return true;
1443  }
1444 
1445  static inline bool areSame(const ParameterList& list1, const ParameterList& list2) {
1446  return compare(list1, list2) && compare(list2, list1);
1447  }
1448 
1449 } // namespace MueLu
1450 
1451 #define MUELU_PARAMETERLISTINTERPRETER_SHORT
1452 #endif /* MUELU_PARAMETERLISTINTERPRETER_DEF_HPP */
Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
This class specifies the default factory that should generate some data on a Level if the data does n...
Teuchos::Array< int > matricesToPrint_
void SetupHierarchy(Hierarchy &H) const
Call the SetupHierarchy routine from the HiearchyManager object.
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
#define MUELU_KOKKOS_FACTORY_NO_DECL(varName, oldFactory, newFactory)
Factory for generating coarse level map. Used by TentativePFactory.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
void SetVerbLevel(const VerbLevel verbLevel)
Set the verbosity level of this object.
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.
int LevenshteinDistance(const char *s, size_t len_s, const char *t, size_t len_t)
Print external lib objects.
Class that encapsulates external library smoothers.
static void DisableMultipleCheckGlobally()
#define MUELU_TEST_AND_SET_VAR(paramList, paramName, paramType, varName)
Factory for building permutation matrix that can be be used to shuffle data (matrices, vectors) among processes.
Print more statistics.
Print additional debugging information.
void AddFactoryManager(int startLevel, int numDesiredLevel, RCP< FactoryManagerBase > manager)
One-liner description of what is happening.
void SetCycle(CycleType Cycle)
Supports VCYCLE and WCYCLE types.
Namespace for MueLu classes and methods.
#define MUELU_KOKKOS_FACTORY(varName, oldFactory, newFactory)
Interface to Zoltan library.
bool IsParamMuemexVariable(const std::string &name)
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Get factory associated with a particular data name.
std::map< std::string, RCP< FactoryManagerBase > > FactoryManagerMap
std::map< std::string, RCP< const FactoryBase > > FactoryMap
Print skeleton for the run, i.e. factory calls and used parameters.
Factory for building tentative prolongator.
#define TEST_MUTUALLY_EXCLUSIVE(arg1, arg2)
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
Factory for coarsening a graph with uncoupled aggregation.
virtual void SetupOperator(Operator &A) const
Setup Operator object.
Prolongator factory performing semi-coarsening.
Factory for building restriction operators using a prolongator factory.
Additional warnings.
int blockSize_
block size of matrix (fixed block size)
Important warning messages (more verbose)
Print statistics that do not involve significant additional computation.
Teuchos::RCP< FactoryFactory > factFact_
Internal factory for factories.
#define MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, paramName, paramType, listWrite)
static bool compare(const ParameterList &list1, const ParameterList &list2)
static CycleType GetDefaultCycle()
Detailed timing information (use Teuchos::TimeMonitor::summarize() to print)
Factory for interacting with Matlab.
Factory for interacting with Matlab.
Factory for building line detection information.
void SetFactoryParameterList(const Teuchos::ParameterList &paramList)
Factory interpreter stuff.
std::pair< std::string, const FactoryBase * > keep_pair
Factory to export aggregation info or visualize aggregates using VTK.
Interface to Zoltan2 library.
Prolongator factory which allows switching between two different prolongator strategies.
By default, enabled timers appears in the teuchos time monitor summary. Use this option if you do not...
RCP< FactoryManagerBase > GetFactoryManager(int levelID) const
Factory for building the constraint operator.
ParameterListInterpreter(Teuchos::ParameterList &paramList, Teuchos::RCP< const Teuchos::Comm< int > > comm=Teuchos::null, Teuchos::RCP< FactoryFactory > factFact=Teuchos::null)
Constructor that accepts a user-provided ParameterList.
Applies permutation to grid transfer operators.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Set Factory.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Print class parameters.
Timers that are enabled (using Timings0/Timings1) will be printed during the execution.
GlobalOrdinal dofOffset_
global offset variable describing offset of DOFs in operator
Performance warnings.
Record timing information level by level. Must be used in combinaison with Timings0/Timings1.
Factory for creating a graph base on a given matrix.
Teuchos::Array< int > prolongatorsToPrint_
void BuildFactoryMap(const Teuchos::ParameterList &paramList, const FactoryMap &factoryMapIn, FactoryMap &factoryMapOut, FactoryManagerMap &factoryManagers) const
Class that encapsulates Matlab smoothers.
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameter list for Parameter list interpreter.
Class for transferring coordinates from a finer level to a coarser one.
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.
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.
std::map< int, std::vector< keep_pair > > keep_
#define TEST_MUTUALLY_EXCLUSIVE_S(arg1, arg2)
void Validate(const Teuchos::ParameterList &paramList) const
#define MUELU_TEST_PARAM_2LIST(paramList, defaultList, paramName, paramType, cmpValue)
Print class parameters (more parameters, more verbose)
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
Description of what is happening (more verbose)
Factory for building coarse matrices.
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Factory for building Smoothed Aggregation prolongators.Input/output of SaPFactory
Xpetra::global_size_t maxCoarseSize_
Factory for building uncoupled aggregates.
static Teuchos::RCP< const Teuchos::ParameterList > List()
Return a "master" list of all valid parameters and their default values.
#define MUELU_SET_VAR_2LIST(paramList, defaultList, paramName, paramType, varName)
Teuchos::Array< int > restrictorsToPrint_
void SetEasyParameterList(const Teuchos::ParameterList &paramList)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Factory for generating nullspace.
long ExtractNonSerializableData(const Teuchos::ParameterList &inList, Teuchos::ParameterList &serialList, Teuchos::ParameterList &nonSerialList)
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
Exception throws to report invalid user entry.
static const RCP< const NoFactory > getRCP()
Static Get() functions.
CycleType Cycle_
multigrid cycle type (V-cycle or W-cycle)