Teko_BlockDiagonalInverseOp.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_BlockDiagonalInverseOp_hpp__
#define __Teko_BlockDiagonalInverseOp_hpp__

#include "Teko_Utilities.hpp"
#include "Teko_BlockImplicitLinearOp.hpp"

namespace Teko {

class BlockDiagonalInverseOp : public BlockImplicitLinearOp {
 public:
  BlockDiagonalInverseOp(BlockedLinearOp &A, const std::vector<LinearOp> &invDiag);

  BlockDiagonalInverseOp()                               = delete;
  BlockDiagonalInverseOp(const BlockDiagonalInverseOp &) = delete;



  virtual VectorSpace range() const { return productRange_; }

  virtual VectorSpace domain() const { return productDomain_; }

  virtual void implicitApply(const BlockedMultiVector &x, BlockedMultiVector &y,
                             const double alpha = 1.0, const double beta = 0.0) const;

  virtual void implicitApply(const Thyra::EOpTransp M_trans, const BlockedMultiVector &x,
                             BlockedMultiVector &y, const double alpha = 1.0,
                             const double beta = 0.0) const;

  virtual void describe(Teuchos::FancyOStream &out_arg,
                        const Teuchos::EVerbosityLevel verbLevel) const;

 protected:
  // fundamental operators to use
  std::vector<LinearOp> invDiag_;  

  Teuchos::RCP<const Thyra::ProductVectorSpaceBase<double> >
      productRange_;  
  Teuchos::RCP<const Thyra::ProductVectorSpaceBase<double> >
      productDomain_;  

  // scratch space...so we don't have to reallocate
  mutable BlockedMultiVector srcScrap_;
  mutable BlockedMultiVector dstScrap_;
  mutable bool allocated = false;
};

inline LinearOp createBlockDiagonalInverseOp(BlockedLinearOp &A,
                                             const std::vector<LinearOp> &invDiag) {
  return Teuchos::rcp(new BlockDiagonalInverseOp(A, invDiag));
}

inline LinearOp createBlockDiagonalInverseOp(BlockedLinearOp &A,
                                             const std::vector<LinearOp> &invDiag,
                                             const std::string &str) {
  Teuchos::RCP<Thyra::LinearOpBase<double> > result =
      Teuchos::rcp(new BlockDiagonalInverseOp(A, invDiag));
  result->setObjectLabel(str);

  return result;
}

}  // end namespace Teko

#endif