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