ROL
ROL_MINRES.hpp
Go to the documentation of this file.
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#pragma once
45#ifndef ROL_MINRES_HPP
46#define ROL_MINRES_HPP
47
48#include <array>
49#include "ROL_Krylov.hpp"
50#include "ROL_VectorClone.hpp"
51
52namespace ROL {
53
60namespace details {
61
62using namespace std;
63
64template<typename Real>
65class MINRES : public Krylov<Real> {
66
67 using V = Vector<Real>;
69
70private:
71
72 // Givens rotation matrix elements
76 array<Real,4> H_;
77 array<Real,2> rhs_;
78
80
81 void givens( Real& c, Real& s, Real& r, Real a, Real b ) const {
82
83 Real zero(0), one(1);
84
85 if( b == zero ) {
86 c = ( a >= zero ? one : -one );
87 s = zero;
88 r = abs(a);
89 }
90 else if( a == zero ) {
91 c = zero;
92 s = ( b >= zero ? one : -one );
93 r = abs(b);
94 }
95 else if( abs(a) > abs(b) ) {
96 auto t = b/a;
97 auto u = copysign(sqrt(one+t*t),a);
98 c = one/u;
99 s = c*t;
100 r = a*u;
101 }
102 else {
103 auto t = a/b;
104 auto u = copysign(sqrt(one+t*t),b);
105 s = 1/u;
106 c = s*t;
107 r = b*u;
108 }
109 } // givens()
110
111public:
112
113 MINRES( Real absTol = 1.e-4, Real relTol = 1.e-2, unsigned maxit = 100, bool useInexact = false) :
114 Krylov<Real>(absTol,relTol,maxit), useInexact_(useInexact),
115 clones_("v_prev","v_curr","v_next","w_prev","w_curr","w_next") { }
116
117 // Note: Preconditioner is not implemented
118 virtual Real run( V &x, OP &A, const V &b, OP &M, int &iter, int &flag ) override {
119
120 auto v_prev = clones_( x, "v_prev" ); v_prev->zero();
121 auto v_curr = clones_( x, "v_curr" ); v_curr->set(b);
122 auto v_next = clones_( x, "v_next" ); v_next->zero();
123 auto w_prev = clones_( x, "w_prev" ); w_prev->zero();
124 auto w_curr = clones_( x, "w_curr" ); w_curr->zero();
125 auto w_next = clones_( x, "w_next" ); w_next->zero();
126
127 Real c_prev{0}, s_prev{0}, c_curr{0}, s_curr{0}, c_next{0}, s_next{0};
128
129 resnorm_ = v_curr->norm();
131 Real itol = sqrt(ROL_EPSILON<Real>());
132
133 for( auto &e: H_ ) e = 0;
134
135 rhs_[0] = resnorm_; rhs_[1] = 0;
136
137 v_curr->scale(1.0/resnorm_);
138
139 for( iter=0; iter < (int)Krylov<Real>::getMaximumIteration(); iter++) {
140 if ( useInexact_ ) {
141 itol = rtol/((Real)Krylov<Real>::getMaximumIteration() * resnorm_);
142 }
143
144 if( resnorm_ < rtol ) break;
145
146 A.apply( *v_next, *v_curr, itol );
147
148 if( iter>0 ) v_next->axpy(-H_[1],*v_prev);
149
150 H_[2] = v_next->dot(*v_curr);
151
152 v_next->axpy(-H_[2],*v_curr);
153
154 H_[3] = v_next->norm();
155
156 v_next->scale(1.0/H_[3]);
157
158 // Rotation on H_[0] and H_[1].
159 if( iter>1 ) {
160 H_[0] = s_prev*H_[1];
161 H_[1] *= c_prev;
162 }
163
164 // Rotation on H_[1] and H_[2]
165 if( iter>0 ) {
166 auto tmp = c_curr*H_[2]-s_curr*H_[1];
167 H_[1] = c_curr*H_[1] + s_curr*H_[2];
168 H_[2] = tmp;
169 }
170
171 // Form rotation coefficients
172 givens( c_next, s_next, H_[2], H_[2], H_[3] );
173 rhs_[1] = -s_next*rhs_[0];
174 rhs_[0] *= c_next;
175
176 w_next->set( *v_curr );
177
178 if( iter>0 ) w_next->axpy( -H_[1], *w_curr );
179 if( iter>1 ) w_next->axpy( -H_[0], *w_prev );
180
181 w_next->scale(1.0/H_[2]);
182
183 x.axpy( rhs_[0], *w_next );
184
185 v_prev->set( *v_curr );
186 v_curr->set( *v_next );
187 w_prev->set( *w_curr );
188 w_curr->set( *w_next );
189
190 c_prev = c_curr;
191 c_curr = c_next;
192 s_prev = s_curr;
193 s_curr = s_next;
194
195 rhs_[0] = rhs_[1];
196
197 H_[1] = H_[3];
198
199 resnorm_ = abs( rhs_[1] );
200
201 } // for (iter)
202
203 if ( iter == (int)Krylov<Real>::getMaximumIteration() ) flag = 1;
204 else iter++;
205
206 return resnorm_;
207 } // run()
208
209}; // class MINRES
210
211} // namespace details
212
213
214using details::MINRES;
215
216
217} // namespace ROL
218
219
220#endif // ROL_MINRES_HPP
221
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:58
Provides the interface to apply a linear operator.
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const =0
Apply linear operator.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
VectorCloneMap< Real > clones_
Definition: ROL_MINRES.hpp:79
array< Real, 4 > H_
Definition: ROL_MINRES.hpp:76
MINRES(Real absTol=1.e-4, Real relTol=1.e-2, unsigned maxit=100, bool useInexact=false)
Definition: ROL_MINRES.hpp:113
virtual Real run(V &x, OP &A, const V &b, OP &M, int &iter, int &flag) override
Definition: ROL_MINRES.hpp:118
void givens(Real &c, Real &s, Real &r, Real a, Real b) const
Definition: ROL_MINRES.hpp:81
array< Real, 2 > rhs_
Definition: ROL_MINRES.hpp:77