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 
00053     template <class S, class LO, class GO, class Node, class Kernel>
00054     class KernelOp : public Tpetra::Operator<S,LO,GO,Node> 
00055     {
00056       protected:
00057         RCP< const Import<LO,GO,Node> > _importer;
00058         RCP< const Export<LO,GO,Node> > _exporter;
00059         mutable RCP< MultiVector<S,LO,GO,Node> >  _importMV,  _exportMV;
00060         RCP< const Map<LO,GO,Node> > _rangeMap, _domainMap;
00061         mutable Kernel _kernel;
00062       public:
00063         KernelOp(Kernel kernel,
00064                  const RCP<const Map<LO,GO,Node> > & domainMap,
00065                  const RCP<const Map<LO,GO,Node> > & rangeMap,
00066                  const RCP<const Import<LO,GO,Node> > & importer,
00067                  const RCP<const Export<LO,GO,Node> > & exporter) 
00068         : _importer(importer), _exporter(exporter)
00069         , _rangeMap(rangeMap), _domainMap(domainMap)
00070         , _kernel(kernel)
00071         {
00072           std::string tfecfFuncName("KernelOp(kernel,domainMap,rangeMap,importer,exporter)");
00073           if (_rangeMap == null) _rangeMap = _domainMap;
00074           TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( _domainMap == null || _rangeMap == null, std::runtime_error,
00075               ":KernelOp(): neither domainMap nor rangeMap may be specified null:\ndomainMap: " << _domainMap << "\nrangeMap: " << _rangeMap << "\n");
00076 #ifdef HAVE_TPETRA_DEBUG
00077           TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( _rangeMap->getNode() != _domainMap->getNode(), std::runtime_error, ": all specified maps must have the same node.");
00078           if (_importer != null) {
00079             TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( !_importer->getSourceMap()->isSameAs(*_domainMap), std::runtime_error, ": domain map is not consistent with importer.");
00080             TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( _importer->getSourceMap()->getNode() != _domainMap->getNode(), std::runtime_error, ": all specified maps must have the same node.");
00081           }
00082           if (_exporter != null) {
00083             TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( !_exporter->getTargetMap()->isSameAs(*_rangeMap), std::runtime_error, ": range map is not consistent with importer.");
00084             TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( _exporter->getTargetMap()->getNode() != _domainMap->getNode(), std::runtime_error, ": all specified maps must have the same node.");
00085           }
00086 #endif
00087         }
00088         //
00089         const RCP<const Map<LO,GO,Node> > & getDomainMap() const { return _domainMap; }
00090         //
00091         const RCP<const Map<LO,GO,Node> > & getRangeMap()  const { return _rangeMap; }
00092         //
00093         void apply(const MultiVector<S,LO,GO,Node> &X, 
00094                          MultiVector<S,LO,GO,Node> &Y,
00095                    Teuchos::ETransp mode = Teuchos::NO_TRANS, 
00096                    S alpha = Teuchos::ScalarTraits<S>::one(), 
00097                    S beta = Teuchos::ScalarTraits<S>::zero()) const
00098         {
00099           const size_t numVectors = X.getNumVectors();
00100           RCP< MultiVector<S,LO,GO,Node> > mvec_inout;
00101           RCP< const MultiVector<S,LO,GO,Node> > mvec_in2;
00102           //
00103           if (_importer != null) {
00104             if (_importMV != null && _importMV->getNumVectors() != numVectors) _importMV = null;
00105             if (_importMV == null) _importMV = createMultiVector<S>(_importer->getTargetMap(), numVectors);
00106             _importMV->doImport( X, *_importer, INSERT );
00107             mvec_in2 = _importMV;
00108           }
00109           else {
00110             mvec_in2 = rcpFromRef(X);
00111           }
00112           //
00113           if (_exporter != null) {
00114             if (_exportMV != null && _exportMV->getNumVectors() != numVectors) _exportMV = null;
00115             if (_exportMV == null) _exportMV = createMultiVector<S>(_exporter->getSourceMap(), numVectors);
00116             mvec_inout = _exportMV;
00117           }
00118           else {
00119             mvec_inout = rcpFromRef(Y);
00120           }
00121           _kernel.setAlphaBeta(alpha,beta);
00122           //
00123           for (size_t j=0; j < numVectors; ++j) 
00124           {
00125             RCP<       Vector<S,LO,GO,Node> > vec_inout = mvec_inout->getVectorNonConst(j);
00126             RCP< const Vector<S,LO,GO,Node> > vec_in2   = mvec_in2->getVector(j);
00127             Tpetra::RTI::detail::binary_transform( *vec_inout, *vec_in2, _kernel );
00128           }
00129           // export
00130           if (_exporter != null) {
00131             Y.doExport(*_exportMV, *_exporter, ADD);
00132           }
00133         }
00134     };
00135 
00137     template <class S, class LO, class GO, class Node, class Kernel> 
00138     RCP< const KernelOp<S,LO,GO,Node,Kernel> >
00139     kernelOp(Kernel kernel, 
00140           const RCP<const Map<LO,GO,Node> > & domainMap,
00141           const RCP<const Map<LO,GO,Node> > & rangeMap = null,
00142           const RCP<const Import<LO,GO,Node> > & importer = null,
00143           const RCP<const Export<LO,GO,Node> > & exporter = null) 
00144     {
00145       return Teuchos::rcp(new KernelOp<S,LO,GO,Node,Kernel>(kernel,domainMap,rangeMap,importer,exporter) );
00146     }
00147 
00149     template <class S, class LO, class GO, class Node, class Op> 
00150     class BinaryOp : public KernelOp<S,LO,GO,Node,Tpetra::RTI::detail::BinaryFunctorAdapterWithAlphaBeta<Op,S> >
00151     {
00152       public:
00153         BinaryOp(Op op, 
00154               const RCP<const Map<LO,GO,Node> > & domainMap,
00155               const RCP<const Map<LO,GO,Node> > & rangeMap,
00156               const RCP<const Import<LO,GO,Node> > & importer,
00157               const RCP<const Export<LO,GO,Node> > & exporter) 
00158         : KernelOp<S,LO,GO,Node,Tpetra::RTI::detail::BinaryFunctorAdapterWithAlphaBeta<Op,S> >( Tpetra::RTI::detail::BinaryFunctorAdapterWithAlphaBeta<Op,S>(op), domainMap, rangeMap, importer, exporter ) {}
00159     };
00160 
00162     template <class S, class LO, class GO, class Node, class Op> 
00163     RCP< const BinaryOp<S,LO,GO,Node,Op> >
00164     binaryOp(Op op, 
00165              const RCP<const Map<LO,GO,Node> > & domainMap,
00166              const RCP<const Map<LO,GO,Node> > & rangeMap = null,
00167              const RCP<const Import<LO,GO,Node> > & importer = null,
00168              const RCP<const Export<LO,GO,Node> > & exporter = null) 
00169     {
00170       return Teuchos::rcp(new BinaryOp<S,LO,GO,Node,Op>(op,domainMap,rangeMap,importer,exporter) );
00171     }
00172 
00173   } // end of namespace RTI
00174 
00175 } // end of namespace Tpetra
00176 
00177 #endif // TPETRA_RTIOP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines