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"
58#include "MueLu_Aggregates.hpp"
65 : type_(type), overlap_(overlap)
78 prec_->SetParameters(
const_cast<ParameterList&
>(this->GetParameterList()));
84 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
85 paramList.setParameters(list);
87 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
89 prec_->SetParameters(*precList);
115 template <
class Node>
117 this->Input(currentLevel,
"A");
119 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
120 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
121 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
122 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
123 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
124 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
125 this->Input(currentLevel,
"CoarseNumZLayers");
126 this->Input(currentLevel,
"LineDetection_VertLineIds");
128 else if (type_ ==
"AGGREGATE")
131 this->Input(currentLevel,
"Aggregates");
136 template <
class Node>
140 this->GetOStream(
Warnings0) <<
"MueLu::IfpackSmoother::Setup(): Setup() has already been called" << std::endl;
142 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
144 double lambdaMax = -1.0;
145 if (type_ ==
"Chebyshev") {
146 std::string maxEigString =
"chebyshev: max eigenvalue";
147 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
150 lambdaMax = Teuchos::getValue<Scalar>(this->GetParameter(maxEigString));
151 this->GetOStream(
Statistics1) << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
153 }
catch (Teuchos::Exceptions::InvalidParameterName&) {
154 lambdaMax = A_->GetMaxEigenvalueEstimate();
156 if (lambdaMax != -1.0) {
157 this->GetOStream(
Statistics1) << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
158 this->SetParameter(maxEigString, ParameterEntry(lambdaMax));
163 const Scalar defaultEigRatio = 20;
165 Scalar ratio = defaultEigRatio;
167 ratio = Teuchos::getValue<Scalar>(this->GetParameter(eigRatioString));
169 }
catch (Teuchos::Exceptions::InvalidParameterName&) {
170 this->SetParameter(eigRatioString, ParameterEntry(ratio));
178 RCP<const Matrix> fineA = currentLevel.
GetPreviousLevel()->Get<RCP<Matrix> >(
"A");
179 size_t nRowsFine = fineA->getGlobalNumRows();
180 size_t nRowsCoarse = A_->getGlobalNumRows();
182 ratio = std::max(ratio, as<Scalar>(nRowsFine)/nRowsCoarse);
184 this->GetOStream(
Statistics1) << eigRatioString <<
" (computed) = " << ratio << std::endl;
185 this->SetParameter(eigRatioString, ParameterEntry(ratio));
189 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
190 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
191 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
192 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
193 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
194 type_ ==
"LINESMOOTHING_BLOCKRELAXATION" ) {
195 ParameterList& myparamList =
const_cast<ParameterList&
>(this->GetParameterList());
197 LO CoarseNumZLayers = currentLevel.
Get<LO>(
"CoarseNumZLayers",
Factory::GetFactory(
"CoarseNumZLayers").get());
198 if (CoarseNumZLayers > 0) {
199 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = currentLevel.
Get< Teuchos::ArrayRCP<LO> >(
"LineDetection_VertLineIds",
Factory::GetFactory(
"LineDetection_VertLineIds").get());
203 for(
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
204 if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
207 size_t numLocalRows = A_->getLocalNumRows();
208 TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0,
Exceptions::RuntimeError,
"MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
210 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
211 myparamList.set(
"partitioner: type",
"user");
212 myparamList.set(
"partitioner: map",&(TVertLineIdSmoo[0]));
213 myparamList.set(
"partitioner: local parts",maxPart+1);
216 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
219 Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
220 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
221 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
222 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
223 myparamList.set(
"partitioner: type",
"user");
224 myparamList.set(
"partitioner: map",&(partitionerMap[0]));
225 myparamList.set(
"partitioner: local parts",maxPart + 1);
228 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
229 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
230 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
231 type_ =
"block relaxation";
233 type_ =
"block relaxation";
236 this->GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
237 myparamList.remove(
"partitioner: type",
false);
238 myparamList.remove(
"partitioner: map",
false);
239 myparamList.remove(
"partitioner: local parts",
false);
240 type_ =
"point relaxation stand-alone";
245 if(type_ ==
"AGGREGATE") {
246 SetupAggregate(currentLevel);
251 ParameterList precList = this->GetParameterList();
252 if(precList.isParameter(
"partitioner: type") && precList.get<std::string>(
"partitioner: type") ==
"linear" &&
253 !precList.isParameter(
"partitioner: local parts")) {
254 precList.set(
"partitioner: local parts", (
int)A_->getLocalNumRows() / A_->GetFixedBlockSize());
261 prec_ = rcp(factory.
Create(type_, &(*epA), overlap_));
262 TEUCHOS_TEST_FOR_EXCEPTION(prec_.is_null(),
Exceptions::RuntimeError,
"Could not create an Ifpack preconditioner with type = \"" << type_ <<
"\"");
269 if (type_ ==
"Chebyshev" && lambdaMax == -1.0) {
270 Teuchos::RCP<Ifpack_Chebyshev> chebyPrec = rcp_dynamic_cast<Ifpack_Chebyshev>(prec_);
271 if (chebyPrec != Teuchos::null) {
272 lambdaMax = chebyPrec->GetLambdaMax();
273 A_->SetMaxEigenvalueEstimate(lambdaMax);
274 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack)" <<
" = " << lambdaMax << std::endl;
276 TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == -1.0,
Exceptions::RuntimeError,
"MueLu::IfpackSmoother::Setup(): no maximum eigenvalue estimate");
279 this->GetOStream(
Statistics1) << description() << std::endl;
283 template <
class Node>
286 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
288 if (this->IsSetup() ==
true) {
289 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2moother::SetupAggregate(): Setup() has already been called" << std::endl;
290 this->GetOStream(
Warnings0) <<
"MueLu::IfpackSmoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
293 this->GetOStream(
Statistics0) <<
"IfpackSmoother: Using Aggregate Smoothing"<<std::endl;
295 RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates> >(currentLevel,
"Aggregates");
296 RCP<const LOMultiVector> vertex2AggId = aggregates->GetVertex2AggId();
297 ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
298 ArrayRCP<LO> dof_ids;
301 if(A_->GetFixedBlockSize() > 1) {
306 LO blocksize = (LO) A_->GetFixedBlockSize();
307 dof_ids.resize(aggregate_ids.size() * blocksize);
308 for(LO i=0; i<(LO)aggregate_ids.size(); i++) {
309 for(LO j=0; j<(LO)blocksize; j++)
310 dof_ids[i*blocksize+j] = aggregate_ids[i];
314 dof_ids = aggregate_ids;
317 paramList.set(
"partitioner: map", dof_ids.getRawPtr());
318 paramList.set(
"partitioner: type",
"user");
319 paramList.set(
"partitioner: overlap", 0);
320 paramList.set(
"partitioner: local parts", (
int)aggregates->GetNumAggregates());
322 paramList.set(
"partitioner: keep singletons",
true);
325 type_ =
"block relaxation stand-alone";
328 prec_ = rcp(factory.
Create(type_, &(*A), overlap_));
329 TEUCHOS_TEST_FOR_EXCEPTION(prec_.is_null(),
Exceptions::RuntimeError,
"Could not create an Ifpack preconditioner with type = \"" << type_ <<
"\"");
332 int rv = prec_->Compute();
333 TEUCHOS_TEST_FOR_EXCEPTION(rv,
Exceptions::RuntimeError,
"Ifpack preconditioner with type = \"" << type_ <<
"\" Compute() call failed.");
338 template <
class Node>
344 Teuchos::ParameterList paramList;
345 bool supportInitialGuess =
false;
346 if (type_ ==
"Chebyshev") {
347 paramList.set(
"chebyshev: zero starting solution", InitialGuessIsZero);
348 supportInitialGuess =
true;
350 }
else if (type_ ==
"point relaxation stand-alone") {
351 paramList.set(
"relaxation: zero starting solution", InitialGuessIsZero);
352 supportInitialGuess =
true;
355 SetPrecParameters(paramList);
358 if (InitialGuessIsZero || supportInitialGuess) {
362 prec_->ApplyInverse(epB, epX);
366 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
371 prec_->ApplyInverse(epB, epX);
373 X.update(1.0, *Correction, 1.0);
377 template <
class Node>
380 smoother->SetParameterList(this->GetParameterList());
381 return Teuchos::rcp_dynamic_cast<MueLu::SmootherPrototype<double, int, int, Node> >(smoother);
384 template <
class Node>
386 std::ostringstream out;
389 if (prec_ == Teuchos::null || this->GetVerbLevel() ==
InterfaceTest) {
391 out <<
"{type = " << type_ <<
"}";
393 out << prec_->Label();
398 template <
class Node>
403 out0 <<
"Prec. type: " << type_ << std::endl;
406 out0 <<
"Parameter list: " << std::endl;
407 Teuchos::OSTab tab2(out);
408 out << this->GetParameterList();
409 out0 <<
"Overlap: " << overlap_ << std::endl;
413 if (prec_ != Teuchos::null) {
414 Teuchos::OSTab tab2(out);
415 out << *prec_ << std::endl;
418 if (verbLevel &
Debug) {
421 <<
"RCP<A_>: " << A_ << std::endl
422 <<
"RCP<prec_>: " << prec_ << std::endl;
426 template <
class Node>
429 return Teuchos::OrdinalTraits<size_t>::invalid();
438#if defined(HAVE_MUELU_EPETRA)
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Ifpack_Preconditioner * Create(const std::string PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
virtual std::string description() const
Return a simple one-line description of this object.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
Class that encapsulates Ifpack smoothers.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
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.
std::string description() const
Return a simple one-line description of this object.
void Setup(Level ¤tLevel)
Set up the smoother.
RCP< SmootherPrototype > Copy() const
void DeclareInput(Level ¤tLevel) const
Input.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void SetupAggregate(Level ¤tLevel)
void SetParameterList(const Teuchos::ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
Class that holds all level-specific information.
RCP< Level > & GetPreviousLevel()
Previous level.
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual void SetParameterList(const Teuchos::ParameterList ¶mList)=0
Set parameters from a parameter list and return with default values.
void declareConstructionOutcome(bool fail, std::string msg)
bool IsSetup() const
Get the state of a smoother prototype.
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
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)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
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.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Statistics1
Print more statistics.
@ External
Print external lib objects.
@ Runtime0
One-liner description of what is happening.
@ Parameters0
Print class parameters.
@ Statistics0
Print statistics that do not involve significant additional computation.
@ Parameters1
Print class parameters (more parameters, more verbose)