1 // @HEADER 2 // ************************************************************************ 3 // 4 // Rapid Optimization Library (ROL) Package 5 // Copyright (2014) Sandia Corporation 6 // 7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 8 // license for use of this work by or on behalf of the U.S. Government. 9 // 10 // Redistribution and use in source and binary forms, with or without 11 // modification, are permitted provided that the following conditions are 12 // met: 13 // 14 // 1. Redistributions of source code must retain the above copyright 15 // notice, this list of conditions and the following disclaimer. 16 // 17 // 2. Redistributions in binary form must reproduce the above copyright 18 // notice, this list of conditions and the following disclaimer in the 19 // documentation and/or other materials provided with the distribution. 20 // 21 // 3. Neither the name of the Corporation nor the names of the 22 // contributors may be used to endorse or promote products derived from 23 // this software without specific prior written permission. 24 // 25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 // 37 // Questions? Contact lead developers: 38 // Drew Kouri (dpkouri@sandia.gov) and 39 // Denis Ridzal (dridzal@sandia.gov) 40 // 41 // ************************************************************************ 42 // @HEADER 43 44 #ifndef ROL_CONJUGATEGRADIENTS_H 45 #define ROL_CONJUGATEGRADIENTS_H 46 47 /** \class ROL::ConjugateGradients 48 \brief Provides definitions of the Conjugate Gradient solver. 49 */ 50 51 #include "ROL_Krylov.hpp" 52 #include "ROL_Types.hpp" 53 54 namespace ROL { 55 56 template<class Real> 57 class ConjugateGradients : public Krylov<Real> { 58 59 bool isInitialized_; 60 bool useInexact_; 61 ROL::Ptr<Vector<Real> > r_; 62 ROL::Ptr<Vector<Real> > v_; 63 ROL::Ptr<Vector<Real> > p_; 64 ROL::Ptr<Vector<Real> > Ap_; 65 66 public: ConjugateGradients(Real absTol=1.e-4,Real relTol=1.e-2,unsigned maxit=100,bool useInexact=false)67 ConjugateGradients(Real absTol = 1.e-4, Real relTol = 1.e-2, unsigned maxit = 100, bool useInexact = false) 68 : Krylov<Real>(absTol,relTol,maxit), isInitialized_(false), useInexact_(useInexact) {} 69 run(Vector<Real> & x,LinearOperator<Real> & A,const Vector<Real> & b,LinearOperator<Real> & M,int & iter,int & flag)70 Real run( Vector<Real> &x, LinearOperator<Real> &A, const Vector<Real> &b, LinearOperator<Real> &M, 71 int &iter, int &flag ) { 72 if ( !isInitialized_ ) { 73 r_ = b.clone(); 74 v_ = x.clone(); 75 p_ = x.clone(); 76 Ap_ = b.clone(); 77 isInitialized_ = true; 78 } 79 80 Real rnorm = b.norm(); 81 Real rtol = std::min(Krylov<Real>::getAbsoluteTolerance(),Krylov<Real>::getRelativeTolerance()*rnorm); 82 Real itol = std::sqrt(ROL_EPSILON<Real>()); 83 84 x.zero(); 85 r_->set(b); 86 87 M.applyInverse(*v_, *r_, itol); 88 p_->set(*v_); 89 90 iter = 0; 91 flag = 0; 92 93 Real kappa(0), beta(0), alpha(0), tmp(0), zero(0); 94 Real gv = v_->dot(r_->dual()); 95 96 for (iter = 0; iter < (int)Krylov<Real>::getMaximumIteration(); iter++) { 97 if ( useInexact_ ) { 98 itol = rtol/((Real)Krylov<Real>::getMaximumIteration() * rnorm); 99 } 100 A.apply(*Ap_, *p_, itol); 101 102 kappa = p_->dot(Ap_->dual()); 103 if ( kappa <= zero ) { 104 flag = 2; 105 break; 106 } 107 alpha = gv/kappa; 108 109 x.axpy(alpha,*p_); 110 111 r_->axpy(-alpha,*Ap_); 112 rnorm = r_->norm(); 113 if ( rnorm < rtol ) { 114 break; 115 } 116 117 itol = std::sqrt(ROL_EPSILON<Real>()); 118 M.applyInverse(*v_, *r_, itol); 119 tmp = gv; 120 gv = v_->dot(r_->dual()); 121 beta = gv/tmp; 122 123 p_->scale(beta); 124 p_->plus(*v_); 125 } 126 if ( iter == (int)Krylov<Real>::getMaximumIteration() ) { 127 flag = 1; 128 } 129 else { 130 iter++; 131 } 132 return rnorm; 133 } 134 }; 135 136 137 } 138 139 #endif 140