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