Teko_TpetraOperatorWrapper.cpp

// @HEADER
//
// ***********************************************************************
//
//      Teko: A package for block and physics based preconditioning
//                  Copyright 2010 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

#include "Teko_TpetraOperatorWrapper.hpp"
#include "Thyra_SpmdVectorBase.hpp"
#include "Thyra_MultiVectorStdOps.hpp"
#include "Teuchos_DefaultSerialComm.hpp"
#include "Thyra_TpetraLinearOp.hpp"
#include "Tpetra_Vector.hpp"
#include "Thyra_TpetraThyraWrappers.hpp"

// #include "Thyra_LinearOperator.hpp"
#include "Thyra_BlockedLinearOpBase.hpp"
#include "Thyra_ProductVectorSpaceBase.hpp"
#include "Thyra_TpetraLinearOp.hpp"

#include "Teko_TpetraThyraConverter.hpp"
#include "Teuchos_Ptr.hpp"

namespace Teko {
namespace TpetraHelpers {

using namespace Teuchos;
using namespace Thyra;

DefaultMappingStrategy::DefaultMappingStrategy(
    const RCP<const Thyra::LinearOpBase<double> >& thyraOp,
    const Teuchos::Comm<Thyra::Ordinal>& comm) {
  RCP<Teuchos::Comm<Thyra::Ordinal> > newComm = comm.duplicate();

  // extract vector spaces from linear operator
  domainSpace_ = thyraOp->domain();
  rangeSpace_  = thyraOp->range();

  domainMap_ = Teko::TpetraHelpers::thyraVSToTpetraMap(*domainSpace_, newComm);
  rangeMap_  = Teko::TpetraHelpers::thyraVSToTpetraMap(*rangeSpace_, newComm);
}

void DefaultMappingStrategy::copyTpetraIntoThyra(
    const Tpetra::MultiVector<ST, LO, GO, NT>& x,
    const Ptr<Thyra::MultiVectorBase<ST> >& thyraVec) const {
  Teko::TpetraHelpers::blockTpetraToThyra(x, thyraVec);
}

void DefaultMappingStrategy::copyThyraIntoTpetra(
    const RCP<const Thyra::MultiVectorBase<ST> >& thyraVec,
    Tpetra::MultiVector<ST, LO, GO, NT>& v) const {
  Teko::TpetraHelpers::blockThyraToTpetra(thyraVec, v);
}

TpetraOperatorWrapper::TpetraOperatorWrapper() {
  useTranspose_ = false;
  mapStrategy_  = Teuchos::null;
  thyraOp_      = Teuchos::null;
  comm_         = Teuchos::null;
  label_        = Teuchos::null;
}

TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<ST> >& thyraOp) {
  SetOperator(thyraOp);
}

TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const Thyra::LinearOpBase<ST> >& thyraOp,
                                             const RCP<const MappingStrategy>& mapStrategy)
    : mapStrategy_(mapStrategy) {
  SetOperator(thyraOp);
}

TpetraOperatorWrapper::TpetraOperatorWrapper(const RCP<const MappingStrategy>& mapStrategy)
    : mapStrategy_(mapStrategy) {
  useTranspose_ = false;
  thyraOp_      = Teuchos::null;
  comm_         = Teuchos::null;
  label_        = Teuchos::null;
}

void TpetraOperatorWrapper::SetOperator(const RCP<const LinearOpBase<ST> >& thyraOp,
                                        bool buildMap) {
  useTranspose_ = false;
  thyraOp_      = thyraOp;
  comm_         = getThyraComm(*thyraOp);
  label_        = thyraOp_->description();
  if (mapStrategy_ == Teuchos::null && buildMap)
    mapStrategy_ = Teuchos::rcp(new DefaultMappingStrategy(thyraOp, *comm_));
}

double TpetraOperatorWrapper::NormInf() const {
  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
                             "TpetraOperatorWrapper::NormInf not implemated");
  return 1.0;
}

void TpetraOperatorWrapper::apply(const Tpetra::MultiVector<ST, LO, GO, NT>& X,
                                  Tpetra::MultiVector<ST, LO, GO, NT>& Y,
                                  Teuchos::ETransp /* mode */, ST alpha, ST beta) const {
  if (!useTranspose_) {
    // allocate space for each vector
    RCP<Thyra::MultiVectorBase<ST> > tX;
    RCP<Thyra::MultiVectorBase<ST> > tY;

    tX = Thyra::createMembers(thyraOp_->domain(), X.getNumVectors());
    tY = Thyra::createMembers(thyraOp_->range(), X.getNumVectors());

    Thyra::assign(tX.ptr(), 0.0);
    Thyra::assign(tY.ptr(), 0.0);

    // copy epetra X into thyra X
    mapStrategy_->copyTpetraIntoThyra(X, tX.ptr());
    mapStrategy_->copyTpetraIntoThyra(
        Y, tY.ptr());  // if this matrix isn't block square, this probably won't work!

    // perform matrix vector multiplication
    thyraOp_->apply(Thyra::NOTRANS, *tX, tY.ptr(), alpha, beta);

    // copy thyra Y into epetra Y
    mapStrategy_->copyThyraIntoTpetra(tY, Y);
  } else {
    TEUCHOS_ASSERT(false);
  }
}

