Teko_MLPreconditionerFactory.cpp

#include "Teko_MLPreconditionerFactory.hpp"

#include "Teko_MLLinearOp.hpp"

#include "ml_include.h"
#include "ml_MultiLevelPreconditioner.h"
#include "ml_epetra_utils.h"
#include "ml_RowMatrix.h"

#include "Thyra_get_Epetra_Operator.hpp"

namespace Teko {

// MLPreconditionerState

MLPreconditionerState::MLPreconditionerState() : isFilled_(false) {}

void MLPreconditionerState::setMLComm(ML_Comm *comm) {
  mlComm_ = Teuchos::rcpWithDealloc(comm, Teuchos::deallocFunctorDelete<ML_Comm>(cleanup_ML_Comm));
}

void MLPreconditionerState::setMLOperator(ML_Operator *op) {
  mlOp_ =
      Teuchos::rcpWithDealloc(op, Teuchos::deallocFunctorDelete<ML_Operator>(cleanup_ML_Operator));
}

void MLPreconditionerState::setIsFilled(bool value) { isFilled_ = value; }

bool MLPreconditionerState::isFilled() const { return isFilled_; }

void MLPreconditionerState::cleanup_ML_Comm(ML_Comm *mlComm) { ML_Comm_Destroy(&mlComm); }

void MLPreconditionerState::cleanup_ML_Operator(ML_Operator *mlOp) { ML_Operator_Destroy(&mlOp); }

void MLPreconditionerState::setAggregationMatrices(const std::vector<Epetra_RowMatrix *> &diags) {
  diagonalOps_ = diags;
}

Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> MLPreconditionerState::constructMLPreconditioner(
    const Teuchos::ParameterList &mainList,
    const std::vector<Teuchos::RCP<const Teuchos::ParameterList> > &coarseningParams) {
  TEUCHOS_ASSERT(isFilled());
  std::vector<Teuchos::ParameterList> cpls(coarseningParams.size());
  for (std::size_t i = 0; i < coarseningParams.size(); i++) cpls[i] = *coarseningParams[i];

  mlPreconditioner_ = rcp(new ML_Epetra::MultiLevelPreconditioner(
      &*mlOp_, mainList, &diagonalOps_[0], &cpls[0], diagonalOps_.size()));

  return mlPreconditioner_;
}

// MLPreconditionerFactory

MLPreconditionerFactory::MLPreconditionerFactory() {}

LinearOp MLPreconditionerFactory::buildPreconditionerOperator(
    BlockedLinearOp &blo, BlockPreconditionerState &state) const {
  MLPreconditionerState &mlState = Teuchos::dyn_cast<MLPreconditionerState>(state);

  if (not mlState.isFilled()) fillMLPreconditionerState(blo, mlState);
  TEUCHOS_ASSERT(mlState.isFilled());

  Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> mlPrecOp =
      mlState.constructMLPreconditioner(mainParams_, coarseningParams_);

  // return Thyra::epetraLinearOp(mlPrecOp);
  return Teuchos::rcp(new MLLinearOp(mlPrecOp));
}

Teuchos::RCP<PreconditionerState> MLPreconditionerFactory::buildPreconditionerState() const {
  return Teuchos::rcp(new MLPreconditionerState());
}

void MLPreconditionerFactory::fillMLPreconditionerState(const BlockedLinearOp &blo,
                                                        MLPreconditionerState &mlState) const {
  TEUCHOS_ASSERT(not mlState.isFilled());

  EpetraExt::CrsMatrix_SolverMap RowMatrixColMapTrans;

  int rowCnt = blockRowCount(blo);
  int colCnt = blockRowCount(blo);

  // construct comm object...add to state
  ML_Comm *mlComm;
  ML_Comm_Create(&mlComm);
  mlState.setMLComm(mlComm);

  ML_Operator *mlBlkMat = ML_Operator_Create(mlComm);
  ML_Operator_BlkMatInit(mlBlkMat, mlComm, rowCnt, colCnt, ML_DESTROY_EVERYTHING);

  std::vector<Epetra_RowMatrix *> aggMats;
  for (int r = 0; r < rowCnt; r++) {
    for (int c = 0; c < colCnt; c++) {
      ML_Operator *tmp       = 0;
      Epetra_RowMatrix *EMat = 0;
      std::stringstream ss;

      ss << "fine_" << r << "_" << c;

      // extract crs matrix
      Teuchos::RCP<Epetra_CrsMatrix> crsMat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(
          Thyra::get_Epetra_Operator(*blo->getNonconstBlock(r, c)));
      EMat = ML_Epetra::ModifyEpetraMatrixColMap(*crsMat, RowMatrixColMapTrans, ss.str().c_str());
      if (r == c)  // setup diagonal scaling matrices
        aggMats.push_back(EMat);

      // extract ml sub operator
      tmp = ML_Operator_Create(mlComm);
      ML_Operator_WrapEpetraMatrix(EMat, tmp);

      // add it to the block
      ML_Operator_BlkMatInsert(mlBlkMat, tmp, r, c);
    }
  }
  ML_Operator_BlkMatFinalize(mlBlkMat);

  // finish setting up state object
  mlState.setMLOperator(mlBlkMat);

  // now set aggregation matrices
  mlState.setAggregationMatrices(aggMats);

  mlState.setIsFilled(true);  // register state object as filled
}

void MLPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList &settings) {
  using Teuchos::RCP;
  using Teuchos::rcp;

  Teko_DEBUG_SCOPE("MLPreconditionerFactory::initializeFromParameterList", 10);

  // clear initial state
  coarseningParams_.clear();
  blockRowCount_ = 0;

  blockRowCount_ = settings.get<int>("Block Row Count");

  // read in main parameter list: with smoothing information
  mainParams_ = settings.sublist("Smoothing Parameters");
  mainParams_.set<Teuchos::RCP<const Teko::InverseLibrary> >("smoother: teko inverse library",
                                                             getInverseLibrary());

  // read in aggregation sub lists
  const Teuchos::ParameterList &aggSublist = settings.sublist("Block Aggregation");

  for (int block = 0; block < blockRowCount_; block++) {
    // write out sub list name: "Block #"
    std::stringstream ss;
    ss << "Block " << block;
    std::string sublistName = ss.str();

    // grab sublist
    RCP<Teuchos::ParameterList> userSpec =
        rcp(new Teuchos::ParameterList(aggSublist.sublist(sublistName)));
    coarseningParams_.push_back(userSpec);
  }
}

}  // end namespace Teko