ROL
ROL_ConjugateGradients.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_CONJUGATEGRADIENTS_H
11 #define ROL_CONJUGATEGRADIENTS_H
12 
17 #include "ROL_Krylov.hpp"
18 #include "ROL_Types.hpp"
19 
20 namespace ROL {
21 
22 template<class Real>
23 class ConjugateGradients : public Krylov<Real> {
24 
27  ROL::Ptr<Vector<Real> > r_;
28  ROL::Ptr<Vector<Real> > v_;
29  ROL::Ptr<Vector<Real> > p_;
30  ROL::Ptr<Vector<Real> > Ap_;
31 
32 public:
33  ConjugateGradients(Real absTol = 1.e-4, Real relTol = 1.e-2, unsigned maxit = 100, bool useInexact = false)
34  : Krylov<Real>(absTol,relTol,maxit), isInitialized_(false), useInexact_(useInexact) {}
35 
37  int &iter, int &flag ) {
38  if ( !isInitialized_ ) {
39  r_ = b.clone();
40  v_ = x.clone();
41  p_ = x.clone();
42  Ap_ = b.clone();
43  isInitialized_ = true;
44  }
45 
46  Real rnorm = b.norm();
48  Real itol = std::sqrt(ROL_EPSILON<Real>());
49 
50  x.zero();
51  r_->set(b);
52 
53  M.applyInverse(*v_, *r_, itol);
54  p_->set(*v_);
55 
56  iter = 0;
57  flag = 0;
58 
59  Real kappa(0), beta(0), alpha(0), tmp(0), zero(0);
60  //Real gv = v_->dot(r_->dual());
61  Real gv = v_->apply(*r_);
62 
63  for (iter = 0; iter < (int)Krylov<Real>::getMaximumIteration(); iter++) {
64  if ( useInexact_ ) {
65  itol = rtol/((Real)Krylov<Real>::getMaximumIteration() * rnorm);
66  }
67  A.apply(*Ap_, *p_, itol);
68 
69  //kappa = p_->dot(Ap_->dual());
70  kappa = p_->apply(*Ap_);
71  if ( kappa <= zero ) {
72  flag = 2;
73  break;
74  }
75  alpha = gv/kappa;
76 
77  x.axpy(alpha,*p_);
78 
79  r_->axpy(-alpha,*Ap_);
80  rnorm = r_->norm();
81  if ( rnorm < rtol ) {
82  break;
83  }
84 
85  itol = std::sqrt(ROL_EPSILON<Real>());
86  M.applyInverse(*v_, *r_, itol);
87  tmp = gv;
88  //gv = v_->dot(r_->dual());
89  gv = v_->apply(*r_);
90  beta = gv/tmp;
91 
92  p_->scale(beta);
93  p_->plus(*v_);
94  }
95  if ( iter == (int)Krylov<Real>::getMaximumIteration() ) {
96  flag = 1;
97  }
98  else {
99  iter++;
100  }
101  return rnorm;
102  }
103 };
104 
105 
106 }
107 
108 #endif
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
Contains definitions of custom data types in ROL.
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const =0
Apply linear operator.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
ROL::Ptr< Vector< Real > > v_
virtual void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
Provides definitions of the Conjugate Gradient solver.
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:24
Provides the interface to apply a linear operator.
ConjugateGradients(Real absTol=1.e-4, Real relTol=1.e-2, unsigned maxit=100, bool useInexact=false)
ROL::Ptr< Vector< Real > > p_
ROL::Ptr< Vector< Real > > r_
virtual Real norm() const =0
Returns where .
ROL::Ptr< Vector< Real > > Ap_
Real run(Vector< Real > &x, LinearOperator< Real > &A, const Vector< Real > &b, LinearOperator< Real > &M, int &iter, int &flag)