void TpetraOperatorWrapper::applyInverse(const Tpetra::MultiVector<ST, LO, GO, NT>& /* X */,
                                         Tpetra::MultiVector<ST, LO, GO, NT>& /* Y */,
                                         Teuchos::ETransp /* mode */, ST /* alpha */,
                                         ST /* beta */) const {
  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
                             "TpetraOperatorWrapper::applyInverse not implemented");
}

RCP<const Teuchos::Comm<Thyra::Ordinal> > TpetraOperatorWrapper::getThyraComm(
    const Thyra::LinearOpBase<ST>& inOp) const {
  RCP<const VectorSpaceBase<ST> > vs = inOp.domain();

  RCP<const SpmdVectorSpaceBase<ST> > spmd;
  RCP<const VectorSpaceBase<ST> > current = vs;
  while (current != Teuchos::null) {
    // try to cast to a product vector space first
    RCP<const ProductVectorSpaceBase<ST> > prod =
        rcp_dynamic_cast<const ProductVectorSpaceBase<ST> >(current);

    // figure out what type it is
    if (prod == Teuchos::null) {
      // hopfully this is a SPMD vector space
      spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<ST> >(current);

      break;
    } else  // get first convenient vector space
      current = prod->getBlock(0);
  }

  TEUCHOS_TEST_FOR_EXCEPTION(spmd == Teuchos::null, std::runtime_error,
                             "TpetraOperatorWrapper requires std::vector space "
                             "blocks to be SPMD std::vector spaces");

  return spmd->getComm();
  /*
    const Thyra::ConstLinearOperator<double> thyraOp = rcpFromRef(inOp);

    RCP<Epetra_Comm> rtn;
    // VectorSpace<double> vs = thyraOp.domain().getBlock(0);
    RCP<const VectorSpaceBase<double> > vs = thyraOp.domain().getBlock(0).constPtr();

    // search for an SpmdVectorSpaceBase object
    RCP<const SpmdVectorSpaceBase<double> > spmd;
    RCP<const VectorSpaceBase<double> > current = vs;
    while(current!=Teuchos::null) {
       // try to cast to a product vector space first
       RCP<const ProductVectorSpaceBase<double> > prod
             = rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(current);

       // figure out what type it is
       if(prod==Teuchos::null) {
          // hopfully this is a SPMD vector space
          spmd = rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(current);

          break;
       }
       else {
          // get first convenient vector space
          current = prod->getBlock(0);
       }
    }

    TEUCHOS_TEST_FOR_EXCEPTION(spmd==Teuchos::null, std::runtime_error,
                       "TpetraOperatorWrapper requires std::vector space "
                       "blocks to be SPMD std::vector spaces");

    const SerialComm<Thyra::Ordinal>* serialComm
      = dynamic_cast<const SerialComm<Thyra::Ordinal>*>(spmd->getComm().get());

  #ifdef HAVE_MPI
    const MpiComm<Thyra::Ordinal>* mpiComm
      = dynamic_cast<const MpiComm<Thyra::Ordinal>*>(spmd->getComm().get());

    TEUCHOS_TEST_FOR_EXCEPTION(mpiComm==0 && serialComm==0, std::runtime_error,
                       "SPMD std::vector space has a communicator that is "
                       "neither a serial comm nor an MPI comm");

    if (mpiComm != 0)
      {
        rtn = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
      }
    else
      {
        rtn = rcp(new Epetra_SerialComm());
      }
  #else
    TEUCHOS_TEST_FOR_EXCEPTION(serialComm==0, std::runtime_error,
                       "SPMD std::vector space has a communicator that is "
                       "neither a serial comm nor an MPI comm");
    rtn = rcp(new Epetra_SerialComm());

  #endif

    TEUCHOS_TEST_FOR_EXCEPTION(rtn.get()==0, std::runtime_error, "null communicator created");
    return rtn;
  */
}

int TpetraOperatorWrapper::GetBlockRowCount() {
  const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp =
      Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());

  return blkOp->productRange()->numBlocks();
}

int TpetraOperatorWrapper::GetBlockColCount() {
  const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp =
      Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());

  return blkOp->productDomain()->numBlocks();
}

Teuchos::RCP<const Tpetra::Operator<ST, LO, GO, NT> > TpetraOperatorWrapper::GetBlock(int i,
                                                                                      int j) const {
  const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp =
      Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());

  RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
      rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(blkOp->getBlock(i, j));

  return tOp->getConstTpetraOperator();
}

}  // namespace TpetraHelpers
}  // namespace Teko