Teko_MLLinearOp.cpp

#include "Teko_MLLinearOp.hpp"

#include "Teko_EpetraOperatorWrapper.hpp"
#include "Teko_mlutils.hpp"
#include "Teko_PreconditionerLinearOp.hpp"

#include "ml_MultiLevelPreconditioner.h"

namespace Teko {

MLLinearOp::MLLinearOp(const Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> &mlPrecOp)
    : mlPrecOp_(mlPrecOp) {
  extractConversionInformation(*mlPrecOp_);
}

void MLLinearOp::extractConversionInformation(ML_Epetra::MultiLevelPreconditioner &mlPrec) {
  const ML *ml                    = mlPrec.GetML();
  const ML_Smoother *preSmoother  = ml->pre_smoother;
  const ML_Smoother *postSmoother = ml->post_smoother;

  // grab data object
  const mlutils::SmootherData *smootherData;
  if (preSmoother != 0)
    smootherData = (const mlutils::SmootherData *)ML_Get_MySmootherData(preSmoother);
  else if (postSmoother != 0)
    smootherData = (const mlutils::SmootherData *)ML_Get_MySmootherData(preSmoother);
  else
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
                               "MLLinearOp::extractConversionInformation pre and post smoother "
                                   << "are both null, cannot build operator");

  Amat_ = Teuchos::rcp_dynamic_cast<Epetra::EpetraOperatorWrapper>(smootherData->Amat);

  // for doing Epetra_Vector -> Thyra vector conversion
  mappingStrategy_ = Amat_->getMapStrategy();

  // to construct vectorspace for this object
  BlockedLinearOp bloA = toBlockedLinearOp(Amat_->getThyraOp());
  productDomain_ =
      bloA->productRange();  // this operator is the inverse of bloA...hence swap range and domain
  productRange_ = bloA->productDomain();
}

void MLLinearOp::implicitApply(const BlockedMultiVector &x, BlockedMultiVector &y,
                               const double alpha, const double beta) const {
  int columns = x->domain()->dim();
  TEUCHOS_ASSERT(columns == y->domain()->dim());  // check for equal columns

  // (re)allocate vectors conditinolly if required size changed
  if (eX_ == Teuchos::null || columns != eX_->NumVectors()) {
    eX_ = Teuchos::rcp(new Epetra_MultiVector(mlPrecOp_->OperatorDomainMap(), x->domain()->dim()));
    eY_ = Teuchos::rcp(new Epetra_MultiVector(mlPrecOp_->OperatorRangeMap(), y->domain()->dim()));
  }

  BlockedMultiVector yCopy;
  if (beta != 0)
    yCopy = deepcopy(y);
  else
    yCopy = y;

  // initialize Epetra vectors
  eY_->PutScalar(0.0);
  mappingStrategy_->copyThyraIntoEpetra(x, *eX_);

  // run multigrid
  mlPrecOp_->ApplyInverse(*eX_, *eY_);

  // fill thyra vectors
  mappingStrategy_->copyEpetraIntoThyra(*eY_, yCopy.ptr());

  // scale result by alpha
  if (beta != 0)
    update(alpha, yCopy, beta, y);  // y = alpha * yCopy + beta * y
  else if (alpha != 1.0)
    scale(alpha, y);  // y = alpha * y
}

void MLLinearOp::describe(Teuchos::FancyOStream &out_arg,
                          const Teuchos::EVerbosityLevel verbLevel) const {
  out_arg << "MLLinearOp";
}

Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner() const {
  return mlPrecOp_;
}

Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> MLLinearOp::getMLPreconditioner() {
  return mlPrecOp_;
}

Teuchos::RCP<const ML_Epetra::MultiLevelPreconditioner> getMLPreconditioner(
    const Teko::LinearOp &lo) {
  Teko::LinearOp precOp = lo;

  // try to pull it of a preconditioner linear op
  Teuchos::RCP<const Teko::PreconditionerLinearOp<double> > plo =
      Teuchos::rcp_dynamic_cast<const Teko::PreconditionerLinearOp<double> >(lo);
  if (plo != Teuchos::null) precOp = plo->getOperator();

  // try to extract the ML operator
  Teuchos::RCP<const MLLinearOp> mlOp = Teuchos::rcp_dynamic_cast<const MLLinearOp>(precOp);

  TEUCHOS_TEST_FOR_EXCEPTION(
      mlOp == Teuchos::null, std::runtime_error,
      "Teko::getMLPreconditioner could not extract a MLLinearOp from the passed in argument");

  return mlOp->getMLPreconditioner();
}

}  // namespace Teko