Tpetra Matrix/Vector Services Version of the Day
Tpetra_RTIOp.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_RTIOP_HPP
00043 #define TPETRA_RTIOP_HPP
00044 
00045 #include "Tpetra_RTI.hpp"
00046 #include "Tpetra_RTI_detail.hpp"
00047 
00048 namespace Tpetra {
00049 
00050   namespace RTI {
00051 
00060     template <class S, class LO, class GO, class Node, class Kernel>
00061     class KernelOp : public Tpetra::Operator<S,LO,GO,Node> {
00062     protected:
00063       RCP<const Import<LO,GO,Node> > _importer;
00064       RCP<const Export<LO,GO,Node> > _exporter;
00065       mutable RCP<MultiVector<S,LO,GO,Node> > _importMV, _exportMV;
00066       RCP<const Map<LO,GO,Node> > _rangeMap, _domainMap;
00067       mutable Kernel _kernel;
00068       
00069     public:
00071       KernelOp (Kernel kernel,
00072     const RCP<const Map<LO,GO,Node> > & domainMap,
00073     const RCP<const Map<LO,GO,Node> > & rangeMap,
00074     const RCP<const Import<LO,GO,Node> > & importer,
00075     const RCP<const Export<LO,GO,Node> > & exporter) 
00076         : _importer (importer), _exporter (exporter)
00077         , _rangeMap (rangeMap), _domainMap (domainMap)
00078         , _kernel (kernel)
00079       {
00080   const char tfecfFuncName[] = "KernelOp(kernel,domainMap,rangeMap,importer,exporter)";
00081   if (_rangeMap == null) {
00082     _rangeMap = _domainMap;
00083   }
00084   TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00085           _domainMap == null || _rangeMap == null, std::runtime_error,
00086     ": neither domainMap nor rangeMap may be specified null:\ndomainMap: " 
00087     << _domainMap << "\nrangeMap: " << _rangeMap << "\n");
00088 #ifdef HAVE_TPETRA_DEBUG
00089         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00090           _rangeMap->getNode() != _domainMap->getNode(), std::runtime_error, 
00091     ": all specified maps must have the same Node instance.");
00092   if (_importer != null) {
00093     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 
00094             ! _importer->getSourceMap ()->isSameAs (*_domainMap), 
00095       std::runtime_error, 
00096       ": domain Map is not consistent with importer.");
00097     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00098             _importer->getSourceMap ()->getNode () != _domainMap->getNode (),
00099       std::runtime_error, 
00100       ": all specified Maps must have the same Node instance.");
00101   }
00102   if (_exporter != null) {
00103     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00104             ! _exporter->getTargetMap ()->isSameAs (*_rangeMap), 
00105       std::runtime_error, ": range Map is not consistent with importer.");
00106           TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00107             _exporter->getTargetMap ()->getNode () != _domainMap->getNode (), 
00108       std::runtime_error, 
00109       ": all specified Maps must have the same Node instance.");
00110   }
00111 #endif // HAVE_TPETRA_DEBUG
00112       }
00113 
00114       RCP<const Map<LO,GO,Node> > getDomainMap () const { return _domainMap; }
00115       RCP<const Map<LO,GO,Node> > getRangeMap ()  const { return _rangeMap; }
00116 
00117       void
00118       apply (const MultiVector<S,LO,GO,Node> &X, 
00119        MultiVector<S,LO,GO,Node> &Y,
00120        Teuchos::ETransp mode = Teuchos::NO_TRANS, 
00121        S alpha = Teuchos::ScalarTraits<S>::one (), 
00122        S beta = Teuchos::ScalarTraits<S>::zero ()) const
00123       {
00124   const size_t numVectors = X.getNumVectors ();
00125   RCP<MultiVector<S,LO,GO,Node> > mvec_inout;
00126   RCP<const MultiVector<S,LO,GO,Node> > mvec_in2;
00127 
00128   if (_importer != null) {
00129     if (_importMV != null && _importMV->getNumVectors () != numVectors) {
00130       _importMV = null;
00131     }
00132     if (_importMV == null) {
00133       _importMV = createMultiVector<S> (_importer->getTargetMap (), numVectors);
00134     }
00135     _importMV->doImport (X, *_importer, INSERT);
00136     mvec_in2 = _importMV;
00137   }
00138   else {
00139     mvec_in2 = rcpFromRef(X);
00140   }
00141 
00142   if (_exporter != null) {
00143     if (_exportMV != null && _exportMV->getNumVectors () != numVectors) {
00144       _exportMV = null;
00145     }
00146     if (_exportMV == null) {
00147       _exportMV = createMultiVector<S> (_exporter->getSourceMap (), numVectors);
00148     }
00149     mvec_inout = _exportMV;
00150   }
00151   else {
00152     mvec_inout = rcpFromRef (Y);
00153   }
00154   _kernel.setAlphaBeta (alpha, beta);
00155   //
00156   for (size_t j=0; j < numVectors; ++j) {
00157     RCP<       Vector<S,LO,GO,Node> > vec_inout = mvec_inout->getVectorNonConst(j);
00158     RCP< const Vector<S,LO,GO,Node> > vec_in2   = mvec_in2->getVector(j);
00159     Tpetra::RTI::detail::binary_transform( *vec_inout, *vec_in2, _kernel );
00160   }
00161   // export
00162   if (_exporter != null) {
00163     Y.doExport (*_exportMV, *_exporter, ADD);
00164   }
00165       }
00166     };
00167 
00169     template <class S, class LO, class GO, class Node, class Kernel> 
00170     RCP<const KernelOp<S,LO,GO,Node,Kernel> >
00171     kernelOp (Kernel kernel, 
00172         const RCP<const Map<LO,GO,Node> > & domainMap,
00173         const RCP<const Map<LO,GO,Node> > & rangeMap = null,
00174         const RCP<const Import<LO,GO,Node> > & importer = null,
00175         const RCP<const Export<LO,GO,Node> > & exporter = null) 
00176     {
00177       return Teuchos::rcp (new KernelOp<S,LO,GO,Node,Kernel> (kernel, domainMap, rangeMap, 
00178                     importer, exporter) );
00179     }
00180 
00182     template <class S, class LO, class GO, class Node, class Op> 
00183     class BinaryOp : 
00184       public KernelOp<S,LO,GO,Node,Tpetra::RTI::detail::BinaryFunctorAdapterWithAlphaBeta<Op,S> >
00185     {
00186     public:
00187       BinaryOp (Op op, 
00188     const RCP<const Map<LO,GO,Node> > & domainMap,
00189     const RCP<const Map<LO,GO,Node> > & rangeMap,
00190     const RCP<const Import<LO,GO,Node> > & importer,
00191     const RCP<const Export<LO,GO,Node> > & exporter) 
00192         : KernelOp<S,LO,GO,Node,Tpetra::RTI::detail::BinaryFunctorAdapterWithAlphaBeta<Op,S> >( Tpetra::RTI::detail::BinaryFunctorAdapterWithAlphaBeta<Op,S>(op), domainMap, rangeMap, importer, exporter ) {}
00193     };
00194 
00196     template <class S, class LO, class GO, class Node, class Op> 
00197     RCP<const BinaryOp<S,LO,GO,Node,Op> >
00198     binaryOp (Op op, 
00199         const RCP<const Map<LO,GO,Node> > & domainMap,
00200         const RCP<const Map<LO,GO,Node> > & rangeMap = null,
00201         const RCP<const Import<LO,GO,Node> > & importer = null,
00202         const RCP<const Export<LO,GO,Node> > & exporter = null) 
00203     {
00204       return Teuchos::rcp (new BinaryOp<S,LO,GO,Node,Op> (op, domainMap, rangeMap, 
00205                 importer, exporter) );
00206     }
00207   } // namespace RTI
00208 } // namespace Tpetra
00209 
00210 #endif // TPETRA_RTIOP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines