45 #ifndef ROL_MINRES_HPP
46 #define ROL_MINRES_HPP
65 template<
typename Real>
82 void givens( Real& c, Real& s, Real& r, Real a, Real b )
const {
87 c = ( a >=
zero ? one : -one );
91 else if( a ==
zero ) {
93 s = ( b >=
zero ? one : -one );
96 else if( abs(a) > abs(b) ) {
98 auto u = copysign(sqrt(one+t*t),a);
105 auto u = copysign(sqrt(one+t*t),b);
114 MINRES( Real absTol = 1.e-4, Real relTol = 1.e-2,
unsigned maxit = 100,
bool useInexact =
false) :
115 Krylov<Real>(absTol,relTol,maxit), useInexact_(useInexact),
116 clones_(
"v_prev",
"v_curr",
"v_next",
"w_prev",
"w_curr",
"w_next") { }
119 virtual Real
run(
V &x,
OP &A,
const V &b,
OP &M,
int &iter,
int &flag )
override {
121 auto v_prev = clones_( x,
"v_prev" ); v_prev->zero();
122 auto v_curr = clones_( x,
"v_curr" ); v_curr->set(b);
123 auto v_next = clones_( x,
"v_next" ); v_next->zero();
124 auto w_prev = clones_( x,
"w_prev" ); w_prev->zero();
125 auto w_curr = clones_( x,
"w_curr" ); w_curr->zero();
126 auto w_next = clones_( x,
"w_next" ); w_next->zero();
128 Real c_prev{0}, s_prev{0}, c_curr{0}, s_curr{0}, c_next{0}, s_next{0};
130 resnorm_ = v_curr->norm();
132 Real itol = sqrt(ROL_EPSILON<Real>());
134 for(
auto &e: H_ ) e = 0;
136 rhs_[0] = resnorm_; rhs_[1] = 0;
138 v_curr->scale(1.0/resnorm_);
140 for( iter=0; iter < (int)Krylov<Real>::getMaximumIteration(); iter++) {
145 if( resnorm_ < rtol )
break;
147 A.
apply( *v_next, *v_curr, itol );
149 if( iter>0 ) v_next->axpy(-H_[1],*v_prev);
151 H_[2] = v_next->dot(*v_curr);
153 v_next->axpy(-H_[2],*v_curr);
155 H_[3] = v_next->norm();
157 v_next->scale(1.0/H_[3]);
161 H_[0] = s_prev*H_[1];
167 auto tmp = c_curr*H_[2]-s_curr*H_[1];
168 H_[1] = c_curr*H_[1] + s_curr*H_[2];
173 givens( c_next, s_next, H_[2], H_[2], H_[3] );
174 rhs_[1] = -s_next*rhs_[0];
177 w_next->set( *v_curr );
179 if( iter>0 ) w_next->axpy( -H_[1], *w_curr );
180 if( iter>1 ) w_next->axpy( -H_[0], *w_prev );
182 w_next->scale(1.0/H_[2]);
184 x.
axpy( rhs_[0], *w_next );
186 v_prev->set( *v_curr );
187 v_curr->set( *v_next );
188 w_prev->set( *w_curr );
189 w_curr->set( *w_next );
200 resnorm_ = abs( rhs_[1] );
221 #endif // ROL_MINRES_HPP
virtual Real run(V &x, OP &A, const V &b, OP &M, int &iter, int &flag) override
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
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.
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides definitions for Krylov solvers.
Provides the interface to apply a linear operator.
void givens(Real &c, Real &s, Real &r, Real a, Real b) const
VectorCloneMap< Real > clones_
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)