Tpetra Matrix/Vector Services Version of the Day
RTIOperatorExample.cpp
Go to the documentation of this file.
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 #include <Teuchos_CommandLineProcessor.hpp>
00047 #include <Teuchos_Tuple.hpp>
00048 #include <Teuchos_TimeMonitor.hpp>
00049 
00050 #include <Kokkos_DefaultNode.hpp>
00051 
00052 #include "Tpetra_Version.hpp"
00053 #include "Tpetra_DefaultPlatform.hpp"
00054 #include "Tpetra_MultiVector.hpp"
00055 #include "Tpetra_CrsMatrix.hpp"
00056 #include "Tpetra_RTIOp.hpp"
00057 
00058 #include <iostream>
00059 #include <iomanip>
00060 #include <functional>
00061 
00067 namespace TpetraExamples {
00068 
00070   template <class S>
00071   class FDStencil : public Tpetra::RTI::detail::StdOpKernel<S> 
00072   {
00073     protected:
00074       int _myImageID, _numImages;
00075       int _n;
00076       S _sub, _diag, _super;
00077     public:
00078       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) {}
00079       inline void execute(int i) const
00080       { 
00081         S res = _diag * this->_vec_in2[i];
00082         if (i >  0 || _myImageID != 0           ) res +=   _sub * this->_vec_in2[i-1];
00083         if (i < _n || _myImageID != _numImages-1) res += _super * this->_vec_in2[i+1];
00084         if (this->_alpha == Teuchos::ScalarTraits<S>::one() && this->_beta == Teuchos::ScalarTraits<S>::zero()) {
00085           this->_vec_inout[i] = res;
00086         }
00087         else {
00088           this->_vec_inout[i] = this->_alpha * res + this->_beta * this->_vec_inout[i];
00089         }
00090       }
00091   };
00092 
00094   template <class S>
00095   class ScaleKernel : public Tpetra::RTI::detail::StdOpKernel<S>
00096   {
00097     protected:
00098       S _gamma;
00099     public:
00100       ScaleKernel(const S & gamma) : _gamma(gamma) {}
00101       inline void execute(int i) { 
00102         if (this->_alpha == Teuchos::ScalarTraits<S>::one() && this->_beta == Teuchos::ScalarTraits<S>::zero()) {
00103           this->_vec_inout[i] = _gamma * this->_vec_in2[i];
00104         }
00105         else {
00106           this->_vec_inout[i] = this->_alpha * _gamma * this->_vec_in2[i] + this->_beta * this->_vec_inout[i];
00107         }
00108       }
00109   };
00110 
00111 }
00112 
00113 int main(int argc, char *argv[])
00114 {
00115   Teuchos::oblackholestream blackhole;
00116   Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);
00117 
00118   using Teuchos::TimeMonitor;
00119   using TpetraExamples::FDStencil;
00120   using TpetraExamples::ScaleKernel;
00121 
00122   // 
00123   // Get the default communicator and node
00124   //
00125   auto &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
00126   auto comm = platform.getComm();
00127   auto node = platform.getNode();
00128   const int myImageID = comm->getRank();
00129   const int numImages = comm->getSize();
00130 
00131   //
00132   // Get example parameters from command-line processor
00133   //  
00134   bool verbose = (myImageID==0);
00135   int numGlobal_user = 100*comm->getSize();
00136   int numTimeTrials = 3;
00137   Teuchos::CommandLineProcessor cmdp(false,true);
00138   cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
00139   cmdp.setOption("global-size",&numGlobal_user,"Global test size.");
00140   cmdp.setOption("num-time-trials",&numTimeTrials,"Number of trials in timing loops.");
00141   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00142     return -1;
00143   }
00144 
00145   // 
00146   // Say hello, print some communicator info
00147   //
00148   if (verbose) {
00149     std::cout << "\n" << Tpetra::version() << std::endl;
00150     std::cout << "Comm info: " << *comm;
00151     std::cout << "Node type: " << Teuchos::typeName(*node) << std::endl;
00152     std::cout << std::endl;
00153   }
00154 
00155   // 
00156   // Create a simple map for domain and range
00157   // 
00158   Tpetra::global_size_t numGlobalRows = numGlobal_user;
00159   auto map = Tpetra::createUniformContigMapWithNode<int,int>(numGlobalRows, comm, node);
00160   // const size_t numLocalRows = map->getNodeNumElements();
00161   auto x = Tpetra::createVector<double>(map),
00162        y = Tpetra::createVector<double>(map);
00163 
00164   // Create a simple diagonal operator using lambda function
00165   auto fTwoOp = Tpetra::RTI::binaryOp<double>( [](double /*y*/, double x) { return 2.0 * x; } , map );
00166   // y = 3*fTwoOp*x + 2*y = 3*2*1 + 2*1 = 8
00167   x->putScalar(1.0);
00168   y->putScalar(1.0);
00169   fTwoOp->apply( *x, *y, Teuchos::NO_TRANS, 3.0, 2.0 );
00170   // check that y == eights
00171   double norm = y->norm1();
00172   if (verbose) {
00173     std::cout << "Tpetra::RTI::binaryOp" << std::endl
00174               << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
00175               << ", expected value is " << numGlobalRows * 8.0 << std::endl;
00176   }
00177 
00178   // Create the same diagonal operator using a Kokkos kernel
00179   auto kTwoOp = Tpetra::RTI::kernelOp<double>( ScaleKernel<double>(2.0), map );
00180   // y = 0.5*kTwop*x + 0.75*y = 0.5*2*1 + 0.75*8 = 7
00181   kTwoOp->apply( *x, *y, Teuchos::NO_TRANS, 0.5, 0.75 );
00182   // check that y == sevens
00183   norm = y->norm1();
00184   if (verbose) {
00185     std::cout << "Tpetra::RTI::kernelOp" << std::endl
00186               << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
00187               << ", expected value is " << numGlobalRows * 7.0 << std::endl;
00188   }
00189 
00190   //
00191   // Create a finite-difference stencil using a Kokkos kernel and non-trivial maps
00192   //
00193   decltype(map) colmap;
00194   if (numImages > 1) {
00195     Teuchos::Array<int>           colElements;
00196     Teuchos::ArrayView<const int> rowElements = map->getNodeElementList();
00197     // This isn't the most efficient Map/Import layout, but it makes for a very straight-forward kernel
00198     if (myImageID != 0) colElements.push_back( map->getMinGlobalIndex() - 1 );
00199     colElements.insert(colElements.end(), rowElements.begin(), rowElements.end());
00200     if (myImageID != numImages-1) colElements.push_back( map->getMaxGlobalIndex() + 1 );
00201     colmap = Tpetra::createNonContigMapWithNode<int,int>(colElements(), comm, node);
00202   }
00203   else {
00204     colmap = map;
00205   }
00206   auto importer = createImport(map,colmap);
00207   // Finite-difference kernel = tridiag(-1, 2, -1)
00208   FDStencil<double> kern(myImageID, numImages, map->getNodeNumElements(), -1.0, 2.0, -1.0);
00209   auto FDStencilOp = Tpetra::RTI::kernelOp<double>( kern, map, map, importer );
00210   // x = ones(), FD(x) = [1 zeros() 1]
00211   auto timeFDStencil = TimeMonitor::getNewTimer("FD RTI Stencil");
00212   {
00213     TimeMonitor lcltimer(*timeFDStencil);
00214     for (int l=0; l != numTimeTrials; ++l) {
00215       FDStencilOp->apply( *x, *y );
00216     }
00217   }
00218   norm = y->norm1();
00219   if (verbose) {
00220     std::cout << std::endl
00221               << "TpetraExamples::FDStencil" << std::endl
00222               << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
00223               << ", expected value is " << 2.0 << std::endl;
00224   }
00225 
00226   //
00227   // Create a finite-difference stencil using a CrsMatrix
00228   // 
00229   auto FDMatrix = Tpetra::createCrsMatrix<double>(map);  
00230   for (int r=map->getMinGlobalIndex(); r <= map->getMaxGlobalIndex(); ++r) {
00231     if (r == map->getMinAllGlobalIndex()) {
00232       FDMatrix->insertGlobalValues(r, Teuchos::tuple<int>(r,r+1), Teuchos::tuple<double>(2.0,-1.0));
00233     }
00234     else if (r == map->getMaxAllGlobalIndex()) {
00235       FDMatrix->insertGlobalValues(r, Teuchos::tuple<int>(r-1,r), Teuchos::tuple<double>(-1.0,2.0));
00236     }
00237     else {
00238       FDMatrix->insertGlobalValues(r, Teuchos::tuple<int>(r-1,r,r+1), Teuchos::tuple<double>(-1.0,2.0,-1.0));
00239     }
00240   }
00241   FDMatrix->fillComplete();
00242   auto timeFDMatrix = TimeMonitor::getNewTimer("FD CrsMatrix");
00243   {
00244     TimeMonitor lcltimer(*timeFDMatrix);
00245     for (int l=0; l != numTimeTrials; ++l) {
00246       FDMatrix->apply(*x, *y);
00247     }
00248   }
00249 
00250   //
00251   // Print timings
00252   //
00253   if (verbose) {
00254     std::cout << std::endl;
00255     TimeMonitor::summarize( std::cout );
00256   }
00257 
00258   if (verbose) {
00259     std::cout << "\nEnd Result: TEST PASSED" << std::endl;
00260   }
00261   return 0;
00262 }
00263 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines