11 #ifndef ROL_MINRES_HPP
12 #define ROL_MINRES_HPP
30 template<
typename Real>
47 void givens( Real& c, Real& s, Real& r, Real a, Real b )
const {
52 c = ( a >=
zero ? one : -one );
56 else if( a ==
zero ) {
58 s = ( b >=
zero ? one : -one );
61 else if( abs(a) > abs(b) ) {
63 auto u = copysign(sqrt(one+t*t),a);
70 auto u = copysign(sqrt(one+t*t),b);
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") { }
84 virtual Real
run(
V &x,
OP &A,
const V &b,
OP &M,
int &iter,
int &flag )
override {
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();
93 Real c_prev{0}, s_prev{0}, c_curr{0}, s_curr{0}, c_next{0}, s_next{0};
95 resnorm_ = v_curr->norm();
97 Real itol = sqrt(ROL_EPSILON<Real>());
99 for(
auto &e: H_ ) e = 0;
101 rhs_[0] = resnorm_; rhs_[1] = 0;
103 v_curr->scale(1.0/resnorm_);
105 for( iter=0; iter < (int)Krylov<Real>::getMaximumIteration(); iter++) {
110 if( resnorm_ < rtol )
break;
112 A.
apply( *v_next, *v_curr, itol );
114 if( iter>0 ) v_next->axpy(-H_[1],*v_prev);
116 H_[2] = v_next->dot(*v_curr);
118 v_next->axpy(-H_[2],*v_curr);
120 H_[3] = v_next->norm();
122 v_next->scale(1.0/H_[3]);
126 H_[0] = s_prev*H_[1];
132 auto tmp = c_curr*H_[2]-s_curr*H_[1];
133 H_[1] = c_curr*H_[1] + s_curr*H_[2];
138 givens( c_next, s_next, H_[2], H_[2], H_[3] );
139 rhs_[1] = -s_next*rhs_[0];
142 w_next->set( *v_curr );
144 if( iter>0 ) w_next->axpy( -H_[1], *w_curr );
145 if( iter>1 ) w_next->axpy( -H_[0], *w_prev );
147 w_next->scale(1.0/H_[2]);
149 x.
axpy( rhs_[0], *w_next );
151 v_prev->set( *v_curr );
152 v_curr->set( *v_next );
153 w_prev->set( *w_curr );
154 w_curr->set( *w_next );
165 resnorm_ = abs( rhs_[1] );
186 #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)