Thyra Version of the Day
Thyra_TpetraLinearOp_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thar: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_TPETRA_LINEAR_OP_HPP
00030 #define THYRA_TPETRA_LINEAR_OP_HPP
00031 
00032 #include "Thyra_TpetraLinearOp_decl.hpp"
00033 #include "Thyra_TpetraVectorSpace.hpp"
00034 
00035 
00036 #ifdef HAVE_THYRA_TPETRA_EPETRA
00037 
00038 
00039 #include "Thyra_EpetraThyraWrappers.hpp"
00040 
00041 
00042 namespace Thyra {
00043 
00044 
00045 // Utilites
00046 
00047 
00049   template<class Scalar, class LocalOrdinal, class GlobalOrdinal>
00050 class GetTpetraEpetraRowMatrixWrapper {
00051 public:
00052   template<class TpetraMatrixType>
00053   static
00054   RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
00055   get(const RCP<TpetraMatrixType> &tpetraMatrix)
00056     {
00057       return Teuchos::null;
00058     }
00059 };
00060 
00061 
00062 // NOTE: We could support other ordinal types, but we have to
00063 // specialize the EpetraRowMatrix
00064 template<>
00065 class GetTpetraEpetraRowMatrixWrapper<double, int, int> {
00066 public:
00067   template<class TpetraMatrixType>
00068   static
00069   RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
00070   get(const RCP<TpetraMatrixType> &tpetraMatrix)
00071     {
00072       return Teuchos::rcp(
00073         new Tpetra::EpetraRowMatrix<TpetraMatrixType>(tpetraMatrix,
00074           *get_Epetra_Comm(
00075             *convertTpetraToThyraComm(tpetraMatrix->getRowMap()->getComm())
00076             )
00077           )
00078         );
00079     }
00080 };
00081 
00082 
00083 #endif // HAVE_THYRA_TPETRA_EPETRA
00084 
00085 
00086 // Constructors/initializers
00087 
00088 
00089 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00090 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TpetraLinearOp()
00091 {}
00092 
00093 
00094 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00095 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::initialize(
00096   const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
00097   const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
00098   const RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
00099   )
00100 {
00101   initializeImpl(rangeSpace, domainSpace, tpetraOperator);
00102 }
00103 
00104 
00105 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00106 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::constInitialize(
00107   const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
00108   const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
00109   const RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
00110   )
00111 {
00112   initializeImpl(rangeSpace, domainSpace, tpetraOperator);
00113 }
00114 
00115 
00116 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00117 RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00118 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getTpetraOperator()
00119 {
00120   return tpetraOperator_.getNonconstObj();
00121 }
00122 
00123 
00124 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00125 RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00126 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraOperator() const
00127 {
00128   return tpetraOperator_;
00129 }
00130 
00131 
00132 // Public Overridden functions from LinearOpBase
00133 
00134 
00135 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00136 RCP<const Thyra::VectorSpaceBase<Scalar> >
00137 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::range() const
00138 {
00139   return rangeSpace_;
00140 }
00141 
00142 
00143 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00144 RCP<const Thyra::VectorSpaceBase<Scalar> >
00145 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::domain() const
00146 {
00147   return domainSpace_;
00148 }
00149 
00150 
00151 // Overridden from EpetraLinearOpBase
00152 
00153 
00154 #ifdef HAVE_THYRA_TPETRA_EPETRA
00155 
00156 
00157 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00158 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstEpetraOpView(
00159   const Ptr<RCP<Epetra_Operator> > &epetraOp,
00160   const Ptr<EOpTransp> &epetraOpTransp,
00161   const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
00162   const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
00163   )
00164 {
00165   TEST_FOR_EXCEPT(true);
00166 }
00167 
00168 
00169 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00170 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getEpetraOpView(
00171   const Ptr<RCP<const Epetra_Operator> > &epetraOp,
00172   const Ptr<EOpTransp> &epetraOpTransp,
00173   const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
00174   const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
00175   ) const
00176 {
00177   using Teuchos::rcp_dynamic_cast;
00178   typedef Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpetraRowMatrix_t;
00179   if (nonnull(tpetraOperator_)) {
00180     if (is_null(epetraOp_)) {
00181       epetraOp_ = GetTpetraEpetraRowMatrixWrapper<Scalar,LocalOrdinal,GlobalOrdinal>::get(
00182         rcp_dynamic_cast<const TpetraRowMatrix_t>(tpetraOperator_.getConstObj(), true));
00183     }
00184     *epetraOp = epetraOp_;
00185     *epetraOpTransp = NOTRANS;
00186     *epetraOpApplyAs = EPETRA_OP_APPLY_APPLY;
00187     *epetraOpAdjointSupport = ( tpetraOperator_->hasTransposeApply()
00188       ? EPETRA_OP_ADJOINT_SUPPORTED : EPETRA_OP_ADJOINT_UNSUPPORTED );
00189   }
00190   else {
00191     *epetraOp = Teuchos::null;
00192   }
00193 }
00194 
00195 
00196 #endif // HAVE_THYRA_TPETRA_EPETRA
00197 
00198 
00199 // Protected Overridden functions from LinearOpBase
00200 
00201 
00202 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00203 bool TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::opSupportedImpl(
00204   Thyra::EOpTransp M_trans) const
00205 {
00206   if (is_null(tpetraOperator_))
00207     return false;
00208   if (M_trans == NOTRANS)
00209     return true;
00210   return tpetraOperator_->hasTransposeApply();
00211 }
00212 
00213 
00214 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00215 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::applyImpl(
00216   const Thyra::EOpTransp M_trans,
00217   const Thyra::MultiVectorBase<Scalar> &X_in,
00218   const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00219   const Scalar alpha,
00220   const Scalar beta
00221   ) const
00222 {
00223   using Teuchos::rcpFromRef;
00224   using Teuchos::rcpFromPtr;
00225   typedef TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>
00226     ConverterT;
00227   typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
00228     TpetraMultiVector_t;
00229 
00230   // Get Tpetra::MultiVector objects for X and Y
00231 
00232   const RCP<const TpetraMultiVector_t> tX =
00233     ConverterT::getConstTpetraMultiVector(rcpFromRef(X_in));
00234 
00235   const RCP<TpetraMultiVector_t> tY =
00236     ConverterT::getTpetraMultiVector(rcpFromPtr(Y_inout));
00237 
00238   const Teuchos::ETransp tTransp =
00239     ( M_trans == NOTRANS ? Teuchos::NO_TRANS : Teuchos::CONJ_TRANS );
00240 
00241   // Apply the operator
00242 
00243   tpetraOperator_->apply(*tX, *tY, tTransp, alpha, beta);
00244 
00245 }
00246 
00247 
00248 // private
00249 
00250 
00251 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00252 template<class TpetraOperator_t>
00253 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::initializeImpl(
00254   const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
00255   const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
00256   const RCP<TpetraOperator_t> &tpetraOperator
00257   )
00258 {
00259 #ifdef THYRA_DEBUG
00260   TEUCHOS_ASSERT(nonnull(rangeSpace));
00261   TEUCHOS_ASSERT(nonnull(domainSpace));
00262   TEUCHOS_ASSERT(nonnull(tpetraOperator));
00263   // ToDo: Assert that spaces are comparible with tpetraOperator
00264 #endif  
00265   rangeSpace_ = rangeSpace;
00266   domainSpace_ = domainSpace;
00267   tpetraOperator_ = tpetraOperator;
00268 }
00269 
00270 
00271 } // namespace Thyra
00272 
00273 
00274 #endif  // THYRA_TPETRA_LINEAR_OP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines