Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_TpetraLinearOp_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: 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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_TPETRA_LINEAR_OP_HPP
00043 #define THYRA_TPETRA_LINEAR_OP_HPP
00044 
00045 #include "Thyra_TpetraLinearOp_decl.hpp"
00046 #include "Thyra_TpetraVectorSpace.hpp"
00047 
00048 #ifdef HAVE_THYRA_TPETRA_EPETRA
00049 #  include "Thyra_EpetraThyraWrappers.hpp"
00050 #endif
00051 
00052 namespace Thyra {
00053 
00054 
00055 #ifdef HAVE_THYRA_TPETRA_EPETRA
00056 
00057 // Utilites
00058 
00059 
00061   template<class Scalar, class LocalOrdinal, class GlobalOrdinal>
00062 class GetTpetraEpetraRowMatrixWrapper {
00063 public:
00064   template<class TpetraMatrixType>
00065   static
00066   RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
00067   get(const RCP<TpetraMatrixType> &tpetraMatrix)
00068     {
00069       return Teuchos::null;
00070     }
00071 };
00072 
00073 
00074 // NOTE: We could support other ordinal types, but we have to
00075 // specialize the EpetraRowMatrix
00076 template<>
00077 class GetTpetraEpetraRowMatrixWrapper<double, int, int> {
00078 public:
00079   template<class TpetraMatrixType>
00080   static
00081   RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
00082   get(const RCP<TpetraMatrixType> &tpetraMatrix)
00083     {
00084       return Teuchos::rcp(
00085         new Tpetra::EpetraRowMatrix<TpetraMatrixType>(tpetraMatrix,
00086           *get_Epetra_Comm(
00087             *convertTpetraToThyraComm(tpetraMatrix->getRowMap()->getComm())
00088             )
00089           )
00090         );
00091     }
00092 };
00093 
00094 
00095 #endif // HAVE_THYRA_TPETRA_EPETRA
00096 
00097 
00098 // Constructors/initializers
00099 
00100 
00101 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00102 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TpetraLinearOp()
00103 {}
00104 
00105 
00106 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00107 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::initialize(
00108   const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
00109   const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
00110   const RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
00111   )
00112 {
00113   initializeImpl(rangeSpace, domainSpace, tpetraOperator);
00114 }
00115 
00116 
00117 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00118 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::constInitialize(
00119   const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
00120   const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
00121   const RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
00122   )
00123 {
00124   initializeImpl(rangeSpace, domainSpace, tpetraOperator);
00125 }
00126 
00127 
00128 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00129 RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00130 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getTpetraOperator()
00131 {
00132   return tpetraOperator_.getNonconstObj();
00133 }
00134 
00135 
00136 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00137 RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00138 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraOperator() const
00139 {
00140   return tpetraOperator_;
00141 }
00142 
00143 
00144 // Public Overridden functions from LinearOpBase
00145 
00146 
00147 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00148 RCP<const Thyra::VectorSpaceBase<Scalar> >
00149 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::range() const
00150 {
00151   return rangeSpace_;
00152 }
00153 
00154 
00155 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00156 RCP<const Thyra::VectorSpaceBase<Scalar> >
00157 TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::domain() const
00158 {
00159   return domainSpace_;
00160 }
00161 
00162 
00163 // Overridden from EpetraLinearOpBase
00164 
00165 
00166 #ifdef HAVE_THYRA_TPETRA_EPETRA
00167 
00168 
00169 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00170 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstEpetraOpView(
00171   const Ptr<RCP<Epetra_Operator> > &epetraOp,
00172   const Ptr<EOpTransp> &epetraOpTransp,
00173   const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
00174   const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
00175   )
00176 {
00177   TEUCHOS_TEST_FOR_EXCEPT(true);
00178 }
00179 
00180 
00181 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00182 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getEpetraOpView(
00183   const Ptr<RCP<const Epetra_Operator> > &epetraOp,
00184   const Ptr<EOpTransp> &epetraOpTransp,
00185   const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
00186   const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
00187   ) const
00188 {
00189   using Teuchos::rcp_dynamic_cast;
00190   typedef Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpetraRowMatrix_t;
00191   if (nonnull(tpetraOperator_)) {
00192     if (is_null(epetraOp_)) {
00193       epetraOp_ = GetTpetraEpetraRowMatrixWrapper<Scalar,LocalOrdinal,GlobalOrdinal>::get(
00194         rcp_dynamic_cast<const TpetraRowMatrix_t>(tpetraOperator_.getConstObj(), true));
00195     }
00196     *epetraOp = epetraOp_;
00197     *epetraOpTransp = NOTRANS;
00198     *epetraOpApplyAs = EPETRA_OP_APPLY_APPLY;
00199     *epetraOpAdjointSupport = ( tpetraOperator_->hasTransposeApply()
00200       ? EPETRA_OP_ADJOINT_SUPPORTED : EPETRA_OP_ADJOINT_UNSUPPORTED );
00201   }
00202   else {
00203     *epetraOp = Teuchos::null;
00204   }
00205 }
00206 
00207 
00208 #endif // HAVE_THYRA_TPETRA_EPETRA
00209 
00210 
00211 // Protected Overridden functions from LinearOpBase
00212 
00213 
00214 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00215 bool TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::opSupportedImpl(
00216   Thyra::EOpTransp M_trans) const
00217 {
00218   if (is_null(tpetraOperator_))
00219     return false;
00220   if (M_trans == NOTRANS)
00221     return true;
00222   return tpetraOperator_->hasTransposeApply();
00223 }
00224 
00225 
00226 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00227 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::applyImpl(
00228   const Thyra::EOpTransp M_trans,
00229   const Thyra::MultiVectorBase<Scalar> &X_in,
00230   const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00231   const Scalar alpha,
00232   const Scalar beta
00233   ) const
00234 {
00235   using Teuchos::rcpFromRef;
00236   using Teuchos::rcpFromPtr;
00237   typedef TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>
00238     ConverterT;
00239   typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
00240     TpetraMultiVector_t;
00241 
00242   // Get Tpetra::MultiVector objects for X and Y
00243 
00244   const RCP<const TpetraMultiVector_t> tX =
00245     ConverterT::getConstTpetraMultiVector(rcpFromRef(X_in));
00246 
00247   const RCP<TpetraMultiVector_t> tY =
00248     ConverterT::getTpetraMultiVector(rcpFromPtr(Y_inout));
00249 
00250   const Teuchos::ETransp tTransp =
00251     ( M_trans == NOTRANS ? Teuchos::NO_TRANS : Teuchos::CONJ_TRANS );
00252 
00253   // Apply the operator
00254 
00255   tpetraOperator_->apply(*tX, *tY, tTransp, alpha, beta);
00256 
00257 }
00258 
00259 
00260 // private
00261 
00262 
00263 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00264 template<class TpetraOperator_t>
00265 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::initializeImpl(
00266   const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
00267   const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
00268   const RCP<TpetraOperator_t> &tpetraOperator
00269   )
00270 {
00271 #ifdef THYRA_DEBUG
00272   TEUCHOS_ASSERT(nonnull(rangeSpace));
00273   TEUCHOS_ASSERT(nonnull(domainSpace));
00274   TEUCHOS_ASSERT(nonnull(tpetraOperator));
00275   // ToDo: Assert that spaces are comparible with tpetraOperator
00276 #endif  
00277   rangeSpace_ = rangeSpace;
00278   domainSpace_ = domainSpace;
00279   tpetraOperator_ = tpetraOperator;
00280 }
00281 
00282 
00283 } // namespace Thyra
00284 
00285 
00286 #endif  // THYRA_TPETRA_LINEAR_OP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines