ROL
ROL_ConjugateResiduals.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_CONJUGATERESIDUALS_H
11 #define ROL_CONJUGATERESIDUALS_H
12 
17 #include "ROL_Krylov.hpp"
18 #include "ROL_Types.hpp"
19 
20 namespace ROL {
21 
22 template<class Real>
23 class ConjugateResiduals : 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  ROL::Ptr<Vector<Real> > MAp_;
32 
33 public:
34  ConjugateResiduals( Real absTol = 1.e-4, Real relTol = 1.e-2, int maxit = 100, bool useInexact = false )
35  : Krylov<Real>(absTol,relTol,maxit), isInitialized_(false), useInexact_(useInexact) {}
36 
37  // Run Krylov Method
39  int &iter, int &flag ) {
40  if ( !isInitialized_ ) {
41  r_ = x.clone();
42  v_ = b.clone();
43  p_ = x.clone();
44  Ap_ = b.clone();
45  MAp_ = x.clone();
46  isInitialized_ = true;
47  }
48 
49  // Initialize
50  Real rnorm = b.norm();
52  Real itol = std::sqrt(ROL_EPSILON<Real>());
53  x.zero();
54 
55  // Apply preconditioner to residual
56  M.applyInverse(*r_,b,itol);
57 
58  // Initialize direction p
59  p_->set(*r_);
60 
61  // Get Hessian tolerance
62  if ( useInexact_ ) {
63  itol = rtol/((Real)Krylov<Real>::getMaximumIteration() * rnorm);
64  }
65 
66  // Apply Hessian to residual
67  A.apply(*v_, *r_, itol);
68 
69  // Apply Hessian to direction p
70  //A.apply(*Ap_, *p_, itol);
71  Ap_->set(*v_);
72 
73  // Initialize scalar quantities
74  iter = 0;
75  flag = 0;
76  Real kappa(0), beta(0), alpha(0), tmp(0);
77  //Real gHg = r_->dot(v_->dual());
78  Real gHg = r_->apply(*v_);
79 
80  for (iter = 0; iter < (int)Krylov<Real>::getMaximumIteration(); iter++) {
81  itol = std::sqrt(ROL_EPSILON<Real>());
82  M.applyInverse(*MAp_, *Ap_, itol);
83  //kappa = MAp_->dot(Ap_->dual());
84  kappa = MAp_->apply(*Ap_);
85  //if ( gHg <= 0.0 || kappa <= 0.0 ) {
86  //flag = 2;
87  //break;
88  //}
89  alpha = gHg/kappa;
90 
91  x.axpy(alpha,*p_);
92 
93  r_->axpy(-alpha,*MAp_);
94  rnorm = r_->norm();
95  if ( rnorm < rtol ) {
96  break;
97  }
98 
99  if ( useInexact_ ) {
100  itol = rtol/((Real)Krylov<Real>::getMaximumIteration() * rnorm);
101  }
102  A.apply(*v_, *r_, itol);
103  tmp = gHg;
104  //gHg = r_->dot(v_->dual());
105  gHg = r_->apply(*v_);
106  beta = gHg/tmp;
107 
108  p_->scale(beta);
109  p_->plus(*r_);
110 
111  Ap_->scale(beta);
112  Ap_->plus(*v_);
113  }
114  if ( iter == (int)Krylov<Real>::getMaximumIteration() ) {
115  flag = 1;
116  }
117  else {
118  iter++;
119  }
120  return rnorm;
121  }
122 };
123 
124 
125 }
126 
127 #endif
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Provides definition of the Conjugate Residual solver.
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
ROL::Ptr< Vector< Real > > r_
ROL::Ptr< Vector< Real > > Ap_
virtual void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
Real run(Vector< Real > &x, LinearOperator< Real > &A, const Vector< Real > &b, LinearOperator< Real > &M, int &iter, int &flag)
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:24
Provides the interface to apply a linear operator.
ROL::Ptr< Vector< Real > > v_
ConjugateResiduals(Real absTol=1.e-4, Real relTol=1.e-2, int maxit=100, bool useInexact=false)
ROL::Ptr< Vector< Real > > p_
ROL::Ptr< Vector< Real > > MAp_
virtual Real norm() const =0
Returns where .