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