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