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