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