ROL
ROL_MINRES.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 #pragma once
11 #ifndef ROL_MINRES_HPP
12 #define ROL_MINRES_HPP
13 
14 #include <array>
15 #include "ROL_Krylov.hpp"
16 #include "ROL_VectorClone.hpp"
17 
18 namespace ROL {
19 
26 namespace details {
27 
28 using namespace std;
29 
30 template<typename Real>
31 class MINRES : public Krylov<Real> {
32 
33  using V = Vector<Real>;
35 
36 private:
37 
38  // Givens rotation matrix elements
39  Real resnorm_;
40  int maxiter_;
42  array<Real,4> H_;
43  array<Real,2> rhs_;
44 
46 
47  void givens( Real& c, Real& s, Real& r, Real a, Real b ) const {
48 
49  Real zero(0), one(1);
50 
51  if( b == zero ) {
52  c = ( a >= zero ? one : -one );
53  s = zero;
54  r = abs(a);
55  }
56  else if( a == zero ) {
57  c = zero;
58  s = ( b >= zero ? one : -one );
59  r = abs(b);
60  }
61  else if( abs(a) > abs(b) ) {
62  auto t = b/a;
63  auto u = copysign(sqrt(one+t*t),a);
64  c = one/u;
65  s = c*t;
66  r = a*u;
67  }
68  else {
69  auto t = a/b;
70  auto u = copysign(sqrt(one+t*t),b);
71  s = 1/u;
72  c = s*t;
73  r = b*u;
74  }
75  } // givens()
76 
77 public:
78 
79  MINRES( Real absTol = 1.e-4, Real relTol = 1.e-2, unsigned maxit = 100, bool useInexact = false) :
80  Krylov<Real>(absTol,relTol,maxit), useInexact_(useInexact),
81  clones_("v_prev","v_curr","v_next","w_prev","w_curr","w_next") { }
82 
83  // Note: Preconditioner is not implemented
84  virtual Real run( V &x, OP &A, const V &b, OP &M, int &iter, int &flag ) override {
85 
86  auto v_prev = clones_( x, "v_prev" ); v_prev->zero();
87  auto v_curr = clones_( x, "v_curr" ); v_curr->set(b);
88  auto v_next = clones_( x, "v_next" ); v_next->zero();
89  auto w_prev = clones_( x, "w_prev" ); w_prev->zero();
90  auto w_curr = clones_( x, "w_curr" ); w_curr->zero();
91  auto w_next = clones_( x, "w_next" ); w_next->zero();
92 
93  Real c_prev{0}, s_prev{0}, c_curr{0}, s_curr{0}, c_next{0}, s_next{0};
94 
95  resnorm_ = v_curr->norm();
97  Real itol = sqrt(ROL_EPSILON<Real>());
98 
99  for( auto &e: H_ ) e = 0;
100 
101  rhs_[0] = resnorm_; rhs_[1] = 0;
102 
103  v_curr->scale(1.0/resnorm_);
104 
105  for( iter=0; iter < (int)Krylov<Real>::getMaximumIteration(); iter++) {
106  if ( useInexact_ ) {
107  itol = rtol/((Real)Krylov<Real>::getMaximumIteration() * resnorm_);
108  }
109 
110  if( resnorm_ < rtol ) break;
111 
112  A.apply( *v_next, *v_curr, itol );
113 
114  if( iter>0 ) v_next->axpy(-H_[1],*v_prev);
115 
116  H_[2] = v_next->dot(*v_curr);
117 
118  v_next->axpy(-H_[2],*v_curr);
119 
120  H_[3] = v_next->norm();
121 
122  v_next->scale(1.0/H_[3]);
123 
124  // Rotation on H_[0] and H_[1].
125  if( iter>1 ) {
126  H_[0] = s_prev*H_[1];
127  H_[1] *= c_prev;
128  }
129 
130  // Rotation on H_[1] and H_[2]
131  if( iter>0 ) {
132  auto tmp = c_curr*H_[2]-s_curr*H_[1];
133  H_[1] = c_curr*H_[1] + s_curr*H_[2];
134  H_[2] = tmp;
135  }
136 
137  // Form rotation coefficients
138  givens( c_next, s_next, H_[2], H_[2], H_[3] );
139  rhs_[1] = -s_next*rhs_[0];
140  rhs_[0] *= c_next;
141 
142  w_next->set( *v_curr );
143 
144  if( iter>0 ) w_next->axpy( -H_[1], *w_curr );
145  if( iter>1 ) w_next->axpy( -H_[0], *w_prev );
146 
147  w_next->scale(1.0/H_[2]);
148 
149  x.axpy( rhs_[0], *w_next );
150 
151  v_prev->set( *v_curr );
152  v_curr->set( *v_next );
153  w_prev->set( *w_curr );
154  w_curr->set( *w_next );
155 
156  c_prev = c_curr;
157  c_curr = c_next;
158  s_prev = s_curr;
159  s_curr = s_next;
160 
161  rhs_[0] = rhs_[1];
162 
163  H_[1] = H_[3];
164 
165  resnorm_ = abs( rhs_[1] );
166 
167  } // for (iter)
168 
169  if ( iter == (int)Krylov<Real>::getMaximumIteration() ) flag = 1;
170  else iter++;
171 
172  return resnorm_;
173  } // run()
174 
175 }; // class MINRES
176 
177 } // namespace details
178 
179 
180 using details::MINRES;
181 
182 
183 } // namespace ROL
184 
185 
186 #endif // ROL_MINRES_HPP
187 
virtual Real run(V &x, OP &A, const V &b, OP &M, int &iter, int &flag) override
Definition: ROL_MINRES.hpp:84
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const =0
Apply linear operator.
array< Real, 4 > H_
Definition: ROL_MINRES.hpp:42
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
array< Real, 2 > rhs_
Definition: ROL_MINRES.hpp:43
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:24
Provides the interface to apply a linear operator.
void givens(Real &c, Real &s, Real &r, Real a, Real b) const
Definition: ROL_MINRES.hpp:47
VectorCloneMap< Real > clones_
Definition: ROL_MINRES.hpp:45
Implements the MINRES algorithm for solving symmetric indefinite systems.
MINRES(Real absTol=1.e-4, Real relTol=1.e-2, unsigned maxit=100, bool useInexact=false)
Definition: ROL_MINRES.hpp:79