Teko Version of the Day
Teko_MLLinearOp.cpp
00001 #include "Teko_MLLinearOp.hpp"
00002 
00003 #include "Teko_EpetraOperatorWrapper.hpp"
00004 #include "Teko_mlutils.hpp"
00005 #include "Teko_PreconditionerLinearOp.hpp"
00006 
00007 #include "ml_MultiLevelPreconditioner.h"
00008 
00009 namespace Teko {
00010 
00011 MLLinearOp::MLLinearOp(const Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> & mlPrecOp)
00012    : mlPrecOp_(mlPrecOp)
00013 {
00014    extractConversionInformation(*mlPrecOp_);
00015 }
00016 
00017 void MLLinearOp::extractConversionInformation(ML_Epetra::MultiLevelPreconditioner & mlPrec)
00018 {
00019    const ML * ml = mlPrec.GetML();
00020    const ML_Smoother * preSmoother = ml->pre_smoother;
00021    const ML_Smoother * postSmoother = ml->post_smoother;
00022 
00023    // grab data object
00024    const mlutils::SmootherData * smootherData;
00025    if(preSmoother!=0)
00026       smootherData = (const mlutils::SmootherData *) ML_Get_MySmootherData(preSmoother);
00027    else if(postSmoother!=0)
00028       smootherData = (const mlutils::SmootherData *) ML_Get_MySmootherData(preSmoother);
00029    else
00030       TEST_FOR_EXCEPTION(true,std::runtime_error,
00031                               "MLLinearOp::extractConversionInformation pre and post smoother " <<
00032                               "are both null, cannot build operator");
00033 
00034    Amat_ = Teuchos::rcp_dynamic_cast<Epetra::EpetraOperatorWrapper>(smootherData->Amat);
00035 
00036    // for doing Epetra_Vector -> Thyra vector conversion
00037    mappingStrategy_ = Amat_->getMapStrategy();
00038 
00039    // to construct vectorspace for this object
00040    BlockedLinearOp bloA = toBlockedLinearOp(Amat_->getThyraOp());
00041    productDomain_ = bloA->productRange(); // this operator is the inverse of bloA...hence swap range and domain
00042    productRange_ = bloA->productDomain();
00043 
00044 }
00045 
00046 void MLLinearOp::implicitApply(const BlockedMultiVector & x, BlockedMultiVector & y,
00047                                const double alpha, const double beta) const
00048 {
00049    int columns = x->domain()->dim();
00050    TEUCHOS_ASSERT(columns==y->domain()->dim()); // check for equal columns
00051 
00052    // (re)allocate vectors conditinolly if required size changed
00053    if(eX_==Teuchos::null || columns!=eX_->NumVectors()) {
00054       eX_ = Teuchos::rcp(new Epetra_MultiVector(mlPrecOp_->OperatorDomainMap(),x->domain()->dim()));
00055       eY_ = Teuchos::rcp(new Epetra_MultiVector(mlPrecOp_->OperatorRangeMap(),y->domain()->dim()));
00056    }
00057 
00058    BlockedMultiVector yCopy;
00059    if(beta!=0)
00060       yCopy = deepcopy(y);
00061    else
00062       yCopy = y;
00063 
00064    // initialize Epetra vectors
00065    eY_->PutScalar(0.0);
00066    mappingStrategy_->copyThyraIntoEpetra(x, *eX_);
00067 
00068    // run multigrid
00069    mlPrecOp_->ApplyInverse(*eX_,*eY_);
00070 
00071    // fill thyra vectors
00072    mappingStrategy_->copyEpetraIntoThyra(*eY_, yCopy.ptr());
00073 
00074    // scale result by alpha
00075    if(beta!=0)
00076       update(alpha,yCopy,beta,y); // y = alpha * yCopy + beta * y 
00077    else if(alpha!=1.0)
00078       scale(alpha,y); // y = alpha * y
00079 }
00080 
00081 void MLLinearOp::describe(Teuchos::FancyOStream & out_arg,
00082                           const Teuchos::EVerbosityLevel verbLevel) const
00083 {
00084    out_arg << "MLLinearOp";
00085 }
00086 
00087 Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner() const
00088 {
00089    return mlPrecOp_;
00090 }
00091 
00092 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner()
00093 {
00094    return mlPrecOp_;
00095 }
00096 
00097 Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> getMLPreconditioner(const Teko::LinearOp & lo)
00098 {
00099    Teko::LinearOp precOp = lo;
00100 
00101    // try to pull it of a preconditioner linear op
00102    Teuchos::RCP<const Teko::PreconditionerLinearOp<double> > plo 
00103          = Teuchos::rcp_dynamic_cast<const Teko::PreconditionerLinearOp<double> >(lo);
00104    if(plo!=Teuchos::null)
00105       precOp = plo->getOperator();
00106    
00107    // try to extract the ML operator
00108    Teuchos::RCP<const MLLinearOp> mlOp = Teuchos::rcp_dynamic_cast<const MLLinearOp>(precOp);
00109 
00110    TEST_FOR_EXCEPTION(mlOp==Teuchos::null,std::runtime_error,
00111                       "Teko::getMLPreconditioner could not extract a MLLinearOp from the passed in argument");
00112 
00113    return mlOp->getMLPreconditioner();
00114 }
00115 
00116 }
 All Classes Files Functions Variables