Tpetra Matrix/Vector Services Version of the Day
RTITutorialOperator.cpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //          Tpetra: Templated Linear Algebra Services Package
00006 //                 Copyright (2008) Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 // @HEADER
00042 */
00043 
00044 #include <Teuchos_GlobalMPISession.hpp>
00045 #include <Teuchos_oblackholestream.hpp>
00046 
00047 #include <Kokkos_DefaultNode.hpp>
00048 #include <Tpetra_DefaultPlatform.hpp>
00049 #include <Tpetra_Version.hpp>
00050 #include <Tpetra_Vector.hpp>
00051 #include <Tpetra_RTIOp.hpp>
00052 
00053 #include <iostream>
00054 #include <iomanip>
00055 
00057 template <class S>
00058 class FDStencil : public Tpetra::RTI::detail::StdOpKernel<S> 
00059 {
00060   protected:
00061     int _myImageID, _numImages;
00062     int _n;
00063     S _sub, _diag, _super;
00064   public:
00065     FDStencil(int myImageID, int numImages, int n, const S &sub, const S &diag, const S &super) : _myImageID(myImageID), _numImages(numImages), _n(n), _sub(sub), _diag(diag), _super(super) {}
00066     inline void execute(int i) const
00067     { 
00068       S res = _diag * this->_vec_in2[i];
00069       if (i >  0 || _myImageID != 0           ) res +=   _sub * this->_vec_in2[i-1];
00070       if (i < _n || _myImageID != _numImages-1) res += _super * this->_vec_in2[i+1];
00071       if (this->_alpha == Teuchos::ScalarTraits<S>::one() && this->_beta == Teuchos::ScalarTraits<S>::zero()) {
00072         this->_vec_inout[i] = res;
00073       }
00074       else {
00075         this->_vec_inout[i] = this->_alpha * res + this->_beta * this->_vec_inout[i];
00076       }
00077     }
00078 };
00079 
00081 template <class S>
00082 class ScaleKernel : public Tpetra::RTI::detail::StdOpKernel<S>
00083 {
00084   protected:
00085     S _gamma;
00086   public:
00087     ScaleKernel(const S & gamma) : _gamma(gamma) {}
00088     inline void execute(int i) { 
00089       if (this->_alpha == Teuchos::ScalarTraits<S>::one() && this->_beta == Teuchos::ScalarTraits<S>::zero()) {
00090         this->_vec_inout[i] = _gamma * this->_vec_in2[i];
00091       }
00092       else {
00093         this->_vec_inout[i] = this->_alpha * _gamma * this->_vec_in2[i] + this->_beta * this->_vec_inout[i];
00094       }
00095     }
00096 };
00097 
00098 int main(int argc, char *argv[])
00099 {
00100   Teuchos::oblackholestream blackhole;
00101   Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);
00102 
00103   // 
00104   // Get the default communicator and node
00105   //
00106   auto &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
00107   auto comm = platform.getComm();
00108   auto node = platform.getNode();
00109   const int myImageID = comm->getRank();
00110   const int numImages = comm->getSize();
00111   const bool verbose = (myImageID==0);
00112 
00113   // 
00114   // Say hello, print some communicator info
00115   //
00116   if (verbose) {
00117     std::cout << "\n" << Tpetra::version() << std::endl;
00118     std::cout << "Comm info: " << *comm;
00119     std::cout << "Node type: " << Teuchos::typeName(*node) << std::endl;
00120     std::cout << std::endl;
00121   }
00122 
00123   // 
00124   // Create a simple map for domain and range
00125   // 
00126   Tpetra::global_size_t numGlobalRows = 1000*numImages;
00127   auto map = Tpetra::createUniformContigMapWithNode<int,int>(numGlobalRows, comm, node);
00128   auto x = Tpetra::createVector<double>(map),
00129        y = Tpetra::createVector<double>(map);
00130 
00131   // Create a simple diagonal operator using lambda function
00132   auto fTwoOp = Tpetra::RTI::binaryOp<double>( [](double /*y*/, double x) { return 2.0 * x; } , map );
00133   // y = 3*fTwoOp*x + 2*y = 3*2*1 + 2*1 = 8
00134   x->putScalar(1.0);
00135   y->putScalar(1.0);
00136   fTwoOp->apply( *x, *y, Teuchos::NO_TRANS, 3.0, 2.0 );
00137   // check that y == eights
00138   double norm = y->norm1();
00139   if (verbose) {
00140     std::cout << "Tpetra::RTI::binaryOp" << std::endl
00141               << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
00142               << ", expected value is " << numGlobalRows * 8.0 << std::endl;
00143   }
00144 
00145   //
00146   // Create a finite-difference stencil using a Kokkos kernel and non-trivial maps
00147   //
00148   decltype(map) colmap;
00149   if (numImages > 1) {
00150     Teuchos::Array<int>           colElements;
00151     Teuchos::ArrayView<const int> rowElements = map->getNodeElementList();
00152     // This isn't the most efficient Map/Import layout, but it makes for a very straight-forward kernel
00153     if (myImageID != 0) colElements.push_back( map->getMinGlobalIndex() - 1 );
00154     colElements.insert(colElements.end(), rowElements.begin(), rowElements.end());
00155     if (myImageID != numImages-1) colElements.push_back( map->getMaxGlobalIndex() + 1 );
00156     colmap = Tpetra::createNonContigMapWithNode<int,int>(colElements(), comm, node);
00157   }
00158   else {
00159     colmap = map;
00160   }
00161   auto importer = createImport(map,colmap);
00162   // Finite-difference kernel = tridiag(-1, 2, -1)
00163   FDStencil<double> kern(myImageID, numImages, map->getNodeNumElements(), -1.0, 2.0, -1.0);
00164   auto FDStencilOp = Tpetra::RTI::kernelOp<double>( kern, map, map, importer );
00165   // x = ones(), FD(x) = [1 zeros() 1]
00166   FDStencilOp->apply( *x, *y );
00167   norm = y->norm1();
00168   if (verbose) {
00169     std::cout << std::endl
00170               << "TpetraExamples::FDStencil" << std::endl
00171               << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
00172               << ", expected value is " << 2.0 << std::endl;
00173   }
00174 
00175   std::cout << "\nEnd Result: TEST PASSED" << std::endl;
00176   return 0;
00177 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines