Tpetra Matrix/Vector Services Version of the Day
RTIInlineCG.hpp
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 <iostream>
00045 #include <functional>
00046 
00047 #include <Teuchos_CommandLineProcessor.hpp>
00048 #include <Teuchos_FancyOStream.hpp>
00049 #include <Teuchos_GlobalMPISession.hpp>
00050 #include <Teuchos_oblackholestream.hpp>
00051 #include <Teuchos_ParameterList.hpp>
00052 #include <Teuchos_TimeMonitor.hpp>
00053 #include <Teuchos_TypeNameTraits.hpp>
00054 #include <Teuchos_XMLParameterListHelpers.hpp>
00055 
00056 
00057 #include <Tpetra_CrsMatrix.hpp>
00058 #include <Tpetra_DefaultPlatform.hpp>
00059 
00060 // no gpu support for RTI
00061 #undef HAVE_KOKKOSCLASSIC_THRUST
00062 #include <Tpetra_HybridPlatform.hpp>
00063 
00064 #include <Tpetra_MatrixIO.hpp>
00065 #include <Tpetra_RTI.hpp>
00066 #include <Tpetra_RTIOp.hpp>
00067 #include <Tpetra_Vector.hpp>
00068 #include <Tpetra_Version.hpp>
00069 
00070 namespace Tpetra {
00071   namespace RTI {
00072     // specialization for pair
00073     template <class T1, class T2>
00074     class ZeroOp<std::pair<T1,T2>> {
00075       public:
00076       static inline std::pair<T1,T2> identity() {
00077         return std::make_pair( Teuchos::ScalarTraits<T1>::zero(), 
00078                                Teuchos::ScalarTraits<T2>::zero() );
00079       }
00080     };
00081   }
00082 }
00083 
00084 namespace TpetraExamples {
00085 
00086   using Teuchos::RCP;
00087   using Teuchos::ParameterList;
00088   using Teuchos::Time;
00089   using Teuchos::null;
00090   using std::binary_function;
00091   using std::pair;
00092   using std::make_pair;
00093   using std::plus;
00094   using std::multiplies;
00095 
00096   template <class T1, class T2, class Op>
00097   class pair_op : public binary_function<pair<T1,T2>,pair<T1,T2>,pair<T1,T2>> {
00098   private:
00099     Op op_;
00100   public:
00101     pair_op(Op op) : op_(op) {}
00102     inline pair<T1,T2> operator()(const pair<T1,T2>& a, const pair<T1,T2>& b) const {
00103       return make_pair(op_(a.first,b.first),op_(a.second,b.second));
00104     }
00105   };
00106 
00107   template <class T1, class T2, class Op>
00108   pair_op<T1,T2,Op> make_pair_op(Op op) { return pair_op<T1,T2,Op>(op); }
00109 
00111   template <class S>
00112   class DiagKernel : public Tpetra::RTI::detail::StdOpKernel<S>
00113   {
00114     protected:
00115       int N_;
00116       S kappa_;
00117     public:
00118       DiagKernel(int N, const S & kappa) : N_(N), kappa_(kappa) {}
00119       inline void execute(int i) { 
00120         this->_vec_inout[i] = this->_alpha * ((kappa_-1.0) * (S)(i)/(S)(N_-1) + 1.0) * this->_vec_in2[i] + this->_beta * this->_vec_inout[i];
00121       }
00122   };
00123 
00124 
00125 
00127   template <class S, class LO, class GO, class Node>      
00128   void RTICG(const RCP<const Tpetra::Operator<S,LO,GO,Node>> &A, RCP<Tpetra::Vector<S,LO,GO,Node>> b, 
00129              const RCP<Teuchos::FancyOStream> &out, ParameterList &plist)
00130   {
00131     using Teuchos::as;
00132     using Tpetra::RTI::ZeroOp;
00133     typedef Tpetra::Vector<S ,LO,GO,Node> Vector;
00134     typedef Tpetra::Operator<S,LO,GO,Node>    Op;
00135     typedef Teuchos::ScalarTraits<S>          ST;
00136     // get objects from level database
00137     const int numIters = plist.get<int>("numIters");
00138     const S tolerance  = plist.get<double>("tolerance");
00139     const int verbose  = plist.get<int>("verbose");
00140     RCP<Time> timer = Teuchos::TimeMonitor::getNewTimer(
00141                         "CG<"+Teuchos::TypeNameTraits<S>::name()+">"
00142                       );
00143     auto r = Tpetra::createVector<S>(A->getRangeMap()),
00144          p = Tpetra::createVector<S>(A->getRangeMap()),
00145         Ap = Tpetra::createVector<S>(A->getRangeMap()); 
00146     auto x = b;
00147 
00148     Teuchos::OSTab tab(out);
00149 
00150     if (verbose) {
00151       *out << "Beginning CG<" << Teuchos::TypeNameTraits<S>::name() << ">" << std::endl;
00152     }
00153     int k;
00154     { // begin timer scop
00155       Teuchos::TimeMonitor lcltimer(*timer);
00156       // initial guess is x=0, so initial residual is r = b - A*x = b
00157       const S r2 = TPETRA_BINARY_PRETRANSFORM_REDUCE(
00158                       r, b,                                         // fused: 
00159                       b,                                            //      : r = x  
00160                       r*r, ZeroOp<S>, plus<S>() );                  //      : sum r'*r
00161       // forget b; refer to it as x now
00162       b = null;
00163       // b comes in, x goes out. now we're done with b, so zero the solution.
00164       TPETRA_UNARY_TRANSFORM( x,  ST::zero() );                     // x = 0
00165       S rr = r2;
00166       TPETRA_BINARY_TRANSFORM( p, r,    r );                        // p = r
00168       for (k=0; k<numIters; ++k) 
00169       {
00170         A->apply(*p,*Ap);                                           // Ap = A*p
00171         S pAp = TPETRA_REDUCE2( p, Ap,     
00172                                 p*Ap, ZeroOp<S>, plus<S>() );       // p'*Ap
00173         const S alpha = rr / pAp;
00174         TPETRA_BINARY_TRANSFORM( x,    p,  x + alpha*p  );          // x = x + alpha*p
00175         S rrold = rr;
00176         rr = TPETRA_BINARY_PRETRANSFORM_REDUCE(
00177                                r, Ap,                               // fused:
00178                                r - alpha*Ap,                        //      : r - alpha*Ap
00179                                r*r, ZeroOp<S>, plus<S>() );         //      : sum r'*r
00180         if (verbose > 1) *out << "|res|/|res_0|: " << ST::squareroot(rr/r2) 
00181                               << std::endl;
00182         if (rr/r2 < tolerance*tolerance) {
00183           if (verbose) {
00184             *out << "Convergence detected!" << std::endl;
00185           }
00186           break;
00187         }
00188         const S beta = rr / rrold;
00189         TPETRA_BINARY_TRANSFORM( p, r,   r + beta*p );               // p = z + beta*p
00190       }
00191     } // end timer scop
00192     if (verbose) {
00193       *out << "Leaving recursiveFPCG<" << Teuchos::TypeNameTraits<S>::name() 
00194            << "> after " << k << " iterations." << std::endl;
00195     }
00196   }
00197 
00198 } // end of namespace TpetraExamples
00199 
00200 template <class S>
00201 class CGDriver {
00202   public:
00203   // input
00204   Teuchos::RCP<Teuchos::FancyOStream>     out; 
00205   Teuchos::RCP<Teuchos::ParameterList> params;
00206   // output
00207   bool                             testPassed;
00208 
00209   template <class Node> 
00210   void run(Teuchos::ParameterList &myMachPL, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, const Teuchos::RCP<Node> &node) 
00211   {
00212     using std::pair;
00213     using std::make_pair;
00214     using std::plus;
00215     using std::endl;
00216     using Teuchos::RCP;
00217     using Teuchos::ParameterList;
00218     using Teuchos::null;
00219     using TpetraExamples::make_pair_op;
00220     using Tpetra::RTI::reductionGlob;
00221     using Tpetra::RTI::ZeroOp;
00222     using Tpetra::RTI::binary_pre_transform_reduce;
00223     using Tpetra::RTI::binary_transform;
00224 
00225     // Static types
00226     typedef int                                    LO;
00227     typedef int                                    GO;
00228     typedef Tpetra::Operator<S,LO,GO,Node>   Operator;
00229     typedef Tpetra::Vector<S,LO,GO,Node>       Vector;
00230     const int        N = params->get<int>("size");
00231     const double kappa = params->get<double>("kappa");
00232 
00233     *out << "Running test with Node==" << Teuchos::typeName(*node) << " on rank " << comm->getRank() << "/" << comm->getSize() << std::endl;
00234 
00235     // create the operator
00236     *out << "Building problem of size " << N << " with condition number " << kappa << std::endl;
00237     RCP<const Operator> A;
00238     {
00239       auto map = Tpetra::createUniformContigMapWithNode<LO,GO,Node>(N,comm,node);
00240       A = Tpetra::RTI::kernelOp<S>( TpetraExamples::DiagKernel<S>(N,kappa), map );
00241     }
00242 
00243     testPassed = true;
00244 
00245     // choose a solution, compute a right-hand-side
00246     auto x = Tpetra::createVector<S>(A->getDomainMap()),
00247          b = Tpetra::createVector<S>(A->getDomainMap()),
00248       xhat = Tpetra::createVector<S>(A->getDomainMap());
00249     x->randomize();
00250     A->apply(*x,*b);
00251     TPETRA_BINARY_TRANSFORM(xhat, b,    b); // xhat = b
00252 
00253     // call the solve
00254     TpetraExamples::RTICG(A,xhat,out,*params);
00255 
00256     // check that residual is as requested
00257     {
00258       auto bhat = Tpetra::createVector<S>(A->getDomainMap());
00259       A->apply(*xhat,*bhat);
00260       // compute bhat-b, while simultaneously computing |bhat-b|^2 and |b|^2
00261       typedef ZeroOp<pair<S,S>> ZERO;
00262       const pair<S,S> nrms = TPETRA_BINARY_PRETRANSFORM_REDUCE(                                               
00263                                bhat, b,                                                                       
00264                                b - bhat,                                                                        
00265                                make_pair(bhat*bhat, b*b), ZERO, (make_pair_op<S,S>(plus<S>()))
00266                              );                                                                               
00267       //auto nrms = binary_pre_transform_reduce(*bhat, *b, 
00268       //                                        reductionGlob<ZeroOp<pair<S,S>>>( 
00269       //                                          [](S bhati, S bi){ return bi-bhati;}, // bhati = bi-bhat
00270       //                                          [](S bhati, S bi){ return make_pair(bhati*bhati, bi*bi); },
00271       //                                          make_pair_op<S,S>(plus<S>())) );
00272       const S enrm = Teuchos::ScalarTraits<S>::squareroot(nrms.first),
00273               bnrm = Teuchos::ScalarTraits<S>::squareroot(nrms.second);
00274       // check that residual is as requested
00275       *out << "|b - A*x|/|b|: " << enrm / bnrm << endl;
00276       const double tolerance = params->get<double>("tolerance");
00277       // give a little slack
00278       if (enrm / bnrm > 5*tolerance) testPassed = false;
00279     }
00280          
00281     // print timings
00282     Teuchos::TimeMonitor::summarize( *out );
00283   }
00284 };
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines