Teko Version of the Day
Loading...
Searching...
No Matches
Teko_MLLinearOp.cpp
1#include "Teko_MLLinearOp.hpp"
2
3#include "Teko_EpetraOperatorWrapper.hpp"
4#include "Teko_mlutils.hpp"
5#include "Teko_PreconditionerLinearOp.hpp"
6
7#include "ml_MultiLevelPreconditioner.h"
8
9namespace Teko {
10
11MLLinearOp::MLLinearOp(const Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> & mlPrecOp)
12 : mlPrecOp_(mlPrecOp)
13{
14 extractConversionInformation(*mlPrecOp_);
15}
16
17void MLLinearOp::extractConversionInformation(ML_Epetra::MultiLevelPreconditioner & mlPrec)
18{
19 const ML * ml = mlPrec.GetML();
20 const ML_Smoother * preSmoother = ml->pre_smoother;
21 const ML_Smoother * postSmoother = ml->post_smoother;
22
23 // grab data object
24 const mlutils::SmootherData * smootherData;
25 if(preSmoother!=0)
26 smootherData = (const mlutils::SmootherData *) ML_Get_MySmootherData(preSmoother);
27 else if(postSmoother!=0)
28 smootherData = (const mlutils::SmootherData *) ML_Get_MySmootherData(preSmoother);
29 else
30 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
31 "MLLinearOp::extractConversionInformation pre and post smoother " <<
32 "are both null, cannot build operator");
33
34 Amat_ = Teuchos::rcp_dynamic_cast<Epetra::EpetraOperatorWrapper>(smootherData->Amat);
35
36 // for doing Epetra_Vector -> Thyra vector conversion
37 mappingStrategy_ = Amat_->getMapStrategy();
38
39 // to construct vectorspace for this object
40 BlockedLinearOp bloA = toBlockedLinearOp(Amat_->getThyraOp());
41 productDomain_ = bloA->productRange(); // this operator is the inverse of bloA...hence swap range and domain
42 productRange_ = bloA->productDomain();
43
44}
45
46void MLLinearOp::implicitApply(const BlockedMultiVector & x, BlockedMultiVector & y,
47 const double alpha, const double beta) const
48{
49 int columns = x->domain()->dim();
50 TEUCHOS_ASSERT(columns==y->domain()->dim()); // check for equal columns
51
52 // (re)allocate vectors conditinolly if required size changed
53 if(eX_==Teuchos::null || columns!=eX_->NumVectors()) {
54 eX_ = Teuchos::rcp(new Epetra_MultiVector(mlPrecOp_->OperatorDomainMap(),x->domain()->dim()));
55 eY_ = Teuchos::rcp(new Epetra_MultiVector(mlPrecOp_->OperatorRangeMap(),y->domain()->dim()));
56 }
57
58 BlockedMultiVector yCopy;
59 if(beta!=0)
60 yCopy = deepcopy(y);
61 else
62 yCopy = y;
63
64 // initialize Epetra vectors
65 eY_->PutScalar(0.0);
66 mappingStrategy_->copyThyraIntoEpetra(x, *eX_);
67
68 // run multigrid
69 mlPrecOp_->ApplyInverse(*eX_,*eY_);
70
71 // fill thyra vectors
72 mappingStrategy_->copyEpetraIntoThyra(*eY_, yCopy.ptr());
73
74 // scale result by alpha
75 if(beta!=0)
76 update(alpha,yCopy,beta,y); // y = alpha * yCopy + beta * y
77 else if(alpha!=1.0)
78 scale(alpha,y); // y = alpha * y
79}
80
81void MLLinearOp::describe(Teuchos::FancyOStream & out_arg,
82 const Teuchos::EVerbosityLevel verbLevel) const
83{
84 out_arg << "MLLinearOp";
85}
86
87Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner() const
88{
89 return mlPrecOp_;
90}
91
92Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner()
93{
94 return mlPrecOp_;
95}
96
97Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> getMLPreconditioner(const Teko::LinearOp & lo)
98{
99 Teko::LinearOp precOp = lo;
100
101 // try to pull it of a preconditioner linear op
102 Teuchos::RCP<const Teko::PreconditionerLinearOp<double> > plo
103 = Teuchos::rcp_dynamic_cast<const Teko::PreconditionerLinearOp<double> >(lo);
104 if(plo!=Teuchos::null)
105 precOp = plo->getOperator();
106
107 // try to extract the ML operator
108 Teuchos::RCP<const MLLinearOp> mlOp = Teuchos::rcp_dynamic_cast<const MLLinearOp>(precOp);
109
110 TEUCHOS_TEST_FOR_EXCEPTION(mlOp==Teuchos::null,std::runtime_error,
111 "Teko::getMLPreconditioner could not extract a MLLinearOp from the passed in argument");
112
113 return mlOp->getMLPreconditioner();
114}
115
116}
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
MultiVector deepcopy(const MultiVector &v)
Perform a deep copy of the vector.
BlockedLinearOp toBlockedLinearOp(LinearOp &clo)
Converts a LinearOp to a BlockedLinearOp.