8#ifndef MUELU_ADAPTIVESAMLPARAMETERLISTINTERPRETER_DEF_HPP_
9#define MUELU_ADAPTIVESAMLPARAMETERLISTINTERPRETER_DEF_HPP_
11#include <Teuchos_XMLParameterListHelpers.hpp>
14#if defined(HAVE_MUELU_ML) && defined(HAVE_MUELU_EPETRA)
15#include <ml_ValidateParameters.h>
18#include <Xpetra_Matrix.hpp>
19#include <Xpetra_MultiVector.hpp>
20#include <Xpetra_MultiVectorFactory.hpp>
21#include <Xpetra_Operator.hpp>
22#include <Xpetra_IO.hpp>
27#include "MueLu_Hierarchy.hpp"
28#include "MueLu_FactoryManager.hpp"
30#include "MueLu_TentativePFactory.hpp"
31#include "MueLu_SaPFactory.hpp"
32#include "MueLu_PgPFactory.hpp"
33#include "MueLu_TransPFactory.hpp"
34#include "MueLu_GenericRFactory.hpp"
35#include "MueLu_SmootherPrototype.hpp"
36#include "MueLu_SmootherFactory.hpp"
37#include "MueLu_TrilinosSmoother.hpp"
39#include "MueLu_DirectSolver.hpp"
40#include "MueLu_HierarchyUtils.hpp"
41#include "MueLu_RAPFactory.hpp"
42#include "MueLu_CoalesceDropFactory.hpp"
43#include "MueLu_CoupledAggregationFactory.hpp"
44#include "MueLu_UncoupledAggregationFactory.hpp"
45#include "MueLu_HybridAggregationFactory.hpp"
46#include "MueLu_NullspaceFactory.hpp"
48#include "MueLu_MLParameterListInterpreter.hpp"
60#define MUELU_READ_PARAM(paramList, paramStr, varType, defaultValue, varName) \
61 varType varName = defaultValue; if (paramList.isParameter(paramStr)) varName = paramList.get<varType>(paramStr);
64#define MUELU_COPY_PARAM(paramList, paramStr, varType, defaultValue, outParamList, outParamStr) \
65 if (paramList.isParameter(paramStr)) \
66 outParamList.set<varType>(outParamStr, paramList.get<varType>(paramStr)); \
67 else outParamList.set<varType>(outParamStr, defaultValue); \
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 Teuchos::RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile(xmlFileName);
82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 Teuchos::ParameterList paramList = paramList_in;
86 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
99 MUELU_READ_PARAM(paramList,
"aggregation: type", std::string,
"Uncoupled", agg_type);
101 MUELU_READ_PARAM(paramList,
"aggregation: damping factor",
double, (
double)4/(
double)3, agg_damping);
103 MUELU_READ_PARAM(paramList,
"aggregation: nodes per aggregate",
int, 1, minPerAgg);
105 MUELU_READ_PARAM(paramList,
"null space: type", std::string,
"default vectors", nullspaceType);
106 MUELU_READ_PARAM(paramList,
"null space: dimension",
int, -1, nullspaceDim);
107 MUELU_READ_PARAM(paramList,
"null space: vectors",
double*, NULL, nullspaceVec);
109 MUELU_READ_PARAM(paramList,
"energy minimization: enable",
bool,
false, bEnergyMinimization);
118 ParameterList paramListWithSubList;
120 paramList = paramListWithSubList;
125 int maxNbrAlreadySelected = 0;
128 this->blksize_ = nDofsPerNode;
131 Teuchos::EVerbosityLevel eVerbLevel = Teuchos::VERB_NONE;
132 if (verbosityLevel == 0) eVerbLevel = Teuchos::VERB_NONE;
133 if (verbosityLevel > 0) eVerbLevel = Teuchos::VERB_LOW;
134 if (verbosityLevel > 4) eVerbLevel = Teuchos::VERB_MEDIUM;
135 if (verbosityLevel > 7) eVerbLevel = Teuchos::VERB_HIGH;
136 if (verbosityLevel > 9) eVerbLevel = Teuchos::VERB_EXTREME;
138 TEUCHOS_TEST_FOR_EXCEPTION(agg_type !=
"Uncoupled" && agg_type !=
"Coupled",
Exceptions::RuntimeError,
"MueLu::MLParameterListInterpreter::Setup(): parameter \"aggregation: type\": only 'Uncoupled' or 'Coupled' aggregation is supported.");
145 RCP<FactoryBase> CoupledAggFact = Teuchos::null;
146 if(agg_type ==
"Uncoupled") {
149 CoupledAggFact2->SetMinNodesPerAggregate(minPerAgg);
150 CoupledAggFact2->SetMaxNeighAlreadySelected(maxNbrAlreadySelected);
151 CoupledAggFact2->SetOrdering(
"natural");
152 CoupledAggFact = CoupledAggFact2;
156 CoupledAggFact2->SetMinNodesPerAggregate(minPerAgg);
157 CoupledAggFact2->SetMaxNeighAlreadySelected(maxNbrAlreadySelected);
158 CoupledAggFact2->SetOrdering(
"natural");
159 CoupledAggFact2->SetPhase3AggCreation(0.5);
160 CoupledAggFact = CoupledAggFact2;
162 if (verbosityLevel > 3) {
163 *out <<
"========================= Aggregate option summary =========================" << std::endl;
164 *out <<
"min Nodes per aggregate : " << minPerAgg << std::endl;
165 *out <<
"min # of root nbrs already aggregated : " << maxNbrAlreadySelected << std::endl;
166 *out <<
"aggregate ordering : natural" << std::endl;
167 *out <<
"=============================================================================" << std::endl;
173 if (agg_damping == 0.0 && bEnergyMinimization ==
false) {
177 }
else if (agg_damping != 0.0 && bEnergyMinimization ==
false) {
179 RCP<SaPFactory> SaPFact = rcp(
new SaPFactory() );
180 SaPFact->SetParameter(
"sa: damping factor", ParameterEntry(agg_damping));
183 }
else if (bEnergyMinimization ==
true) {
189 RCP<RAPFactory> AcFact = rcp(
new RAPFactory() );
190 for (
size_t i = 0; i<TransferFacts_.size(); i++) {
191 AcFact->AddTransferFactory(TransferFacts_[i]);
201 if (nullspaceType !=
"default vectors") {
202 TEUCHOS_TEST_FOR_EXCEPTION(nullspaceType !=
"pre-computed",
Exceptions::RuntimeError,
"MueLu::MLParameterListInterpreter: no valid nullspace (no pre-computed null space). error.");
203 TEUCHOS_TEST_FOR_EXCEPTION(nullspaceDim == -1,
Exceptions::RuntimeError,
"MueLu::MLParameterListInterpreter: no valid nullspace (nullspace dim == -1). error.");
204 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\'");
206 nullspaceDim_ = nullspaceDim;
207 nullspace_ = nullspaceVec;
210 Teuchos::RCP<NullspaceFactory> nspFact = Teuchos::rcp(
new NullspaceFactory());
211 nspFact->SetFactory(
"Nullspace", PtentFact);
219 this->numDesiredLevel_ = maxLevels;
220 this->maxCoarseSize_ = maxCoarseSize;
223 RCP<SmootherFactory> initSmootherFact = Teuchos::null;
224 if(paramList.isSublist(
"init smoother")) {
225 ParameterList& initList = paramList.sublist(
"init smoother");
228 std::string ifpackType =
"RELAXATION";
229 Teuchos::ParameterList smootherParamList;
230 smootherParamList.set(
"relaxation: type",
"symmetric Gauss-Seidel");
231 smootherParamList.set(
"smoother: sweeps", 1);
232 smootherParamList.set(
"smoother: damping factor", 1.0);
233 RCP<SmootherPrototype> smooProto = rcp(
new TrilinosSmoother(ifpackType, smootherParamList, 0) );
236 initSmootherFact->SetSmootherPrototypes(smooProto, smooProto);
242 ParameterList& coarseList = paramList.sublist(
"coarse: list");
258 for (
int levelID=0; levelID < maxLevels; levelID++) {
275 ParameterList levelSmootherParam =
GetMLSubList(paramList,
"smoother", levelID);
282 manager->SetFactory(
"Smoother", smootherFact);
283 smootherFact->DisableMultipleCallCheck();
285 initmanager->SetFactory(
"Smoother", initSmootherFact);
286 initmanager->SetFactory(
"CoarseSolver", initSmootherFact);
287 initSmootherFact->DisableMultipleCallCheck();
295 Teuchos::rcp_dynamic_cast<PFactory>(PFact)->DisableMultipleCallCheck();
296 Teuchos::rcp_dynamic_cast<PFactory>(PtentFact)->DisableMultipleCallCheck();
297 Teuchos::rcp_dynamic_cast<TwoLevelFactoryBase>(RFact)->DisableMultipleCallCheck();
298 Teuchos::rcp_dynamic_cast<SingleLevelFactoryBase>(coarseFact)->DisableMultipleCallCheck();
299 Teuchos::rcp_dynamic_cast<SingleLevelFactoryBase>(dropFact)->DisableMultipleCallCheck();
300 Teuchos::rcp_dynamic_cast<SingleLevelFactoryBase>(CoupledAggFact)->DisableMultipleCallCheck();
301 Teuchos::rcp_dynamic_cast<TwoLevelFactoryBase>(AcFact)->DisableMultipleCallCheck();
302 Teuchos::rcp_dynamic_cast<SingleLevelFactoryBase>(nspFact)->DisableMultipleCallCheck();
304 manager->SetFactory(
"CoarseSolver", coarseFact);
305 manager->SetFactory(
"Graph", dropFact);
306 manager->SetFactory(
"Aggregates", CoupledAggFact);
307 manager->SetFactory(
"DofsPerNode", dropFact);
308 manager->SetFactory(
"A", AcFact);
309 manager->SetFactory(
"P", PFact);
310 manager->SetFactory(
"Ptent", PtentFact);
311 manager->SetFactory(
"R", RFact);
312 manager->SetFactory(
"Nullspace", nspFact);
315 initmanager->SetFactory(
"Graph", dropFact);
316 initmanager->SetFactory(
"Aggregates", CoupledAggFact);
317 initmanager->SetFactory(
"DofsPerNode", dropFact);
318 initmanager->SetFactory(
"A", AcFact);
319 initmanager->SetFactory(
"P", PtentFact);
320 initmanager->SetFactory(
"Ptent", PtentFact);
321 initmanager->SetFactory(
"R", RFact);
322 initmanager->SetFactory(
"Nullspace", nspFact);
324 this->AddFactoryManager(levelID, 1, manager);
325 this->AddInitFactoryManager(levelID, 1, initmanager);
329 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
334 RCP<Operator> Op = l->Get<RCP<Operator> >(
"A");
342 int lastLevelID = this->numDesiredLevel_ - 1;
343 bool isLastLevel =
false;
345 while(!isLastLevel) {
346 bool r = H.
Setup(levelID,
347 InitLvlMngr(levelID-1, lastLevelID),
348 InitLvlMngr(levelID, lastLevelID),
349 InitLvlMngr(levelID+1, lastLevelID));
351 isLastLevel = r || (levelID == lastLevelID);
356 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
362 if (this->nullspace_ != NULL) {
363 RCP<Level> fineLevel = H.
GetLevel(0);
364 const RCP<const Map> rowMap = fineLevel->Get< RCP<Matrix> >(
"A")->getRowMap();
365 RCP<MultiVector> nullspace = MultiVectorFactory::Build(rowMap, nullspaceDim_,
true);
367 for (
size_t i=0; i < Teuchos::as<size_t>(nullspaceDim_); i++) {
368 Teuchos::ArrayRCP<Scalar> nullspacei = nullspace->getDataNonConst(i);
369 const size_t myLength = nullspace->getLocalLength();
371 for (
size_t j = 0; j < myLength; j++) {
372 nullspacei[j] = nullspace_[i*myLength + j];
376 fineLevel->Set(
"Nullspace", nullspace);
385 SetupInitHierarchy(H);
389 Teuchos::RCP<MueLu::Level> Finest = H.
GetLevel(0);
390 Teuchos::RCP<MultiVector> nspVector2 = Finest->Get<Teuchos::RCP<MultiVector> >(
"Nullspace");
392 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"orig_nsp.vec", *nspVector2);
394 RCP<Matrix> Op = Finest->Get<RCP<Matrix> >(
"A");
395 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"A.mat", *Op);
398 Teuchos::RCP<MultiVector> homogRhsVec = MultiVectorFactory::Build(nspVector2->getMap(),nspVector2->getNumVectors(),
true);
399 homogRhsVec->putScalar(0.0);
404 H.
Iterate(*homogRhsVec, *nspVector2, 1,
false);
407 Finest->Set(
"Nullspace",nspVector2);
409 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"new_nsp.vec", *nspVector2);
434 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
437 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!");
438 TransferFacts_.push_back(factory);
441 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
443 return TransferFacts_.size();
446 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
449 Matrix& A =
dynamic_cast<Matrix&
>(Op);
450 if (A.IsFixedBlockSizeSet() && (A.GetFixedBlockSize() != blksize_))
451 this->GetOStream(
Warnings0) <<
"Setting matrix block size to " << blksize_ <<
" (value of the parameter in the list) "
452 <<
"instead of " << A.GetFixedBlockSize() <<
" (provided matrix)." << std::endl;
454 A.SetFixedBlockSize(blksize_);
456 }
catch (std::bad_cast& e) {
457 this->GetOStream(
Warnings0) <<
"Skipping setting block size as the operator is not a matrix" << std::endl;
#define MUELU_READ_PARAM(paramList, paramStr, varType, defaultValue, varName)
size_t NumTransferFactories() const
Returns number of transfer factories.
void SetParameterList(const Teuchos::ParameterList ¶mList)
virtual void SetupOperator(Operator &Op) const
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
void SetupInitHierarchy(Hierarchy &H) const
void AddTransferFactory(const RCP< FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories for RAPFactory.
AdaptiveSaMLParameterListInterpreter()
Constructor.
Factory for creating a graph based on a given matrix.
Factory for coarsening a graph with uncoupled aggregation.
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.
RCP< FactoryManagerBase > GetFactoryManager(int levelID) const
size_t getNumFactoryManagers() const
returns number of factory managers stored in levelManagers_ vector.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
bool Setup(int coarseLevelID, const RCP< const FactoryManagerBase > fineLevelManager, const RCP< const FactoryManagerBase > coarseLevelManager, const RCP< const FactoryManagerBase > nextLevelManager=Teuchos::null)
Multi-level setup phase: build a new level of the hierarchy.
ConvergenceStatus Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
void SetMaxCoarseSize(Xpetra::global_size_t maxCoarseSize)
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
static RCP< SmootherFactory > GetSmootherFactory(const Teuchos::ParameterList ¶mList, const RCP< FactoryBase > &AFact=Teuchos::null)
Read smoother options and build the corresponding smoother factory.
Factory for generating nullspace.
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Factory for building coarse matrices.
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 ¶mList, const std::string &str)
@ Warnings0
Important warning messages (one line)
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 ¶mList, const std::string &type, int levelID)
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
Translate Teuchos verbosity level to MueLu verbosity level.