ROL
ROL_GMRES.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_GMRES_H
11 #define ROL_GMRES_H
12 
17 #include "ROL_Krylov.hpp"
18 #include "ROL_Types.hpp"
19 #include "ROL_LAPACK.hpp"
20 #include "ROL_LinearAlgebra.hpp"
21 
22 
23 namespace ROL {
24 
25 template<class Real>
26 class GMRES : public Krylov<Real> {
27 
28  typedef LA::Matrix<Real> SDMatrix;
29  typedef LA::Vector<Real> SDVector;
30 
31 private:
32 
33  ROL::Ptr<Vector<Real> > r_;
34  ROL::Ptr<Vector<Real> > z_;
35  ROL::Ptr<Vector<Real> > w_;
36 
37  ROL::Ptr<SDMatrix> H_; // quasi-Hessenberg matrix
38  ROL::Ptr<SDVector> cs_; // Givens Rotations cosine components
39  ROL::Ptr<SDVector> sn_; // Givens Rotations sine components
40  ROL::Ptr<SDVector> s_;
41  ROL::Ptr<SDVector> y_;
42  ROL::Ptr<SDVector> cnorm_;
43 
44  ROL::Ptr<std::vector<Real> > res_;
45 
48  bool useInitialGuess_; // If false, inital x will be ignored and zero vec used
50  ROL::Ptr<std::ostream> outStream_;
51 
52  ROL::LAPACK<int,Real> lapack_;
53 
54 public:
55 
56  GMRES( ROL::ParameterList &parlist ) : Krylov<Real>(parlist), isInitialized_(false), printIters_(false) {
57 
58  using std::vector;
59 
60  Real zero(0);
61 
62  ROL::ParameterList &gList = parlist.sublist("General");
63  ROL::ParameterList &kList = gList.sublist("Krylov");
64 
65  useInexact_ = gList.get("Inexact Hessian-Times-A-Vector",false);
66  useInitialGuess_ = kList.get("Use Initial Guess",false);
68 
69  H_ = ROL::makePtr<SDMatrix>( maxit+1, maxit );
70  cs_ = ROL::makePtr<SDVector>( maxit );
71  sn_ = ROL::makePtr<SDVector>( maxit );
72  s_ = ROL::makePtr<SDVector>( maxit+1 );
73  y_ = ROL::makePtr<SDVector>( maxit+1 );
74  cnorm_ = ROL::makePtr<SDVector>( maxit );
75  res_ = ROL::makePtr<std::vector<Real>>(maxit+1,zero);
76 
77  }
78 
80  LinearOperator<Real> &M, int &iter, int &flag ) {
81 
82  Real absTol = Krylov<Real>::getAbsoluteTolerance();
83  Real relTol = Krylov<Real>::getRelativeTolerance();
85 
86  flag = 0;
87 
88  Real zero(0), one(1);
89 
90  if ( !isInitialized_ ) {
91  r_ = b.clone();
92  w_ = b.clone();
93  z_ = x.clone();
94 
95  isInitialized_ = true;
96  }
97 
98  Real itol = std::sqrt(ROL_EPSILON<Real>());
99 
100  // Compute initial residual
101  if(useInitialGuess_) {
102 
103  A.apply(*r_,x,itol);
104  r_->scale(-one);
105  r_->plus(b); // r = b-Ax
106 
107  }
108  else {
109  x.zero();
110  r_->set(b);
111  }
112 
113  Real temp = 0;
114 
115  std::vector<ROL::Ptr<Vector<Real > > > V;
116  std::vector<ROL::Ptr<Vector<Real > > > Z;
117 
118  (*res_)[0] = r_->norm();
119 
120  if (printIters_) {
121  *outStream_ << "GMRES Iteration " << 0 << ", Residual = " << (*res_)[0] << "\n";
122  }
123 
124  // This should be a tolerance check
125  Real rtol = std::min(absTol,relTol*(*res_)[0]);
126  if ((*res_)[0] <= rtol) {
127  iter = 0;
128  flag = 0;
129  return (*res_)[0];
130  }
131 
132  V.push_back(b.clone());
133  (V[0])->set(*r_);
134  (V[0])->scale(one/(*res_)[0]);
135 
136  (*s_)(0) = (*res_)[0];
137 
138  for( iter=0; iter<maxit; ++iter ) {
139 
140  if( useInexact_ ) {
141  itol = rtol/(maxit*(*res_)[iter]);
142  }
143 
144  Z.push_back(x.clone());
145 
146  // Apply right preconditioner
147  M.applyInverse(*(Z[iter]),*(V[iter]),itol);
148 
149  // Apply operator
150  A.apply(*w_,*(Z[iter]),itol);
151 
152  // Evaluate coefficients and orthogonalize using Gram-Schmidt
153  for( int k=0; k<=iter; ++k ) {
154  (*H_)(k,iter) = w_->dot(*(V[k]));
155  w_->axpy( -(*H_)(k,iter), *(V[k]) );
156  }
157 
158  (*H_)(iter+1,iter) = w_->norm();
159 
160  V.push_back( b.clone() );
161  (V[iter+1])->set(*w_);
162  (V[iter+1])->scale(one/((*H_)(iter+1,iter)));
163 
164  // Apply Givens rotations
165  for( int k=0; k<=iter-1; ++k ) {
166  temp = (*cs_)(k)*(*H_)(k,iter) + (*sn_)(k)*(*H_)(k+1,iter);
167  (*H_)(k+1,iter) = -(*sn_)(k)*(*H_)(k,iter) + (*cs_)(k)*(*H_)(k+1,iter);
168  (*H_)(k,iter) = temp;
169  }
170 
171  // Form i-th rotation matrix
172  if( (*H_)(iter+1,iter) == zero ) {
173  (*cs_)(iter) = one;
174  (*sn_)(iter) = zero;
175  }
176  else if ( std::abs((*H_)(iter+1,iter)) > std::abs((*H_)(iter,iter)) ) {
177  temp = (*H_)(iter,iter) / (*H_)(iter+1,iter);
178  (*sn_)(iter) = one / std::sqrt( one + temp*temp );
179  (*cs_)(iter) = temp*(*sn_)(iter);
180  }
181  else {
182  temp = (*H_)(iter+1,iter) / (*H_)(iter,iter);
183  (*cs_)(iter) = one / std::sqrt( one + temp*temp );
184  (*sn_)(iter) = temp*(*cs_)(iter);
185  }
186 
187  // Approximate residual norm
188  temp = (*cs_)(iter)*(*s_)(iter);
189  (*s_)(iter+1) = -(*sn_)(iter)*(*s_)(iter);
190  (*s_)(iter) = temp;
191  (*H_)(iter,iter) = (*cs_)(iter)*(*H_)(iter,iter) + (*sn_)(iter)*(*H_)(iter+1,iter);
192  (*H_)(iter+1,iter) = zero;
193  (*res_)[iter+1] = std::abs((*s_)(iter+1));
194 
195  if (printIters_) {
196  *outStream_ << "GMRES Iteration " << iter+1 << ", Residual = " << (*res_)[iter+1] << "\n";
197  }
198 
199  // Update solution approximation.
200  const char uplo = 'U';
201  const char trans = 'N';
202  const char diag = 'N';
203  const char normin = 'N';
204  Real scaling = zero;
205  int info = 0;
206  *y_ = *s_;
207  lapack_.LATRS(uplo, trans, diag, normin, iter+1, H_->values(), maxit+1, y_->values(), &scaling, cnorm_->values(), &info);
208 
209  z_->zero();
210 
211  for( int k=0; k<=iter;++k ) {
212  z_->axpy((*y_)(k),*(Z[k]));
213  }
214 
215  if( (*res_)[iter+1] <= rtol ) {
216  // Update solution vector
217  x.plus(*z_);
218  break;
219  }
220 
221  } // loop over iter
222 
223  if(iter == maxit) {
224  flag = 1;
225  x.plus(*z_);
226  return (*res_)[iter];
227  }
228 
229  return (*res_)[iter+1];
230  }
231 
232  void enableOutput(std::ostream & outStream) {
233  printIters_ = true;
234  outStream_ = ROL::makePtrFromRef(outStream);;
235  }
236 
237  void disableOutput() {printIters_ = false;}
238 
239 }; // class GMRES
240 
241 } // namespace ROL
242 
243 #endif // ROL_GMRES_H
244 
ROL::Ptr< SDVector > cs_
Definition: ROL_GMRES.hpp:38
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Real getRelativeTolerance(void) const
Definition: ROL_Krylov.hpp:60
Real getAbsoluteTolerance(void) const
Definition: ROL_Krylov.hpp:57
bool isInitialized_
Definition: ROL_GMRES.hpp:46
virtual void plus(const Vector &x)=0
Compute , where .
GMRES(ROL::ParameterList &parlist)
Definition: ROL_GMRES.hpp:56
bool useInexact_
Definition: ROL_GMRES.hpp:47
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.
unsigned getMaximumIteration(void) const
Definition: ROL_Krylov.hpp:63
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_
Definition: ROL_GMRES.hpp:33
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Vector< Real > V
bool printIters_
Definition: ROL_GMRES.hpp:49
ROL::Ptr< std::ostream > outStream_
Definition: ROL_GMRES.hpp:50
LA::Matrix< Real > SDMatrix
Definition: ROL_GMRES.hpp:28
Real run(Vector< Real > &x, LinearOperator< Real > &A, const Vector< Real > &b, LinearOperator< Real > &M, int &iter, int &flag)
Definition: ROL_GMRES.hpp:79
ROL::Ptr< std::vector< Real > > res_
Definition: ROL_GMRES.hpp:44
Preconditioned GMRES solver.
Definition: ROL_GMRES.hpp:26
void disableOutput()
Definition: ROL_GMRES.hpp:237
bool useInitialGuess_
Definition: ROL_GMRES.hpp:48
virtual void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
ROL::Ptr< SDVector > y_
Definition: ROL_GMRES.hpp:41
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:24
Provides the interface to apply a linear operator.
ROL::Ptr< SDVector > cnorm_
Definition: ROL_GMRES.hpp:42
void enableOutput(std::ostream &outStream)
Definition: ROL_GMRES.hpp:232
ROL::Ptr< SDVector > sn_
Definition: ROL_GMRES.hpp:39
ROL::Ptr< Vector< Real > > w_
Definition: ROL_GMRES.hpp:35
ROL::Ptr< SDMatrix > H_
Definition: ROL_GMRES.hpp:37
LA::Vector< Real > SDVector
Definition: ROL_GMRES.hpp:29
ROL::Ptr< Vector< Real > > z_
Definition: ROL_GMRES.hpp:34
ROL::Ptr< SDVector > s_
Definition: ROL_GMRES.hpp:40
ROL::LAPACK< int, Real > lapack_
Definition: ROL_GMRES.hpp:52