Teko_NeumannSeriesPreconditionerFactory.hpp

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef __Teko_NeumannSeriesPreconditionerFactory_hpp__
#define __Teko_NeumannSeriesPreconditionerFactory_hpp__
 
#include "Teko_NeumannSeriesPreconditionerFactoryDecl.hpp"
 
#include "Thyra_DefaultPreconditioner.hpp"
#include "Thyra_DefaultPreconditioner.hpp"
#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
#include "Thyra_DefaultAddedLinearOp.hpp"
#include "Thyra_DefaultMultipliedLinearOp.hpp"
#include "Thyra_DefaultIdentityLinearOp.hpp"
 
#include "Teuchos_Array.hpp"
#include "Teuchos_StandardParameterEntryValidators.hpp"
 
namespace Teko {
 
using Teuchos::RCP;
using Teuchos::rcp;
 
static RCP<Teuchos::StringToIntegralParameterEntryValidator<Teko::DiagonalType> > scalingTypeVdtor;
 
template <typename ScalarT>
NeumannSeriesPreconditionerFactory<ScalarT>::NeumannSeriesPreconditionerFactory()
    : numberOfTerms_(1), scalingType_(Teko::NotDiag) {}

template <typename ScalarT>
bool NeumannSeriesPreconditionerFactory<ScalarT>::isCompatible(
    const Thyra::LinearOpSourceBase<ScalarT> & /* fwdOpSrc */) const {
  return true;
}

template <typename ScalarT>
RCP<Thyra::PreconditionerBase<ScalarT> > NeumannSeriesPreconditionerFactory<ScalarT>::createPrec()
    const {
  return rcp(new Thyra::DefaultPreconditioner<ScalarT>());
}

template <typename ScalarT>
void NeumannSeriesPreconditionerFactory<ScalarT>::initializePrec(
    const RCP<const Thyra::LinearOpSourceBase<ScalarT> > &fwdOpSrc,
    Thyra::PreconditionerBase<ScalarT> *prec,
    const Thyra::ESupportSolveUse /* supportSolveUse */) const {
  using Thyra::add;
  using Thyra::multiply;
  using Thyra::scale;
 
  RCP<const Thyra::LinearOpBase<ScalarT> > M;  // left-preconditioner
  RCP<const Thyra::LinearOpBase<ScalarT> > A = fwdOpSrc->getOp();
  if (scalingType_ != Teko::NotDiag) {
    M = Teko::getInvDiagonalOp(A, scalingType_);
    A = Thyra::multiply(M, A);
  }
 
  RCP<const Thyra::LinearOpBase<ScalarT> > id   = Thyra::identity<ScalarT>(A->range());  // I
  RCP<const Thyra::LinearOpBase<ScalarT> > idMA = add(id, scale(-1.0, A));               // I - A
 
  RCP<const Thyra::LinearOpBase<ScalarT> > precOp;
  if (numberOfTerms_ == 1) {
    // no terms requested, just return identity matrix
    precOp = id;
  } else {
    int iters = numberOfTerms_ - 1;
    // use Horner's method to computed higher order polynomials
    precOp = add(scale(2.0, id), scale(-1.0, A));                              // I + (I - A)
    for (int i = 0; i < iters; i++) precOp = add(id, multiply(idMA, precOp));  // I + (I - A) * p
  }
 
  // multiply by the preconditioner if it exists
  if (M != Teuchos::null) precOp = Thyra::multiply(precOp, M);
 
  // must first cast that to be initialized
  Thyra::DefaultPreconditioner<ScalarT> &dPrec =
      Teuchos::dyn_cast<Thyra::DefaultPreconditioner<ScalarT> >(*prec);
 
  // this const-cast is unfortunate...needs to be fixed (larger than it seems!) ECC FIXME!
  dPrec.initializeUnspecified(Teuchos::rcp_const_cast<Thyra::LinearOpBase<ScalarT> >(precOp));
}

template <typename ScalarT>
void NeumannSeriesPreconditionerFactory<ScalarT>::uninitializePrec(
    Thyra::PreconditionerBase<ScalarT> *prec,
    RCP<const Thyra::LinearOpSourceBase<ScalarT> > * /* fwdOpSrc */,
    Thyra::ESupportSolveUse * /* supportSolveUse */) const {
  Thyra::DefaultPreconditioner<ScalarT> &dPrec =
      Teuchos::dyn_cast<Thyra::DefaultPreconditioner<ScalarT> >(*prec);
 
  // wipe out any old preconditioner
  dPrec.uninitialize();
}
 
// for ParameterListAcceptor

template <typename ScalarT>
void NeumannSeriesPreconditionerFactory<ScalarT>::setParameterList(
    const RCP<Teuchos::ParameterList> &paramList) {
  TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
 
  // check the parameters are correct
  paramList->validateParametersAndSetDefaults(*getValidParameters(), 0);
 
  // store the parameter list
  paramList_ = paramList;
 
  numberOfTerms_ = paramList_->get<int>("Number of Terms");
 
  // get the scaling type
  scalingType_                         = Teko::NotDiag;
  const Teuchos::ParameterEntry *entry = paramList_->getEntryPtr("Scaling Type");
  if (entry != NULL) scalingType_ = scalingTypeVdtor->getIntegralValue(*entry);
}

template <typename ScalarT>
RCP<const Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getValidParameters()
    const {
  static RCP<Teuchos::ParameterList> validPL;
 
  // only fill valid parameter list once
  if (validPL == Teuchos::null) {
    RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
 
    // build the validator for scaling type
    scalingTypeVdtor = Teuchos::stringToIntegralParameterEntryValidator<DiagonalType>(
        Teuchos::tuple<std::string>("Diagonal", "Lumped", "AbsRowSum", "None"),
        Teuchos::tuple<Teko::DiagonalType>(Teko::Diagonal, Teko::Lumped, Teko::AbsRowSum,
                                           Teko::NotDiag),
        "Scaling Type");
 
    pl->set<int>("Number of Terms", 1,
                 "The number of terms to use in the Neumann series expansion.");
    pl->set("Scaling Type", "None", "The number of terms to use in the Neumann series expansion.",
            scalingTypeVdtor);
 
    validPL = pl;
  }
 
  return validPL;
}

template <typename ScalarT>
RCP<Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::unsetParameterList() {
  Teuchos::RCP<Teuchos::ParameterList> oldList = paramList_;
  paramList_                                   = Teuchos::null;
  return oldList;
}

template <typename ScalarT>
RCP<const Teuchos::ParameterList> NeumannSeriesPreconditionerFactory<ScalarT>::getParameterList()
    const {
  return paramList_;
}

template <typename ScalarT>
RCP<Teuchos::ParameterList>
NeumannSeriesPreconditionerFactory<ScalarT>::getNonconstParameterList() {
  return paramList_;
}
 
template <typename ScalarT>
std::string NeumannSeriesPreconditionerFactory<ScalarT>::description() const {
  std::ostringstream oss;
  oss << "Teko::NeumannSeriesPreconditionerFactory";
  return oss.str();
}
 
}  // end namespace Teko
 
#endif