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