MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_MLParameterListInterpreter_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_MLPARAMETERLISTINTERPRETER_DEF_HPP
47 #define MUELU_MLPARAMETERLISTINTERPRETER_DEF_HPP
48 
50 
51 #include "MueLu_ConfigDefs.hpp"
52 #if defined(HAVE_MUELU_ML)
53 #include <ml_ValidateParameters.h>
54 #endif
55 
56 #include <Xpetra_Matrix.hpp>
57 #include <Xpetra_MatrixUtils.hpp>
58 #include <Xpetra_MultiVector.hpp>
59 #include <Xpetra_MultiVectorFactory.hpp>
60 #include <Xpetra_Operator.hpp>
61 
63 
64 #include "MueLu_Level.hpp"
65 #include "MueLu_Hierarchy.hpp"
66 #include "MueLu_FactoryManager.hpp"
67 
68 #include "MueLu_TentativePFactory.hpp"
69 #include "MueLu_SaPFactory.hpp"
70 #include "MueLu_PgPFactory.hpp"
71 #include "MueLu_AmalgamationFactory.hpp"
72 #include "MueLu_TransPFactory.hpp"
73 #include "MueLu_GenericRFactory.hpp"
74 #include "MueLu_SmootherPrototype.hpp"
75 #include "MueLu_SmootherFactory.hpp"
76 #include "MueLu_TrilinosSmoother.hpp"
77 #include "MueLu_IfpackSmoother.hpp"
78 #include "MueLu_DirectSolver.hpp"
79 #include "MueLu_HierarchyUtils.hpp"
80 #include "MueLu_RAPFactory.hpp"
81 #include "MueLu_CoalesceDropFactory.hpp"
82 #include "MueLu_UncoupledAggregationFactory.hpp"
83 #include "MueLu_NullspaceFactory.hpp"
85 
86 #include "MueLu_CoalesceDropFactory_kokkos.hpp"
87 // #include "MueLu_CoordinatesTransferFactory_kokkos.hpp"
88 #include "MueLu_TentativePFactory_kokkos.hpp"
89 
90 #if defined(HAVE_MUELU_ISORROPIA) && defined(HAVE_MPI)
91 #include "MueLu_IsorropiaInterface.hpp"
92 #include "MueLu_RepartitionHeuristicFactory.hpp"
93 #include "MueLu_RepartitionFactory.hpp"
94 #include "MueLu_RebalanceTransferFactory.hpp"
95 #include "MueLu_RepartitionInterface.hpp"
96 #include "MueLu_RebalanceAcFactory.hpp"
97 //#include "MueLu_RebalanceMapFactory.hpp"
98 #endif
99 
100 // Note: do not add options that are only recognized by MueLu.
101 
102 // TODO: this parameter list interpreter should force MueLu to use default ML parameters
103 // - Ex: smoother sweep=2 by default for ML
104 
105 // Read a parameter value from a parameter list and store it into a variable named 'varName'
106 #define MUELU_READ_PARAM(paramList, paramStr, varType, defaultValue, varName) \
107  varType varName = defaultValue; \
108  if (paramList.isParameter(paramStr)) varName = paramList.get<varType>(paramStr);
109 
110 // Read a parameter value from a paraeter list and copy it into a new parameter list (with another parameter name)
111 #define MUELU_COPY_PARAM(paramList, paramStr, varType, defaultValue, outParamList, outParamStr) \
112  if (paramList.isParameter(paramStr)) \
113  outParamList.set(outParamStr, paramList.get<varType>(paramStr)); \
114  else \
115  outParamList.set(outParamStr, static_cast<varType>(defaultValue));
116 
117 namespace MueLu {
118 
119 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121  : nullspace_(NULL)
122  , xcoord_(NULL)
123  , ycoord_(NULL)
124  , zcoord_(NULL)
125  , TransferFacts_(factoryList)
126  , blksize_(1) {
127  if (paramList.isParameter("xml parameter file")) {
128  std::string filename = paramList.get("xml parameter file", "");
129  if (filename.length() != 0) {
130  TEUCHOS_TEST_FOR_EXCEPTION(comm.is_null(), Exceptions::RuntimeError, "xml parameter file requires a valid comm");
131  Teuchos::ParameterList paramList2 = paramList;
132  Teuchos::updateParametersFromXmlFileAndBroadcast(filename, Teuchos::Ptr<Teuchos::ParameterList>(&paramList2), *comm);
133  paramList2.remove("xml parameter file");
134  SetParameterList(paramList2);
135  } else
136  SetParameterList(paramList);
137  } else
138  SetParameterList(paramList);
139 }
140 
141 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
143  : nullspace_(NULL)
144  , TransferFacts_(factoryList)
145  , blksize_(1) {
146  Teuchos::RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile(xmlFileName);
147  SetParameterList(*paramList);
148 }
149 
150 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152  Teuchos::ParameterList paramList = paramList_in;
153 
154  //
155  // Read top-level of the parameter list
156  //
157 
158  // hard-coded default values == ML defaults according to the manual
159  MUELU_READ_PARAM(paramList, "ML output", int, 0, verbosityLevel);
160  MUELU_READ_PARAM(paramList, "max levels", int, 10, maxLevels);
161  MUELU_READ_PARAM(paramList, "PDE equations", int, 1, nDofsPerNode);
162 
163  MUELU_READ_PARAM(paramList, "coarse: max size", int, 128, maxCoarseSize);
164 
165  MUELU_READ_PARAM(paramList, "aggregation: type", std::string, "Uncoupled", agg_type);
166  // MUELU_READ_PARAM(paramList, "aggregation: threshold", double, 0.0, agg_threshold);
167  MUELU_READ_PARAM(paramList, "aggregation: damping factor", double, (double)4 / (double)3, agg_damping);
168  // MUELU_READ_PARAM(paramList, "aggregation: smoothing sweeps", int, 1, agg_smoothingsweeps);
169  MUELU_READ_PARAM(paramList, "aggregation: nodes per aggregate", int, 1, minPerAgg);
170  MUELU_READ_PARAM(paramList, "aggregation: keep Dirichlet bcs", bool, false, bKeepDirichletBcs); // This is a MueLu specific extension that does not exist in ML
171  MUELU_READ_PARAM(paramList, "aggregation: max neighbours already aggregated", int, 0, maxNbrAlreadySelected); // This is a MueLu specific extension that does not exist in M
172  MUELU_READ_PARAM(paramList, "aggregation: aux: enable", bool, false, agg_use_aux);
173  MUELU_READ_PARAM(paramList, "aggregation: aux: threshold", double, false, agg_aux_thresh);
174 
175  MUELU_READ_PARAM(paramList, "null space: type", std::string, "default vectors", nullspaceType);
176  MUELU_READ_PARAM(paramList, "null space: dimension", int, -1, nullspaceDim); // TODO: ML default not in documentation
177  MUELU_READ_PARAM(paramList, "null space: vectors", double*, NULL, nullspaceVec); // TODO: ML default not in documentation
178 
179  MUELU_READ_PARAM(paramList, "energy minimization: enable", bool, false, bEnergyMinimization);
180 
181  MUELU_READ_PARAM(paramList, "RAP: fix diagonal", bool, false, bFixDiagonal); // This is a MueLu specific extension that does not exist in ML
182 
183  MUELU_READ_PARAM(paramList, "x-coordinates", double*, NULL, xcoord);
184  MUELU_READ_PARAM(paramList, "y-coordinates", double*, NULL, ycoord);
185  MUELU_READ_PARAM(paramList, "z-coordinates", double*, NULL, zcoord);
186 
187  //
188  // Move smoothers/aggregation/coarse parameters to sublists
189  //
190 
191  // ML allows to have level-specific smoothers/aggregation/coarse parameters at the top level of the list or/and defined in sublists:
192  // See also: ML Guide section 6.4.1, MueLu::CreateSublists, ML_CreateSublists
193  ParameterList paramListWithSubList;
194  MueLu::CreateSublists(paramList, paramListWithSubList);
195  paramList = paramListWithSubList; // swap
196 
197  // pull out "use kokkos refactor"
198  bool setKokkosRefactor = false;
199  bool useKokkosRefactor = !Node::is_serial;
200  if (paramList.isType<bool>("use kokkos refactor")) {
201  useKokkosRefactor = paramList.get<bool>("use kokkos refactor");
202  setKokkosRefactor = true;
203  paramList.remove("use kokkos refactor");
204  }
205 
206  //
207  // Validate parameter list
208  //
209 
210  {
211  bool validate = paramList.get("ML validate parameter list", true); /* true = default in ML */
212  if (validate) {
213 #if defined(HAVE_MUELU_ML) && defined(HAVE_MUELU_EPETRA)
214  // Validate parameter list using ML validator
215  int depth = paramList.get("ML validate depth", 5); /* 5 = default in ML */
216  TEUCHOS_TEST_FOR_EXCEPTION(!ML_Epetra::ValidateMLPParameters(paramList, depth), Exceptions::RuntimeError,
217  "ERROR: ML's Teuchos::ParameterList contains incorrect parameter!");
218 #else
219  // If no validator available: issue a warning and set parameter value to false in the output list
220  this->GetOStream(Warnings0) << "Warning: MueLu_ENABLE_ML=OFF. The parameter list cannot be validated." << std::endl;
221  paramList.set("ML validate parameter list", false);
222 
223 #endif // HAVE_MUELU_ML
224  } // if(validate)
225  } // scope
226 
227  // Matrix option
228  blksize_ = nDofsPerNode;
229 
230  // Translate verbosity parameter
231 
232  // Translate verbosity parameter
233  MsgType eVerbLevel = None;
234  if (verbosityLevel == 0) eVerbLevel = None;
235  if (verbosityLevel >= 1) eVerbLevel = Low;
236  if (verbosityLevel >= 5) eVerbLevel = Medium;
237  if (verbosityLevel >= 10) eVerbLevel = High;
238  if (verbosityLevel >= 11) eVerbLevel = Extreme;
239  if (verbosityLevel >= 42) eVerbLevel = Test;
240  if (verbosityLevel >= 43) eVerbLevel = InterfaceTest;
241  this->verbosity_ = eVerbLevel;
242 
243  TEUCHOS_TEST_FOR_EXCEPTION(agg_type != "Uncoupled", Exceptions::RuntimeError,
244  "MueLu::MLParameterListInterpreter::SetParameterList(): parameter \"aggregation: type\": only 'Uncoupled' aggregation is supported.");
245 
246  // Create MueLu factories
247  RCP<Factory> dropFact;
248  if (useKokkosRefactor)
249  dropFact = rcp(new CoalesceDropFactory_kokkos());
250  else
251  dropFact = rcp(new CoalesceDropFactory());
252 
253  if (agg_use_aux) {
254  dropFact->SetParameter("aggregation: drop scheme", Teuchos::ParameterEntry(std::string("distance laplacian")));
255  dropFact->SetParameter("aggregation: drop tol", Teuchos::ParameterEntry(agg_aux_thresh));
256  }
257 
258  // Uncoupled aggregation
259  RCP<Factory> AggFact = Teuchos::null;
260  AggFact = rcp(new UncoupledAggregationFactory());
261 
262  AggFact->SetFactory("Graph", dropFact);
263  AggFact->SetFactory("DofsPerNode", dropFact);
264  AggFact->SetParameter("aggregation: preserve Dirichlet points", Teuchos::ParameterEntry(bKeepDirichletBcs));
265  AggFact->SetParameter("aggregation: ordering", Teuchos::ParameterEntry(std::string("natural")));
266  AggFact->SetParameter("aggregation: max selected neighbors", Teuchos::ParameterEntry(maxNbrAlreadySelected));
267  AggFact->SetParameter("aggregation: min agg size", Teuchos::ParameterEntry(minPerAgg));
268 
269  if (verbosityLevel > 3) {
270  std::ostringstream oss;
271  oss << "========================= Aggregate option summary  =========================" << std::endl;
272  oss << "min Nodes per aggregate :              " << minPerAgg << std::endl;
273  oss << "min # of root nbrs already aggregated : " << maxNbrAlreadySelected << std::endl;
274  oss << "aggregate ordering :                    natural" << std::endl;
275  oss << "=============================================================================" << std::endl;
276  this->GetOStream(Runtime1) << oss.str();
277  }
278 
279  RCP<Factory> PFact;
280  RCP<Factory> RFact;
281  RCP<Factory> PtentFact;
282  if (useKokkosRefactor)
283  PtentFact = rcp(new TentativePFactory_kokkos());
284  else
285  PtentFact = rcp(new TentativePFactory());
286  if (agg_damping == 0.0 && bEnergyMinimization == false) {
287  // tentative prolongation operator (PA-AMG)
288  PFact = PtentFact;
289  RFact = rcp(new TransPFactory());
290  } else if (agg_damping != 0.0 && bEnergyMinimization == false) {
291  // smoothed aggregation (SA-AMG)
292  RCP<Factory> SaPFact = rcp(new SaPFactory());
293  SaPFact->SetParameter("sa: damping factor", ParameterEntry(agg_damping));
294  PFact = SaPFact;
295  RFact = rcp(new TransPFactory());
296  } else if (bEnergyMinimization == true) {
297  // Petrov Galerkin PG-AMG smoothed aggregation (energy minimization in ML)
298  PFact = rcp(new PgPFactory());
299  RFact = rcp(new GenericRFactory());
300  }
301 
302  RCP<RAPFactory> AcFact = rcp(new RAPFactory());
303  AcFact->SetParameter("RepairMainDiagonal", Teuchos::ParameterEntry(bFixDiagonal));
304  for (size_t i = 0; i < TransferFacts_.size(); i++) {
305  AcFact->AddTransferFactory(TransferFacts_[i]);
306  }
307 
308  //
309  // introduce rebalancing
310  //
311 #if defined(HAVE_MUELU_ISORROPIA) && defined(HAVE_MPI)
312  Teuchos::RCP<Factory> RebalancedPFact = Teuchos::null;
313  Teuchos::RCP<Factory> RebalancedRFact = Teuchos::null;
314  Teuchos::RCP<Factory> RepartitionFact = Teuchos::null;
315  Teuchos::RCP<RebalanceAcFactory> RebalancedAFact = Teuchos::null;
316 
317  MUELU_READ_PARAM(paramList, "repartition: enable", int, 0, bDoRepartition);
318  if (bDoRepartition == 1) {
319  // The Factory Manager will be configured to return the rebalanced versions of P, R, A by default.
320  // Everytime we want to use the non-rebalanced versions, we need to explicitly define the generating factory.
321  RFact->SetFactory("P", PFact);
322  //
323  AcFact->SetFactory("P", PFact);
324  AcFact->SetFactory("R", RFact);
325 
326  // define rebalancing factory for coarse matrix
328  rebAmalgFact->SetFactory("A", AcFact);
329 
330  MUELU_READ_PARAM(paramList, "repartition: max min ratio", double, 1.3, maxminratio);
331  MUELU_READ_PARAM(paramList, "repartition: min per proc", int, 512, minperproc);
332 
333  // Repartitioning heuristic
335  {
336  Teuchos::ParameterList paramListRepFact;
337  paramListRepFact.set("repartition: min rows per proc", minperproc);
338  paramListRepFact.set("repartition: max imbalance", maxminratio);
339  RepartitionHeuristicFact->SetParameterList(paramListRepFact);
340  }
341  RepartitionHeuristicFact->SetFactory("A", AcFact);
342 
343  // create "Partition"
345  isoInterface->SetFactory("A", AcFact);
346  isoInterface->SetFactory("number of partitions", RepartitionHeuristicFact);
347  isoInterface->SetFactory("UnAmalgamationInfo", rebAmalgFact);
348 
349  // create "Partition" by unamalgamtion
351  repInterface->SetFactory("A", AcFact);
352  repInterface->SetFactory("number of partitions", RepartitionHeuristicFact);
353  repInterface->SetFactory("AmalgamatedPartition", isoInterface);
354  // repInterface->SetFactory("UnAmalgamationInfo", rebAmalgFact); // not necessary?
355 
356  // Repartitioning (creates "Importer" from "Partition")
357  RepartitionFact = Teuchos::rcp(new RepartitionFactory());
358  RepartitionFact->SetFactory("A", AcFact);
359  RepartitionFact->SetFactory("number of partitions", RepartitionHeuristicFact);
360  RepartitionFact->SetFactory("Partition", repInterface);
361 
362  // Reordering of the transfer operators
363  RebalancedPFact = Teuchos::rcp(new RebalanceTransferFactory());
364  RebalancedPFact->SetParameter("type", Teuchos::ParameterEntry(std::string("Interpolation")));
365  RebalancedPFact->SetFactory("P", PFact);
366  RebalancedPFact->SetFactory("Nullspace", PtentFact);
367  RebalancedPFact->SetFactory("Importer", RepartitionFact);
368 
369  RebalancedRFact = Teuchos::rcp(new RebalanceTransferFactory());
370  RebalancedRFact->SetParameter("type", Teuchos::ParameterEntry(std::string("Restriction")));
371  RebalancedRFact->SetFactory("R", RFact);
372  RebalancedRFact->SetFactory("Importer", RepartitionFact);
373 
374  // Compute Ac from rebalanced P and R
375  RebalancedAFact = Teuchos::rcp(new RebalanceAcFactory());
376  RebalancedAFact->SetFactory("A", AcFact);
377  }
378 #else // #ifdef HAVE_MUELU_ISORROPIA
379  // Get rid of [-Wunused] warnings
380  //(void)
381  //
382  // ^^^ FIXME (mfh 17 Nov 2013) That definitely doesn't compile.
383 #endif
384 
385  //
386  // Nullspace factory
387  //
388 
389  // Set fine level nullspace
390  // extract pre-computed nullspace from ML parameter list
391  // store it in nullspace_ and nullspaceDim_
392  if (nullspaceType != "default vectors") {
393  TEUCHOS_TEST_FOR_EXCEPTION(nullspaceType != "pre-computed", Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: no valid nullspace (no pre-computed null space). error.");
394  TEUCHOS_TEST_FOR_EXCEPTION(nullspaceDim == -1, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: no valid nullspace (nullspace dim == -1). error.");
395  TEUCHOS_TEST_FOR_EXCEPTION(nullspaceVec == NULL, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: no valid nullspace (nullspace == NULL). You have to provide a valid fine-level nullspace in \'null space: vectors\'");
396 
397  nullspaceDim_ = nullspaceDim;
398  nullspace_ = nullspaceVec;
399  }
400 
402  nspFact->SetFactory("Nullspace", PtentFact);
403 
404  // Stash coordinates
405  xcoord_ = xcoord;
406  ycoord_ = ycoord;
407  zcoord_ = zcoord;
408 
409  //
410  // Hierarchy + FactoryManager
411  //
412 
413  // Hierarchy options
414  this->numDesiredLevel_ = maxLevels;
415  this->maxCoarseSize_ = maxCoarseSize;
416 
417  //
418  // Coarse Smoother
419  //
420  ParameterList& coarseList = paramList.sublist("coarse: list");
421  // check whether coarse solver is set properly. If not, set default coarse solver.
422  if (!coarseList.isParameter("smoother: type"))
423  coarseList.set("smoother: type", "Amesos-KLU"); // set default coarse solver according to ML 5.0 guide
424  RCP<SmootherFactory> coarseFact = GetSmootherFactory(coarseList, Teuchos::null);
425 
426  // Smoothers Top Level Parameters
427 
428  RCP<ParameterList> topLevelSmootherParam = ExtractSetOfParameters(paramList, "smoother");
429 
430  //
431 
432  // Prepare factory managers
433  // TODO: smootherFact can be reuse accross level if same parameters/no specific parameterList
434 
435  for (int levelID = 0; levelID < maxLevels; levelID++) {
436  //
437  // Level FactoryManager
438  //
439 
440  RCP<FactoryManager> manager = rcp(new FactoryManager());
441  if (setKokkosRefactor)
442  manager->SetKokkosRefactor(useKokkosRefactor);
443 
444  //
445  // Smoothers
446  //
447 
448  {
449  // Merge level-specific parameters with global parameters. level-specific parameters takes precedence.
450  // TODO: unit-test this part alone
451 
452  ParameterList levelSmootherParam = GetMLSubList(paramList, "smoother", levelID); // copy
453  MergeParameterList(*topLevelSmootherParam, levelSmootherParam, false); /* false = do no overwrite levelSmootherParam parameters by topLevelSmootherParam parameters */
454  // std::cout << std::endl << "Merged List for level " << levelID << std::endl;
455  // std::cout << levelSmootherParam << std::endl;
456 
457  RCP<SmootherFactory> smootherFact = GetSmootherFactory(levelSmootherParam, Teuchos::null); // TODO: missing AFact input arg.
458 
459  manager->SetFactory("Smoother", smootherFact);
460  }
461 
462  //
463  // Misc
464  //
465 
466  manager->SetFactory("CoarseSolver", coarseFact); // TODO: should not be done in the loop
467  manager->SetFactory("Graph", dropFact);
468  manager->SetFactory("Aggregates", AggFact);
469  manager->SetFactory("DofsPerNode", dropFact);
470  manager->SetFactory("Ptent", PtentFact);
471 
472 #if defined(HAVE_MUELU_ISORROPIA) && defined(HAVE_MPI)
473  if (bDoRepartition == 1) {
474  manager->SetFactory("A", RebalancedAFact);
475  manager->SetFactory("P", RebalancedPFact);
476  manager->SetFactory("R", RebalancedRFact);
477  manager->SetFactory("Nullspace", RebalancedPFact);
478  manager->SetFactory("Importer", RepartitionFact);
479  } else {
480 #endif // #ifdef HAVE_MUELU_ISORROPIA
481  manager->SetFactory("Nullspace", nspFact); // use same nullspace factory throughout all multigrid levels
482  manager->SetFactory("A", AcFact); // same RAP factory for all levels
483  manager->SetFactory("P", PFact); // same prolongator and restrictor factories for all levels
484  manager->SetFactory("R", RFact); // same prolongator and restrictor factories for all levels
485 #if defined(HAVE_MUELU_ISORROPIA) && defined(HAVE_MPI)
486  }
487 #endif
488 
489  this->AddFactoryManager(levelID, 1, manager);
490  } // for (level loop)
491 }
492 
493 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
495  // if nullspace_ has already been extracted from ML parameter list
496  // make nullspace available for MueLu
497  if (nullspace_ != NULL) {
498  RCP<Level> fineLevel = H.GetLevel(0);
499  RCP<Operator> Op = fineLevel->Get<RCP<Operator> >("A");
500  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Op);
501  if (!A.is_null()) {
502  const RCP<const Map> rowMap = fineLevel->Get<RCP<Matrix> >("A")->getRowMap();
503  RCP<MultiVector> nullspace = MultiVectorFactory::Build(rowMap, nullspaceDim_, true);
504 
505  for (size_t i = 0; i < Teuchos::as<size_t>(nullspaceDim_); i++) {
506  Teuchos::ArrayRCP<Scalar> nullspacei = nullspace->getDataNonConst(i);
507  const size_t myLength = nullspace->getLocalLength();
508 
509  for (size_t j = 0; j < myLength; j++) {
510  nullspacei[j] = nullspace_[i * myLength + j];
511  }
512  }
513 
514  fineLevel->Set("Nullspace", nullspace);
515  }
516  }
517 
518  // Do the same for coordinates
519  size_t num_coords = 0;
520  double* coordPTR[3];
521  if (xcoord_) {
522  coordPTR[0] = xcoord_;
523  num_coords++;
524  if (ycoord_) {
525  coordPTR[1] = ycoord_;
526  num_coords++;
527  if (zcoord_) {
528  coordPTR[2] = zcoord_;
529  num_coords++;
530  }
531  }
532  }
533  if (num_coords) {
534  Teuchos::RCP<Level> fineLevel = H.GetLevel(0);
535  Teuchos::RCP<Operator> Op = fineLevel->Get<RCP<Operator> >("A");
536  Teuchos::RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Op);
537  if (!A.is_null()) {
538  const Teuchos::RCP<const Map> rowMap = fineLevel->Get<RCP<Matrix> >("A")->getRowMap();
539  Teuchos::RCP<MultiVector> coordinates = MultiVectorFactory::Build(rowMap, num_coords, true);
540 
541  for (size_t i = 0; i < num_coords; i++) {
542  Teuchos::ArrayRCP<Scalar> coordsi = coordinates->getDataNonConst(i);
543  const size_t myLength = coordinates->getLocalLength();
544  for (size_t j = 0; j < myLength; j++) {
545  coordsi[j] = coordPTR[i][j];
546  }
547  }
548  fineLevel->Set("Coordinates", coordinates);
549  }
550  }
551 
553 }
554 
555 // TODO: code factorization with MueLu_ParameterListInterpreter.
556 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
560  const RCP<FactoryBase>& AFact) {
561  typedef Teuchos::ScalarTraits<Scalar> STS;
562  SC one = STS::one();
563 
564  std::string type = "symmetric Gauss-Seidel"; // default
565 
566  //
567  // Get 'type'
568  //
569 
570  // //TODO: fix defaults!!
571 
572  // // Default coarse grid smoother
573  // std::string type;
574  // if ("smoother" == "coarse") {
575  // #if (defined(HAVE_MUELU_EPETRA) && defined( HAVE_MUELU_AMESOS)) || (defined(HAVE_MUELU_AMESOS2)) // FIXME: test is wrong (ex: compiled with Epetra&&Tpetra&&Amesos2 but without Amesos => error running Epetra problem)
576  // type = ""; // use default defined by AmesosSmoother or Amesos2Smoother
577  // #else
578  // type = "symmetric Gauss-Seidel"; // use a sym Gauss-Seidel (with no damping) as fallback "coarse solver" (TODO: needs Ifpack(2))
579  // #endif
580  // } else {
581  // // TODO: default smoother?
582  // type = "";
583  // }
584 
585  if (paramList.isParameter("smoother: type")) type = paramList.get<std::string>("smoother: type");
586  TEUCHOS_TEST_FOR_EXCEPTION(type.empty(), Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: no \"smoother: type\" in the smoother parameter list" << std::endl
587  << paramList);
588 
589  //
590  // Create the smoother prototype
591  //
592 
593  RCP<SmootherPrototype> smooProto;
594  std::string ifpackType;
595  Teuchos::ParameterList smootherParamList;
596 
597  if (type == "Jacobi" || type == "Gauss-Seidel" || type == "symmetric Gauss-Seidel") {
598  if (type == "symmetric Gauss-Seidel") type = "Symmetric Gauss-Seidel"; // FIXME
599 
600  ifpackType = "RELAXATION";
601  smootherParamList.set("relaxation: type", type);
602 
603  MUELU_COPY_PARAM(paramList, "smoother: sweeps", int, 2, smootherParamList, "relaxation: sweeps");
604  MUELU_COPY_PARAM(paramList, "smoother: damping factor", Scalar, one, smootherParamList, "relaxation: damping factor");
605 
606  smooProto = rcp(new TrilinosSmoother(ifpackType, smootherParamList, 0));
607  smooProto->SetFactory("A", AFact);
608 
609  } else if (type == "Chebyshev" || type == "MLS") {
610  ifpackType = "CHEBYSHEV";
611 
612  MUELU_COPY_PARAM(paramList, "smoother: sweeps", int, 2, smootherParamList, "chebyshev: degree");
613  if (paramList.isParameter("smoother: MLS alpha")) {
614  MUELU_COPY_PARAM(paramList, "smoother: MLS alpha", double, 20, smootherParamList, "chebyshev: ratio eigenvalue");
615  } else {
616  MUELU_COPY_PARAM(paramList, "smoother: Chebyshev alpha", double, 20, smootherParamList, "chebyshev: ratio eigenvalue");
617  }
618 
619  smooProto = rcp(new TrilinosSmoother(ifpackType, smootherParamList, 0));
620  smooProto->SetFactory("A", AFact);
621 
622  } else if (type == "Hiptmair") {
623  ifpackType = "HIPTMAIR";
624  std::string subSmootherType = "Chebyshev";
625  if (paramList.isParameter("subsmoother: type"))
626  subSmootherType = paramList.get<std::string>("subsmoother: type");
627  std::string subSmootherIfpackType;
628  if (subSmootherType == "Chebyshev")
629  subSmootherIfpackType = "CHEBYSHEV";
630  else if (subSmootherType == "Jacobi" || subSmootherType == "Gauss-Seidel" || subSmootherType == "symmetric Gauss-Seidel") {
631  if (subSmootherType == "symmetric Gauss-Seidel") subSmootherType = "Symmetric Gauss-Seidel"; // FIXME
632  subSmootherIfpackType = "RELAXATION";
633  } else
634  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: unknown smoother type. '" << subSmootherType << "' not supported by MueLu.");
635 
636  smootherParamList.set("hiptmair: smoother type 1", subSmootherIfpackType);
637  smootherParamList.set("hiptmair: smoother type 2", subSmootherIfpackType);
638 
639  auto smoother1ParamList = smootherParamList.sublist("hiptmair: smoother list 1");
640  auto smoother2ParamList = smootherParamList.sublist("hiptmair: smoother list 2");
641 
642  if (subSmootherType == "Chebyshev") {
643  MUELU_COPY_PARAM(paramList, "subsmoother: edge sweeps", int, 2, smoother1ParamList, "chebyshev: degree");
644  MUELU_COPY_PARAM(paramList, "subsmoother: node sweeps", int, 2, smoother2ParamList, "chebyshev: degree");
645 
646  MUELU_COPY_PARAM(paramList, "subsmoother: Chebyshev", double, 20, smoother1ParamList, "chebyshev: ratio eigenvalue");
647  MUELU_COPY_PARAM(paramList, "subsmoother: Chebyshev", double, 20, smoother2ParamList, "chebyshev: ratio eigenvalue");
648  } else {
649  MUELU_COPY_PARAM(paramList, "subsmoother: edge sweeps", int, 2, smoother1ParamList, "relaxation: sweeps");
650  MUELU_COPY_PARAM(paramList, "subsmoother: node sweeps", int, 2, smoother2ParamList, "relaxation: sweeps");
651 
652  MUELU_COPY_PARAM(paramList, "subsmoother: SGS damping factor", double, 0.8, smoother2ParamList, "relaxation: damping factor");
653  }
654 
655  smooProto = rcp(new TrilinosSmoother(ifpackType, smootherParamList, 0));
656  smooProto->SetFactory("A", AFact);
657 
658  } else if (type == "IFPACK") { // TODO: this option is not described in the ML Guide v5.0
659 
660 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_IFPACK)
661  ifpackType = paramList.get<std::string>("smoother: ifpack type");
662 
663  if (ifpackType == "ILU") {
664  // TODO fix this (type mismatch double vs. int)
665  // MUELU_COPY_PARAM(paramList, "smoother: ifpack level-of-fill", double /*int*/, 0.0 /*2*/, smootherParamList, "fact: level-of-fill");
666  if (paramList.isParameter("smoother: ifpack level-of-fill"))
667  smootherParamList.set("fact: level-of-fill", Teuchos::as<int>(paramList.get<double>("smoother: ifpack level-of-fill")));
668  else
669  smootherParamList.set("fact: level-of-fill", as<int>(0));
670 
671  MUELU_COPY_PARAM(paramList, "smoother: ifpack overlap", int, 2, smootherParamList, "partitioner: overlap");
672 
673  // TODO change to TrilinosSmoother as soon as Ifpack2 supports all preconditioners from Ifpack
674  smooProto =
675  MueLu::GetIfpackSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>(ifpackType,
676  smootherParamList,
677  paramList.get<int>("smoother: ifpack overlap"));
678  smooProto->SetFactory("A", AFact);
679  } else {
680  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: unknown ML smoother type " + type + " (IFPACK) not supported by MueLu. Only ILU is supported.");
681  }
682 #else
683  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: MueLu compiled without Ifpack support");
684 #endif
685 
686  } else if (type.length() > strlen("Amesos") && type.substr(0, strlen("Amesos")) == "Amesos") { /* catch Amesos-* */
687  std::string solverType = type.substr(strlen("Amesos") + 1); /* ("Amesos-KLU" -> "KLU") */
688 
689  // Validator: following upper/lower case is what is allowed by ML
690  bool valid = false;
691  const int validatorSize = 5;
692  std::string validator[validatorSize] = {"Superlu", "Superludist", "KLU", "UMFPACK", "MUMPS"}; /* TODO: should "" be allowed? */
693  for (int i = 0; i < validatorSize; i++) {
694  if (validator[i] == solverType) valid = true;
695  }
696  TEUCHOS_TEST_FOR_EXCEPTION(!valid, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: unknown smoother type. '" << type << "' not supported.");
697 
698  // FIXME: MueLu should accept any Upper/Lower case. Not the case for the moment
699  std::transform(solverType.begin() + 1, solverType.end(), solverType.begin() + 1, ::tolower);
700 
701  smooProto = Teuchos::rcp(new DirectSolver(solverType, Teuchos::ParameterList()));
702  smooProto->SetFactory("A", AFact);
703 
704  } else {
705  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: unknown smoother type. '" << type << "' not supported by MueLu.");
706  }
707  TEUCHOS_TEST_FOR_EXCEPTION(smooProto == Teuchos::null, Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter: no smoother prototype. fatal error.");
708 
709  //
710  // Create the smoother factory
711  //
712 
713  RCP<SmootherFactory> SmooFact = rcp(new SmootherFactory());
714 
715  // Set parameters of the smoother factory
716  MUELU_READ_PARAM(paramList, "smoother: pre or post", std::string, "both", preOrPost);
717  if (preOrPost == "both") {
718  SmooFact->SetSmootherPrototypes(smooProto, smooProto);
719  } else if (preOrPost == "pre") {
720  SmooFact->SetSmootherPrototypes(smooProto, Teuchos::null);
721  } else if (preOrPost == "post") {
722  SmooFact->SetSmootherPrototypes(Teuchos::null, smooProto);
723  }
724 
725  return SmooFact;
726 }
727 
728 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
730  // check if it's a TwoLevelFactoryBase based transfer factory
731  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast, "Transfer factory is not derived from TwoLevelFactoryBase. Since transfer factories will be handled by the RAPFactory they have to be derived from TwoLevelFactoryBase!");
732  TransferFacts_.push_back(factory);
733 }
734 
735 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
737  return TransferFacts_.size();
738 }
739 
740 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
742  try {
743  Matrix& A = dynamic_cast<Matrix&>(Op);
744  if (A.IsFixedBlockSizeSet() && (A.GetFixedBlockSize() != blksize_))
745  this->GetOStream(Warnings0) << "Setting matrix block size to " << blksize_ << " (value of the parameter in the list) "
746  << "instead of " << A.GetFixedBlockSize() << " (provided matrix)." << std::endl;
747 
748  A.SetFixedBlockSize(blksize_);
749 
750 #ifdef HAVE_MUELU_DEBUG
751  MatrixUtils::checkLocalRowMapMatchesColMap(A);
752 #endif // HAVE_MUELU_DEBUG
753 
754  } catch (std::bad_cast&) {
755  this->GetOStream(Warnings0) << "Skipping setting block size as the operator is not a matrix" << std::endl;
756  }
757 }
758 
759 } // namespace MueLu
760 
761 #define MUELU_MLPARAMETERLISTINTERPRETER_SHORT
762 #endif /* MUELU_MLPARAMETERLISTINTERPRETER_DEF_HPP */
763 
764 // TODO: see if it can be factorized with ML interpreter (ex: generation of Ifpack param list)
Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Exception indicating invalid cast attempted.
This class specifies the default factory that should generate some data on a Level if the data does n...
void MergeParameterList(const Teuchos::ParameterList &source, Teuchos::ParameterList &dest, bool overWrite)
: merge two parameter lists
Factory for determing the number of partitions for rebalancing.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
T & get(const std::string &name, T def_value)
Class that encapsulates external library smoothers.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Interface to IsorropiaInterface to Isorropia allowing to access other rebalancing/repartitioning algo...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Factory for building permutation matrix that can be be used to shuffle data (matrices, vectors) among processes.
void SetKokkosRefactor(const bool useKokkos)
std::string tolower(const std::string &str)
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
#define MUELU_READ_PARAM(paramList, paramStr, varType, defaultValue, varName)
const Teuchos::ParameterList & GetMLSubList(const Teuchos::ParameterList &paramList, const std::string &type, int levelID)
void AddTransferFactory(const RCP< FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories for RAPFactory.
Factory for building tentative prolongator.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
Factory for building restriction operators using a prolongator factory.
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
bool isParameter(const std::string &name) const
void CreateSublists(const ParameterList &List, ParameterList &newList)
bool remove(std::string const &name, bool throwIfNotExists=true)
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
void SetSmootherPrototypes(RCP< SmootherPrototype > preAndPostSmootherPrototype)
Set smoother prototypes.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
static RCP< SmootherFactory > GetSmootherFactory(const Teuchos::ParameterList &paramList, const RCP< FactoryBase > &AFact=Teuchos::null)
Read smoother options and build the corresponding smoother factory.
virtual void SetupOperator(Operator &Op) const
Setup Operator object.
Teuchos::RCP< Teuchos::ParameterList > ExtractSetOfParameters(const Teuchos::ParameterList &paramList, const std::string &str)
size_t NumTransferFactories() const
Returns number of transfer factories.
AmalgamationFactory for subblocks of strided map based amalgamation data.
Applies permutation to grid transfer operators.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void SetParameter(const std::string &name, const ParameterEntry &entry)
Set a parameter directly as a ParameterEntry.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Set Factory.
Factory for creating a graph based on a given matrix.
Helper class which transforms an &quot;AmalgamatedPartition&quot; array to an unamalgamated &quot;Partition&quot;...
void SetParameterList(const Teuchos::ParameterList &paramList)
virtual void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)=0
Configuration.
Scalar SC
bool isType(const std::string &name) const
Factory for building restriction operators.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Factory for creating a graph based on a given matrix.
Factory for building coarse matrices.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
Factory for building coarse matrices.
#define MUELU_COPY_PARAM(paramList, paramStr, varType, defaultValue, outParamList, outParamStr)
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Factory for building Smoothed Aggregation prolongators.Input/output of SaPFactory
Factory for building uncoupled aggregates.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Factory for generating nullspace.
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
bool is_null() const