Thyra_ScaledIdentityLinearOpWithSolve.hpp

Go to the documentation of this file.

// @HEADER
// ****************************************************************************
//                Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER

#ifndef Thyra_ScaledIdentityLinearOpWithSolve_hpp
#define Thyra_ScaledIdentityLinearOpWithSolve_hpp

#include "Thyra_LinearOpWithSolveBase.hpp"
#include "Thyra_MultiVectorStdOps.hpp"

namespace Thyra {

template <class Scalar>
class ScaledIdentityLinearOpWithSolve
  : virtual public LinearOpWithSolveBase<Scalar> {
 public:

  ScaledIdentityLinearOpWithSolve() {}

  void initialize(const RCP<const VectorSpaceBase<Scalar> >& space,
                  const Scalar& s)
  {
    validateInitialize(space);
    space_ = space;
    s_     = s;
  }

  void uninitialize() { space_ = Teuchos::null; }

  RCP<const VectorSpaceBase<Scalar> > space() const { return space_; }
  Scalar scale() const { return s_; }
  void setScale(const Scalar& s) { s_ = s; }



  RCP<const VectorSpaceBase<Scalar> > range() const { return space_; }

  RCP<const VectorSpaceBase<Scalar> > domain() const { return space_; }

  RCP<const LinearOpBase<Scalar> > clone() const
  {
    RCP<ScaledIdentityLinearOpWithSolve<Scalar> > op =
        rcp(new ScaledIdentityLinearOpWithSolve<Scalar>());
    op->initialize(space_, s_);
    return op;
  }

 protected:
  bool opSupportedImpl(EOpTransp /* M_trans */) const { return true; }

  void applyImpl(const EOpTransp M_trans, const MultiVectorBase<Scalar>& X,
                 const Ptr<MultiVectorBase<Scalar> >& Y, const Scalar alpha,
                 const Scalar beta) const
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    Thyra::scale(beta, Y);
    if (M_trans == CONJ || M_trans == CONJTRANS)
      V_StVpV(Y, ST::conjugate(s_) * alpha, X, *Y);
    else
      V_StVpV(Y, s_ * alpha, X, *Y);
  }

  bool solveSupportsImpl(EOpTransp /* M_trans */) const { return true; }

  bool solveSupportsNewImpl(
      EOpTransp /* M_trans */,
      const Ptr<const SolveCriteria<Scalar> > /* solveCriteria */) const
  {
    return true;
  }

  bool solveSupportsSolveMeasureTypeImpl(
      EOpTransp /* M_trans */,
      const SolveMeasureType& /* solveMeasureType */) const
  {
    return true;
  }

  SolveStatus<Scalar> solveImpl(
      const EOpTransp M_trans, const MultiVectorBase<Scalar>& B,
      const Ptr<MultiVectorBase<Scalar> >& X,
      const Ptr<const SolveCriteria<Scalar> > /* solveCriteria */) const
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;
    assign(X, ST::zero());
    if (M_trans == CONJ || M_trans == CONJTRANS)
      V_StVpV(X, ST::one() / ST::conjugate(s_), B, *X);
    else
      V_StVpV(X, ST::one() / s_, B, *X);
    SolveStatus<Scalar> solveStatus;
    solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
    return solveStatus;
  }

 private:
  // //////////////////////////////
  // Private data members

  RCP<const VectorSpaceBase<Scalar> > space_;
  Scalar s_;

  // //////////////////////////////
  // Private member functions

  static void validateInitialize(
      const RCP<const VectorSpaceBase<Scalar> >& space)
  {
#ifdef TEUCHOS_DEBUG
    TEUCHOS_TEST_FOR_EXCEPT(is_null(space));
#else
    (void)space;
#endif
  }
};

template <class Scalar>
RCP<ScaledIdentityLinearOpWithSolve<Scalar> > scaledIdentity()
{
  return Teuchos::rcp(new ScaledIdentityLinearOpWithSolve<Scalar>());
}

template <class Scalar>
RCP<ScaledIdentityLinearOpWithSolve<Scalar> > scaledIdentity(
    const RCP<const VectorSpaceBase<Scalar> >& space, const Scalar& s)
{
  RCP<ScaledIdentityLinearOpWithSolve<Scalar> > op =
      Teuchos::rcp(new ScaledIdentityLinearOpWithSolve<Scalar>());
  op->initialize(space, s);
  return op;
}

}  // end namespace Thyra

#endif