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_HierarchyHelpers.hpp" 41 #include "MueLu_RAPFactory.hpp" 42 #include "MueLu_CoalesceDropFactory.hpp" 43 #include "MueLu_CoupledAggregationFactory.hpp" 44 #include "MueLu_UncoupledAggregationFactory.hpp" 45 #include "MueLu_NullspaceFactory.hpp" 47 #include "MueLu_MLParameterListInterpreter.hpp" 59 #define MUELU_READ_PARAM(paramList, paramStr, varType, defaultValue, varName) \ 60 varType varName = defaultValue; if (paramList.isParameter(paramStr)) varName = paramList.get<varType>(paramStr); 63 #define MUELU_COPY_PARAM(paramList, paramStr, varType, defaultValue, outParamList, outParamStr) \ 64 if (paramList.isParameter(paramStr)) \ 65 outParamList.set<varType>(outParamStr, paramList.get<varType>(paramStr)); \ 66 else outParamList.set<varType>(outParamStr, defaultValue); \ 70 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 Teuchos::RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile(xmlFileName);
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 Teuchos::ParameterList paramList = paramList_in;
85 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
98 MUELU_READ_PARAM(paramList,
"aggregation: type", std::string,
"Uncoupled", agg_type);
100 MUELU_READ_PARAM(paramList,
"aggregation: damping factor",
double, (
double)4/(
double)3, agg_damping);
102 MUELU_READ_PARAM(paramList,
"aggregation: nodes per aggregate",
int, 1, minPerAgg);
104 MUELU_READ_PARAM(paramList,
"null space: type", std::string,
"default vectors", nullspaceType);
105 MUELU_READ_PARAM(paramList,
"null space: dimension",
int, -1, nullspaceDim);
106 MUELU_READ_PARAM(paramList,
"null space: vectors",
double*, NULL, nullspaceVec);
108 MUELU_READ_PARAM(paramList,
"energy minimization: enable",
bool,
false, bEnergyMinimization);
117 ParameterList paramListWithSubList;
119 paramList = paramListWithSubList;
124 int maxNbrAlreadySelected = 0;
130 Teuchos::EVerbosityLevel eVerbLevel = Teuchos::VERB_NONE;
131 if (verbosityLevel == 0) eVerbLevel = Teuchos::VERB_NONE;
132 if (verbosityLevel > 0) eVerbLevel = Teuchos::VERB_LOW;
133 if (verbosityLevel > 4) eVerbLevel = Teuchos::VERB_MEDIUM;
134 if (verbosityLevel > 7) eVerbLevel = Teuchos::VERB_HIGH;
135 if (verbosityLevel > 9) eVerbLevel = Teuchos::VERB_EXTREME;
137 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.");
144 RCP<FactoryBase> CoupledAggFact = Teuchos::null;
145 if(agg_type ==
"Uncoupled") {
148 CoupledAggFact2->SetMinNodesPerAggregate(minPerAgg);
149 CoupledAggFact2->SetMaxNeighAlreadySelected(maxNbrAlreadySelected);
150 CoupledAggFact2->SetOrdering(
"natural");
151 CoupledAggFact = CoupledAggFact2;
155 CoupledAggFact2->SetMinNodesPerAggregate(minPerAgg);
156 CoupledAggFact2->SetMaxNeighAlreadySelected(maxNbrAlreadySelected);
157 CoupledAggFact2->SetOrdering(
"natural");
158 CoupledAggFact2->SetPhase3AggCreation(0.5);
159 CoupledAggFact = CoupledAggFact2;
161 if (verbosityLevel > 3) {
162 *out <<
"========================= Aggregate option summary =========================" << std::endl;
163 *out <<
"min Nodes per aggregate : " << minPerAgg << std::endl;
164 *out <<
"min # of root nbrs already aggregated : " << maxNbrAlreadySelected << std::endl;
165 *out <<
"aggregate ordering : natural" << std::endl;
166 *out <<
"=============================================================================" << std::endl;
172 if (agg_damping == 0.0 && bEnergyMinimization ==
false) {
176 }
else if (agg_damping != 0.0 && bEnergyMinimization ==
false) {
178 RCP<SaPFactory> SaPFact = rcp(
new SaPFactory() );
179 SaPFact->SetParameter(
"sa: damping factor", ParameterEntry(agg_damping));
182 }
else if (bEnergyMinimization ==
true) {
188 RCP<RAPFactory> AcFact = rcp(
new RAPFactory() );
200 if (nullspaceType !=
"default vectors") {
201 TEUCHOS_TEST_FOR_EXCEPTION(nullspaceType !=
"pre-computed",
Exceptions::RuntimeError,
"MueLu::MLParameterListInterpreter: no valid nullspace (no pre-computed null space). error.");
202 TEUCHOS_TEST_FOR_EXCEPTION(nullspaceDim == -1,
Exceptions::RuntimeError,
"MueLu::MLParameterListInterpreter: no valid nullspace (nullspace dim == -1). error.");
203 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\'");
209 Teuchos::RCP<NullspaceFactory> nspFact = Teuchos::rcp(
new NullspaceFactory());
210 nspFact->SetFactory(
"Nullspace", PtentFact);
222 RCP<SmootherFactory> initSmootherFact = Teuchos::null;
223 if(paramList.isSublist(
"init smoother")) {
224 ParameterList& initList = paramList.sublist(
"init smoother");
227 std::string ifpackType =
"RELAXATION";
228 Teuchos::ParameterList smootherParamList;
229 smootherParamList.set(
"relaxation: type",
"symmetric Gauss-Seidel");
230 smootherParamList.set(
"smoother: sweeps", 1);
231 smootherParamList.set(
"smoother: damping factor", 1.0);
232 RCP<SmootherPrototype> smooProto = rcp(
new TrilinosSmoother(ifpackType, smootherParamList, 0) );
235 initSmootherFact->SetSmootherPrototypes(smooProto, smooProto);
241 ParameterList& coarseList = paramList.sublist(
"coarse: list");
257 for (
int levelID=0; levelID < maxLevels; levelID++) {
274 ParameterList levelSmootherParam =
GetMLSubList(paramList,
"smoother", levelID);
281 manager->SetFactory(
"Smoother", smootherFact);
282 smootherFact->DisableMultipleCallCheck();
284 initmanager->SetFactory(
"Smoother", initSmootherFact);
285 initmanager->SetFactory(
"CoarseSolver", initSmootherFact);
286 initSmootherFact->DisableMultipleCallCheck();
303 manager->SetFactory(
"CoarseSolver", coarseFact);
304 manager->SetFactory(
"Graph", dropFact);
305 manager->SetFactory(
"Aggregates", CoupledAggFact);
306 manager->SetFactory(
"DofsPerNode", dropFact);
307 manager->SetFactory(
"A", AcFact);
308 manager->SetFactory(
"P", PFact);
309 manager->SetFactory(
"Ptent", PtentFact);
310 manager->SetFactory(
"R", RFact);
311 manager->SetFactory(
"Nullspace", nspFact);
314 initmanager->SetFactory(
"Graph", dropFact);
315 initmanager->SetFactory(
"Aggregates", CoupledAggFact);
316 initmanager->SetFactory(
"DofsPerNode", dropFact);
317 initmanager->SetFactory(
"A", AcFact);
318 initmanager->SetFactory(
"P", PtentFact);
319 initmanager->SetFactory(
"Ptent", PtentFact);
320 initmanager->SetFactory(
"R", RFact);
321 initmanager->SetFactory(
"Nullspace", nspFact);
328 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
333 RCP<Operator> Op = l->Get<RCP<Operator> >(
"A");
342 bool isLastLevel =
false;
344 while(!isLastLevel) {
345 bool r = H.
Setup(levelID,
350 isLastLevel = r || (levelID == lastLevelID);
355 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
362 RCP<Level> fineLevel = H.
GetLevel(0);
363 const RCP<const Map> rowMap = fineLevel->Get< RCP<Matrix> >(
"A")->getRowMap();
364 RCP<MultiVector> nullspace = MultiVectorFactory::Build(rowMap,
nullspaceDim_,
true);
366 for (
size_t i=0; i < Teuchos::as<size_t>(
nullspaceDim_); i++) {
367 Teuchos::ArrayRCP<Scalar> nullspacei = nullspace->getDataNonConst(i);
368 const size_t myLength = nullspace->getLocalLength();
370 for (
size_t j = 0; j < myLength; j++) {
375 fineLevel->Set(
"Nullspace", nullspace);
388 Teuchos::RCP<MueLu::Level> Finest = H.
GetLevel(0);
389 Teuchos::RCP<MultiVector> nspVector2 = Finest->Get<Teuchos::RCP<MultiVector> >(
"Nullspace");
391 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"orig_nsp.vec", *nspVector2);
393 RCP<Matrix> Op = Finest->Get<RCP<Matrix> >(
"A");
394 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"A.mat", *Op);
397 Teuchos::RCP<MultiVector> homogRhsVec = MultiVectorFactory::Build(nspVector2->getMap(),nspVector2->getNumVectors(),
true);
398 homogRhsVec->putScalar(0.0);
403 H.
Iterate(*homogRhsVec, *nspVector2, 1,
false);
406 Finest->Set(
"Nullspace",nspVector2);
408 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"new_nsp.vec", *nspVector2);
433 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
436 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!");
440 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
445 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
448 Matrix& A =
dynamic_cast<Matrix&
>(Op);
449 if (A.GetFixedBlockSize() !=
blksize_)
451 <<
"instead of " << A.GetFixedBlockSize() <<
" (provided matrix)." << std::endl;
455 }
catch (std::bad_cast& e) {
456 this->
GetOStream(
Warnings0) <<
"Skipping setting block size as the operator is not a matrix" << std::endl;
virtual void SetupExtra(Hierarchy &H) const
Setup extra data.
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
int blksize_
Matrix configuration storage.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
void SetVerbLevel(const VerbLevel verbLevel)
Set the verbosity level of this object.
void SetParameterList(const Teuchos::ParameterList ¶mList)
Class that encapsulates external library smoothers.
void SetMaxCoarseSize(Xpetra::global_size_t maxCoarseSize)
int nullspaceDim_
nullspace can be embedded in the ML parameter list
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
AdaptiveSaMLParameterListInterpreter()
Constructor.
void AddFactoryManager(int startLevel, int numDesiredLevel, RCP< FactoryManagerBase > manager)
Base class for factories that use two levels (fineLevel and coarseLevel).
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
Namespace for MueLu classes and methods.
const Teuchos::ParameterList & GetMLSubList(const Teuchos::ParameterList ¶mList, const std::string &type, int levelID)
Teuchos::RCP< FactoryManagerBase > InitLvlMngr(int levelID, int lastLevelID) const
Factory for building tentative prolongator.
Factory for coarsening a graph with uncoupled aggregation.
virtual void SetupOperator(Operator &Op) const
Factory for building restriction operators using a prolongator factory.
void CreateSublists(const ParameterList &List, ParameterList &newList)
size_t getNumFactoryManagers() const
returns number of factory managers stored in levelManagers_ vector.
void AddTransferFactory(const RCP< FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories for RAPFactory.
size_t NumTransferFactories() const
Returns number of transfer factories.
static RCP< SmootherFactory > GetSmootherFactory(const Teuchos::ParameterList ¶mList, const RCP< FactoryBase > &AFact=Teuchos::null)
Read smoother options and build the corresponding smoother factory.
Teuchos::RCP< Teuchos::ParameterList > ExtractSetOfParameters(const Teuchos::ParameterList ¶mList, const std::string &str)
RCP< FactoryManagerBase > GetFactoryManager(int levelID) const
void SetupInitHierarchy(Hierarchy &H) const
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Factory for creating a graph base on a given matrix.
std::vector< RCP< FactoryBase > > TransferFacts_
An exception safe way to call the method TwoLevelFactoryBase::DisableMultipleCallCheck.
Factory for building restriction operators.
Exception throws to report errors in the internal logical of the program.
Factory that provides an interface for a concrete implementation of a prolongation operator...
ReturnType Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
Factory for building coarse matrices.
Factory for building Petrov-Galerkin Smoothed Aggregation prolongators.
Factory for building Smoothed Aggregation prolongators.Input/output of SaPFactory
Xpetra::global_size_t maxCoarseSize_
#define MUELU_READ_PARAM(paramList, paramStr, varType, defaultValue, varName)
Factory for building uncoupled aggregates.
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.
void AddInitFactoryManager(int startLevel, int numDesiredLevel, RCP< FactoryManagerBase > manager)
internal routine to add a new factory manager used for the initialization phase
Base class for factories that use one level (currentLevel).
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Factory for generating nullspace.
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.