48 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_IFPACK) 50 #include <Ifpack_Chebyshev.h> 51 #include "Xpetra_MultiVectorFactory.hpp" 56 #include "MueLu_Utilities.hpp" 63 : type_(type), overlap_(overlap)
81 ParameterList& paramList =
const_cast<ParameterList&
>(this->
GetParameterList());
82 paramList.setParameters(list);
86 prec_->SetParameters(*precList);
112 template <
class Node>
114 this->
Input(currentLevel,
"A");
116 if (
type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
117 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
118 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
119 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
120 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
121 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
122 this->
Input(currentLevel,
"CoarseNumZLayers");
123 this->
Input(currentLevel,
"LineDetection_VertLineIds");
127 template <
class Node>
131 this->
GetOStream(
Warnings0) <<
"MueLu::IfpackSmoother::Setup(): Setup() has already been called";
133 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
135 double lambdaMax = -1.0;
136 if (
type_ ==
"Chebyshev") {
137 std::string maxEigString =
"chebyshev: max eigenvalue";
138 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
141 lambdaMax = Teuchos::getValue<Scalar>(this->
GetParameter(maxEigString));
142 this->
GetOStream(
Statistics1) << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
144 }
catch (Teuchos::Exceptions::InvalidParameterName) {
145 lambdaMax =
A_->GetMaxEigenvalueEstimate();
147 if (lambdaMax != -1.0) {
148 this->
GetOStream(
Statistics1) << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
149 this->
SetParameter(maxEigString, ParameterEntry(lambdaMax));
154 const Scalar defaultEigRatio = 20;
156 Scalar ratio = defaultEigRatio;
158 ratio = Teuchos::getValue<Scalar>(this->
GetParameter(eigRatioString));
160 }
catch (Teuchos::Exceptions::InvalidParameterName) {
161 this->
SetParameter(eigRatioString, ParameterEntry(ratio));
164 if (currentLevel.GetLevelID()) {
169 RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >(
"A");
170 size_t nRowsFine = fineA->getGlobalNumRows();
171 size_t nRowsCoarse =
A_->getGlobalNumRows();
173 ratio = std::max(ratio, as<Scalar>(nRowsFine)/nRowsCoarse);
176 this->
SetParameter(eigRatioString, ParameterEntry(ratio));
180 if (
type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
181 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
182 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
183 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
184 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
185 type_ ==
"LINESMOOTHING_BLOCKRELAXATION" ) {
186 ParameterList& myparamList =
const_cast<ParameterList&
>(this->
GetParameterList());
188 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,
"CoarseNumZLayers");
189 if (CoarseNumZLayers > 0) {
190 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel,
"LineDetection_VertLineIds");
194 for(
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
195 if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
198 size_t numLocalRows =
A_->getNodeNumRows();
199 TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0,
Exceptions::RuntimeError,
"MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
201 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
202 myparamList.set(
"partitioner: type",
"user");
203 myparamList.set(
"partitioner: map",&(TVertLineIdSmoo[0]));
204 myparamList.set(
"partitioner: local parts",maxPart+1);
207 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
210 Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
211 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
212 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
213 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
214 myparamList.set(
"partitioner: type",
"user");
215 myparamList.set(
"partitioner: map",&(partitionerMap[0]));
216 myparamList.set(
"partitioner: local parts",maxPart + 1);
219 if (
type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
220 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
221 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
222 type_ =
"block relaxation";
224 type_ =
"block relaxation";
227 this->
GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
228 myparamList.remove(
"partitioner: type",
false);
229 myparamList.remove(
"partitioner: map",
false);
230 myparamList.remove(
"partitioner: local parts",
false);
231 type_ =
"point relaxation stand-alone";
246 if (
type_ ==
"Chebyshev" && lambdaMax == -1.0) {
248 if (chebyPrec != Teuchos::null) {
250 A_->SetMaxEigenvalueEstimate(lambdaMax);
251 this->
GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack)" <<
" = " << lambdaMax << std::endl;
253 TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == -1.0,
Exceptions::RuntimeError,
"MueLu::IfpackSmoother::Setup(): no maximum eigenvalue estimate");
259 template <
class Node>
264 Teuchos::ParameterList paramList;
265 bool supportInitialGuess =
false;
266 if (
type_ ==
"Chebyshev") {
267 paramList.set(
"chebyshev: zero starting solution", InitialGuessIsZero);
268 supportInitialGuess =
true;
270 }
else if (
type_ ==
"point relaxation stand-alone") {
271 paramList.set(
"relaxation: zero starting solution", InitialGuessIsZero);
272 supportInitialGuess =
true;
278 if (InitialGuessIsZero || supportInitialGuess) {
282 prec_->ApplyInverse(epB, epX);
286 RCP<MultiVector> Correction = MultiVectorFactory::Build(
A_->getDomainMap(), X.getNumVectors());
291 prec_->ApplyInverse(epB, epX);
293 X.update(1.0, *Correction, 1.0);
297 template <
class Node>
304 template <
class Node>
306 std::ostringstream out;
311 out <<
"{type = " <<
type_ <<
"}";
313 out <<
prec_->Label();
318 template <
class Node>
323 out0 <<
"Prec. type: " <<
type_ << std::endl;
326 out0 <<
"Parameter list: " << std::endl;
327 Teuchos::OSTab tab2(out);
329 out0 <<
"Overlap: " <<
overlap_ << std::endl;
333 if (
prec_ != Teuchos::null) {
334 Teuchos::OSTab tab2(out);
335 out << *
prec_ << std::endl;
338 if (verbLevel &
Debug) {
341 <<
"RCP<A_>: " <<
A_ << std::endl
342 <<
"RCP<prec_>: " <<
prec_ << std::endl;
351 #if defined(HAVE_MUELU_EPETRA) void Setup(Level ¤tLevel)
Set up the smoother.
Important warning messages (one line)
RCP< SmootherPrototype > Copy() const
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Print external lib objects.
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print additional debugging information.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
virtual const Teuchos::ParameterList & GetParameterList() const
Ifpack_Preconditioner * Create(const std::string PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
Print skeleton for the run, i.e. factory calls and used parameters.
IfpackSmoother(std::string const &type, Teuchos::ParameterList const ¶mList=Teuchos::ParameterList(), LO const &overlap=0)
Constructor.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
VerbLevel GetVerbLevel() const
Get the verbosity level.
Print statistics that do not involve significant additional computation.
RCP< ParameterList > RemoveFactoriesFromList(const ParameterList &list) const
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
void DeclareInput(Level ¤tLevel) const
Input.
RCP< Ifpack_Preconditioner > prec_
pointer to Ifpack solver object
std::string type_
ifpack-specific key phrase that denote smoother type
virtual void SetParameterList(const ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
Class that holds all level-specific information.
LO overlap_
overlap when using the smoother in additive Schwarz mode
bool IsSetup() const
Get the state of a smoother prototype.
RCP< Matrix > A_
Matrix. Not used directly, but held inside of prec_. So we have to keep an RCP pointer to it! ...
virtual double GetLambdaMax()
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
void SetParameter(const std::string &name, const ParameterEntry &entry)
Set a parameter directly as a ParameterEntry.
void SetParameterList(const Teuchos::ParameterList ¶mList)
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
std::string description() const
Return a simple one-line description of this object.
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
const ParameterEntry & GetParameter(const std::string &name) const
Retrieves a const entry with the name name.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
void Input(Level &level, const std::string &varName) const
virtual std::string description() const
Return a simple one-line description of this object.
Class that encapsulates Ifpack smoothers.