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