Tpetra Matrix/Vector Services Version of the Day
IRTR_solvers.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 #ifndef MULTIPRECCG_HPP_
00045 #define MULTIPRECCG_HPP_
00046 
00047 #include <Teuchos_TimeMonitor.hpp>
00048 #include <Teuchos_TypeNameTraits.hpp>
00049 #include <Teuchos_ParameterList.hpp>
00050 #include <Teuchos_XMLParameterListHelpers.hpp>
00051 #include <Teuchos_FancyOStream.hpp>
00052 
00053 #include <Tpetra_CrsMatrix.hpp>
00054 #include <Tpetra_Vector.hpp>
00055 #include <Tpetra_RTI.hpp>
00056 #include <Tpetra_MatrixIO.hpp>
00057 
00058 #include <iostream>
00059 #include <functional>
00060 
00061 #ifdef HAVE_TPETRA_QD
00062 # include <qd/qd_real.h>
00063 #endif
00064 
00065 #define XFORM1     TPETRA_UNARY_TRANSFORM
00066 #define XFORM2     TPETRA_BINARY_TRANSFORM
00067 #define XFORM3     TPETRA_TERTIARY_TRANSFORM
00068 
00069 #define REDUCE1(in, gexp) \
00070   Tpetra::RTI::reduce( *in,                                                           \
00071     Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((in)->meanValue())>>(     \
00072                                [=]( decltype((in)->meanValue()) in ){ return gexp; }, \
00073                                std::plus<decltype((in)->meanValue())>() ) )
00074 
00075 #define REDUCE2(in1, in2, gexp) \
00076   Tpetra::RTI::reduce( *in1, *in2,                                                      \
00077     Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((in1)->meanValue())>>(      \
00078                                [=]( decltype((in1)->meanValue()) in1,                   \
00079                                     decltype((in2)->meanValue()) in2 ) { return gexp; },\
00080                                std::plus<decltype((in1)->meanValue())>() ) )
00081 
00082 #define XFORM2RED(out,in,texp,gexp) \
00083   Tpetra::RTI::binary_pre_transform_reduce( *out, *in,                                  \
00084     Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((out)->meanValue())>>(      \
00085                                [=]( decltype((out)->meanValue()) out,                   \
00086                                     decltype((in)->meanValue()) in )  { return texp; }, \
00087                                [=]( decltype((out)->meanValue()) out,                   \
00088                                     decltype((in)->meanValue()) in ) { return gexp; },  \
00089                                std::plus<decltype((out)->meanValue())>() ) )
00090 
00091 #define XFORM3RED(out,in2,in3, texp, gexp)                                          \
00092   Tpetra::RTI::tertiary_pre_transform_reduce( *out, *in2, *in3,                     \
00093     Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((out)->meanValue())>>(  \
00094                                [=]( decltype((out)->meanValue()) out,               \
00095                                     decltype((in2)->meanValue()) in2,               \
00096                                     decltype((in3)->meanValue()) in3 )              \
00097                                     { return texp; },                               \
00098                                [=]( decltype((out)->meanValue()) out,               \
00099                                     decltype((in2)->meanValue()) in2,               \
00100                                     decltype((in3)->meanValue()) in3 )              \
00101                                     { return gexp; },                               \
00102                                std::plus<decltype((out)->meanValue())>() ) )
00103 
00104 namespace TpetraExamples {
00105 
00106   using Teuchos::RCP;
00107   using Tpetra::createVector;
00108   using Teuchos::ParameterList;
00109 
00111   template <class Tout, class Tin, class LO, class GO, class Node> 
00112   struct convertHelp {
00113     static RCP<const Tpetra::CrsMatrix<Tout,LO,GO,Node>> doit(const RCP<const Tpetra::CrsMatrix<Tin,LO,GO,Node>> &A)
00114     {
00115       return A->template convert<Tout>();
00116     }
00117   };
00118 
00120   template <class T, class LO, class GO, class Node> 
00121   struct convertHelp<T,T,LO,GO,Node> {
00122     static RCP<const Tpetra::CrsMatrix<T,LO,GO,Node>> doit(const RCP<const Tpetra::CrsMatrix<T,LO,GO,Node>> &A)
00123     {
00124       return A;
00125     }
00126   };
00127 
00128 
00129   /***************************************************
00130   *   Implicit Riemannian Trust-Region Eigensolver
00131   *   Baker, Absil, Gallivan 2008
00132   ***************************************************/
00133 
00134   template <class S, class Souter, class LO, class GO, class Node>
00135   S diagPrecond(const RCP<const Tpetra::Vector<S,LO,GO,Node> > &D,
00136                 const RCP<const Tpetra::Vector<S,LO,GO,Node> > &x,
00137                 const RCP<const Tpetra::Vector<S,LO,GO,Node> > &iDx,
00138                 const S &xiDx,
00139                 const RCP<const Tpetra::Vector<S,LO,GO,Node> > &r,
00140                 const RCP<      Tpetra::Vector<S,LO,GO,Node> > &z)
00141   {
00142     /*
00143        Solve 
00144           P D P z = r
00145        for solution
00146           r = P_{u,x} inv(D) r
00147             = (I - u inv(x'*u) x') inv(D) r
00148             = z - u inv(x'*u) x'*z
00149        for u = inv(D) x = iDx
00150            z = inv(D) r
00151        
00152        for our purposes,
00153     1)        z = r
00154     2) fuse   z  \= D        
00155         and   xz  = x'*z     
00156     3) fuse   z = z - iDx ( xt / xiDx )
00157         and   return z'*r
00158     */
00159     XFORM2(z,r,  r);
00160     const S alpha = XFORM3RED(z,D,x,  z/D, x*z )  /  xiDx;
00161     return XFORM3RED(z,iDx,r,  z - iDx*alpha, z*r );
00162   }
00163 
00164   template <class S, class LO, class GO, class Node>
00165   struct tCG_workspace {
00166     RCP<Tpetra::Vector<S,LO,GO,Node>> r, z, iDx, d, Hd;
00167   };
00168 
00169   /***************************************************
00170   *   Trust-region Model Subproblem CG Solver
00171   *   Given x, compute eta that minimizes the
00172   *   quadratic Taylor expansion of f around x, 
00173   *   subject to  rho_x(eta) >= rho_prime, 
00174   *   or, equiv., eta'*eta <= ONE/rho_prime - ONE
00175   ***************************************************/
00177   template <class S, class Souter, class LO, class GO, class Node>
00178   void truncateCG(const RCP<Teuchos::FancyOStream> &out, int verbose, 
00179                   const RCP<const Tpetra::Operator<S,LO,GO,Node> > &A,
00180                   const RCP<const Tpetra::Vector<S,LO,GO,Node> > &diag,
00181                   const RCP<const Tpetra::Vector<Souter,LO,GO,Node> > &x_outer,
00182                   const RCP<const Tpetra::Vector<Souter,LO,GO,Node> > &grad,
00183                   const Souter &lx_outer, 
00184                   const RCP<Tpetra::Vector<S,LO,GO,Node> > &eta,
00185                   const tCG_workspace<S,LO,GO,Node> &workspace,
00186                   double rho_prime, double kappa, double theta, int maxItersCG, double outertol)
00187   {
00188     using Teuchos::as;
00189     using Teuchos::NO_TRANS;
00190     typedef Teuchos::ScalarTraits<S> ST;
00191     typedef Teuchos::ScalarTraits<S> STO;
00192     Teuchos::OSTab tab(out,2);
00193     const S ZERO = ST::zero();
00194     const S  ONE = ST::one();
00195     const S   D2 = ONE/rho_prime - ONE;
00196     const S lx = as<S>(lx_outer);
00197     decltype(eta) r = workspace.r,
00198                   z = workspace.z,
00199                   d = workspace.d,
00200                 iDx = workspace.iDx,
00201                  Hd = workspace.Hd;
00202     static RCP<Teuchos::Time> timer;
00203     if (timer == Teuchos::null) {
00204       std::string methName("truncatedCG<"+Teuchos::TypeNameTraits<S>::name()+">");
00205       timer = Teuchos::TimeMonitor::getNewTimer(methName);
00206     }
00207     timer->start();
00208     RCP<const Tpetra::Vector<S,LO,GO,Node>> x;
00209     {
00210       auto xinit = createVector<S>(r->getMap());
00211       XFORM2(xinit,x_outer, as<S>(x_outer));
00212       x = xinit;
00213     }
00215     // init precond. compute:
00216     //     iDx = inv(D)*x 
00217     //    xiDx = x'*iDx
00218     const S xiDx = XFORM3RED(iDx,diag,x,  x/diag,  x*iDx );
00220     // init solver
00221     XFORM1( eta,      ZERO );                            // eta = 0
00222     S r_r = XFORM2RED( r,grad,   as<S>(grad),  r*r );    // r = grad
00223     S e_e = ZERO;
00224     S z_r = diagPrecond<S,Souter,LO,GO,Node>(diag,x,iDx,xiDx,r,z);
00225     XFORM2( d,z, -z );                            // d = -z
00226     const S rnrm0 = ST::squareroot(r_r);
00227     const S kconv = rnrm0*kappa;
00228     const S tconv = ST::pow(rnrm0, (theta+1));
00229     const S futile_tol = TEUCHOS_MAX(outertol, 100*rnrm0*ST::eps());
00230     const S convtol = TEUCHOS_MAX(TEUCHOS_MIN(kconv,tconv), futile_tol);
00231     std::string convstr;
00232     if (kconv < tconv) convstr = "non-local convergence";
00233     else               convstr = "local convergence";
00235     // iterate
00236     if (verbose > 2) *out << "starting tCG with |res| =  " << rnrm0 << std::endl;
00237     if (verbose > 2) *out << "targeting tol = " << convtol << std::endl;
00238     int k = 0;
00239     std::string reason;
00240     while (true)
00241     {
00242       // apply Hessian Hd = A*d - lx*d
00243       A->apply(*d,*Hd);
00244       const S xHd = XFORM3RED(Hd,x,d,  Hd - d*lx,     x*Hd );
00245       const S dHd = XFORM3RED(Hd,x,d,  Hd - x*xHd,    d*Hd );
00246       if (verbose > 19) *out << "curvature: " << dHd << std::endl;
00247 
00248       /* 
00249          m_x(e + alpha * d) 
00250          == (e + alpha * d)'*r + 1/2 * (e + alpha * d)' * H * (e + alpha * d)
00251          == e'*r + alpha*d'*r + 1/2*e'*H*e + 1/2*alpha^2*d'*H*d
00252                    ----------                ------------------
00253          derivative w.r.t. alpha, set to zero:
00254            0 = d'*r + alpha*d'*H*d
00255          so that 
00256            alpha = -d'*r / d'*H*d
00257          substituting d = (gram-schmidt) -z
00258            alpha = z'*r / d'*H*d
00259        */
00260       // compute step size
00261       const S alpha = z_r / dHd;
00262       if (verbose > 19) *out << "alpha: " << alpha << std::endl;
00263       if (verbose > 19) *out << "z_r: " << z_r << std::endl;
00264 
00265       // anticipate e'*e for this step size
00266       const S e_d = REDUCE2(eta,d,   eta*d),
00267               d_d = REDUCE1(d,       d*d  );
00268       const S new_e_e = e_e + 2.0*alpha*e_d + alpha*alpha*d_d;
00269 
00270       // check for truncation: positive curvature or implicit trust-region boundary
00271       if (dHd <= ZERO || new_e_e > D2) 
00272       {
00273         // find truncated step to boundary
00274         // find alpha_prime 
00275         // eta_new'*eta_new = (eta + alpha_prime*d)'*(eta + alpha_prime*d)
00276         //                  = e_e + 2*alpha_prime*e_d + alpha_prime^2*d_d
00277         // choose alpha_prime > 0 such that eta_new'*eta_new == D2
00278         // solve for alpha_prime
00279         // 
00280         const S alpha_prime = (-e_d + ST::squareroot(e_d*e_d + d_d*(D2-e_e))) / d_d;
00281         if (verbose > 19) *out << "alpha_prime: " << alpha_prime << std::endl;
00282         XFORM2(eta,d,  eta + alpha_prime*d );
00283         if (dHd <= ZERO) reason="negative curvature";
00284         else             reason="exceeded trust-region";
00285         //
00286         break;
00287       }
00288 
00289       // compute new optimal point
00290       XFORM2(eta,d,  eta + alpha*d );
00291       e_e = new_e_e;    
00292 
00293       // update the residual
00294       const S x_r = XFORM3RED(r,Hd,x,   r + alpha*Hd,  x*r );
00295       // re-tangentialize r
00296       r_r = XFORM2RED(r,x,              r - x*x_r,     r*r );
00297       // compute norm
00298       S rnrm = ST::squareroot(r_r);
00299       if (verbose > 9) *out << "|r| = " << rnrm << std::endl;
00300 
00301       // check stopping criteria
00302       if (rnrm <= convtol) {reason = convstr; break;}
00303 
00304       // save the old z'*r
00305       const S zold_rold = z_r;
00306       // z = diagPrec(r)
00307       z_r = diagPrecond<S,Souter,LO,GO,Node>(diag,x,iDx,xiDx,r,z);
00308       // compute new search direction
00309       const S beta = z_r / zold_rold;
00310       XFORM2(d,z,   -z + beta*d );
00311 
00312       k += 1;
00313       if (k >= maxItersCG) {reason="maximum number of inner iterations";break;}
00314     }
00315     timer->stop();
00316     if (verbose > 1) {
00317       *out << "leaving tCG after " << k << " iterations: " << reason << std::endl;
00318     }
00319   }
00320 
00322   template <class Sinner, class S, class LO, class GO, class Node>      
00323   S IRTR(const RCP<Teuchos::FancyOStream> &out, ParameterList &params,
00324                     const RCP<const Tpetra::CrsMatrix<S,LO,GO,Node> > &A,
00325                     const RCP<Tpetra::Vector<S,LO,GO,Node> > &x)
00326   {
00327     using Teuchos::NO_TRANS;
00328     Teuchos::OSTab tab(out,2);
00329     std::string methName("IRTR<"+Teuchos::TypeNameTraits<S>::name()+","+Teuchos::TypeNameTraits<Sinner>::name()+">");
00330     typedef Tpetra::Vector<S,LO,GO,Node>     Vec;
00331     typedef Tpetra::CrsMatrix<Sinner,LO,GO,Node>  MatInner;
00332     typedef Teuchos::ScalarTraits<S>          ST;
00333     const Sinner ONE = Teuchos::ScalarTraits<Sinner>::one();
00334     // get objects from level database
00335     const int    maxIters   = params.get<int   >("maxIters",  100);
00336     const int    maxItersCG = params.get<int   >("maxItersCG",A->getDomainMap()->getGlobalNumElements());
00337     const int    verbose    = params.get<int   >("verbose",0);
00338     const double tolerance  = params.get<double>("tolerance");
00339     const double kappa      = params.get<double>("kappa",0.1);
00340     const double theta      = params.get<double>("theta",2.0);
00341     const double rho_prime  = params.get<double>("rho_prime",.1);
00342     const bool   precond    = params.get<int>   ("precond",1);
00343     if (verbose > 10) *out << params << std::endl;
00344     if (verbose > 10) *out << "maxIters: " << maxIters << std::endl;
00345     TEUCHOS_TEST_FOR_EXCEPTION(rho_prime >= 1.0, std::runtime_error, "IRTR: rho_prime must be less than 1");
00346     RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::getNewTimer(methName);
00347     // x, ax, gradfx, eta
00348     auto map = A->getDomainMap();
00349     auto Ax = createVector<S>(map),
00350           r = createVector<S>(map),
00351         eta = createVector<Sinner>(map),
00352        diag = createVector<Sinner>(map);
00353     tCG_workspace<Sinner,LO,GO,Node> workspace;
00354     workspace.r   = createVector<Sinner>(map);
00355     workspace.z   = createVector<Sinner>(map);
00356     workspace.iDx = createVector<Sinner>(map);
00357     workspace.d   = createVector<Sinner>(map);
00358     workspace.Hd  = createVector<Sinner>(map);
00359     RCP<const MatInner> Ainner = convertHelp<Sinner,S,LO,GO,Node>::doit(A);
00360     if (precond) {
00361       Ainner->getLocalDiagCopy(*diag);
00362     }
00363     else {
00364       XFORM1(diag, ONE);
00365     }
00366     timer->start(); 
00368     // init the solver
00369     S x_x = REDUCE1(x,         x*x );                           // x'*x
00370     S xnrm  = ST::squareroot(x_x);
00371     XFORM1( x,  x/xnrm);                                        // normalize x
00372     A->apply(*x,*Ax);                                           // Ax = A*x
00373     S lambda = REDUCE2(x,Ax,  x*Ax );                           // x'*A*x
00374     S rr = XFORM3RED(r,x,Ax,  Ax - x*lambda,  r*r );            // r = (I - x*x')Ax = Ax - x*lambda and compute its norm
00375     S rnrm = ST::squareroot(rr);
00376     if (verbose) {
00377       *out << "Beginning " << methName << " with λ = " << lambda << std::endl;
00378     }
00380     // iterate
00381     int k=0;
00382     while (k < maxIters && rnrm > tolerance) {
00384       // solve the trust-region subproblem for update eta
00385       truncateCG<Sinner,S,LO,GO,Node>(out,verbose,Ainner,diag,
00386                                       x,r,lambda,
00387                                       eta, workspace,
00388                                       rho_prime,kappa,theta,maxItersCG,tolerance);
00390       // update iterate, compute new eigenvalue
00391       x_x = XFORM2RED(x,eta,   x + eta,  x*x );                 // x = x+eta and x'*x
00392       xnrm = ST::squareroot(x_x);
00393       XFORM1(x,  x/xnrm);                                       // normalize x
00394       A->apply(*x,*Ax);                                         // Ax = A*x
00395       lambda = REDUCE2(x,Ax,    x*Ax );                         // lambda = x'*A*x
00396       rr = XFORM3RED(r,x,Ax,  Ax - x*lambda, r*r );             // r = A*x - x*lambda and r'*r
00397       rnrm = ST::squareroot(rr) / ST::magnitude(lambda);
00398       k += 1;
00399       if (verbose) {
00400         *out << "After iteration " << k << ", λ = " << std::setprecision(2) << std::scientific << lambda << ", |Ax-xλ|/|λ| = " << rnrm << std::endl;
00401       }
00402     }
00403     timer->stop(); 
00404     if (verbose) {
00405       *out << "Leaving " << methName << " after " << k << " iterations." << std::endl;
00406     }
00407     return lambda;
00408   }
00409 
00410 } // end of namespace TpetraExamples
00411 
00412 #endif // MULTIPRECCG_HPP_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines