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