1 #include "Teko_MLLinearOp.hpp"
3 #include "Teko_EpetraOperatorWrapper.hpp"
4 #include "Teko_mlutils.hpp"
5 #include "Teko_PreconditionerLinearOp.hpp"
7 #include "ml_MultiLevelPreconditioner.h"
11 MLLinearOp::MLLinearOp(
const Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> &mlPrecOp)
12 : mlPrecOp_(mlPrecOp) {
13 extractConversionInformation(*mlPrecOp_);
16 void MLLinearOp::extractConversionInformation(ML_Epetra::MultiLevelPreconditioner &mlPrec) {
17 const ML *ml = mlPrec.GetML();
18 const ML_Smoother *preSmoother = ml->pre_smoother;
19 const ML_Smoother *postSmoother = ml->post_smoother;
22 const mlutils::SmootherData *smootherData;
24 smootherData = (
const mlutils::SmootherData *)ML_Get_MySmootherData(preSmoother);
25 else if (postSmoother != 0)
26 smootherData = (
const mlutils::SmootherData *)ML_Get_MySmootherData(preSmoother);
28 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
29 "MLLinearOp::extractConversionInformation pre and post smoother "
30 <<
"are both null, cannot build operator");
32 Amat_ = Teuchos::rcp_dynamic_cast<Epetra::EpetraOperatorWrapper>(smootherData->Amat);
35 mappingStrategy_ = Amat_->getMapStrategy();
38 BlockedLinearOp bloA = toBlockedLinearOp(Amat_->getThyraOp());
41 productRange_ = bloA->productDomain();
44 void MLLinearOp::implicitApply(
const BlockedMultiVector &x, BlockedMultiVector &y,
45 const double alpha,
const double beta)
const {
46 int columns = x->domain()->dim();
47 TEUCHOS_ASSERT(columns == y->domain()->dim());
50 if (eX_ == Teuchos::null || columns != eX_->NumVectors()) {
51 eX_ = Teuchos::rcp(
new Epetra_MultiVector(mlPrecOp_->OperatorDomainMap(), x->domain()->dim()));
52 eY_ = Teuchos::rcp(
new Epetra_MultiVector(mlPrecOp_->OperatorRangeMap(), y->domain()->dim()));
55 BlockedMultiVector yCopy;
63 mappingStrategy_->copyThyraIntoEpetra(x, *eX_);
66 mlPrecOp_->ApplyInverse(*eX_, *eY_);
69 mappingStrategy_->copyEpetraIntoThyra(*eY_, yCopy.ptr());
73 update(alpha, yCopy, beta, y);
74 else if (alpha != 1.0)
78 void MLLinearOp::describe(Teuchos::FancyOStream &out_arg,
79 const Teuchos::EVerbosityLevel verbLevel)
const {
80 out_arg <<
"MLLinearOp";
83 Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner()
const {
87 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner() {
91 Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> getMLPreconditioner(
92 const Teko::LinearOp &lo) {
93 Teko::LinearOp precOp = lo;
96 Teuchos::RCP<const Teko::PreconditionerLinearOp<double> > plo =
98 if (plo != Teuchos::null) precOp = plo->
getOperator();
101 Teuchos::RCP<const MLLinearOp> mlOp = Teuchos::rcp_dynamic_cast<
const MLLinearOp>(precOp);
103 TEUCHOS_TEST_FOR_EXCEPTION(
104 mlOp == Teuchos::null, std::runtime_error,
105 "Teko::getMLPreconditioner could not extract a MLLinearOp from the passed in argument");
107 return mlOp->getMLPreconditioner();
Teko::LinearOp getOperator() const
Get teko linear operator.
Class that wraps a PreconditionerBase object it makes it behave like a linear operator.