19 #include "ROL_LAPACK.hpp"
20 #include "ROL_LinearAlgebra.hpp"
33 ROL::Ptr<Vector<Real> >
r_;
34 ROL::Ptr<Vector<Real> >
z_;
35 ROL::Ptr<Vector<Real> >
w_;
37 ROL::Ptr<SDMatrix>
H_;
38 ROL::Ptr<SDVector>
cs_;
39 ROL::Ptr<SDVector>
sn_;
40 ROL::Ptr<SDVector>
s_;
41 ROL::Ptr<SDVector>
y_;
44 ROL::Ptr<std::vector<Real> >
res_;
62 ROL::ParameterList &gList = parlist.sublist(
"General");
63 ROL::ParameterList &kList = gList.sublist(
"Krylov");
65 useInexact_ = gList.get(
"Inexact Hessian-Times-A-Vector",
false);
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);
98 Real itol = std::sqrt(ROL_EPSILON<Real>());
115 std::vector<ROL::Ptr<Vector<Real > > >
V;
116 std::vector<ROL::Ptr<Vector<Real > > > Z;
118 (*res_)[0] =
r_->norm();
121 *
outStream_ <<
"GMRES Iteration " << 0 <<
", Residual = " << (*res_)[0] <<
"\n";
125 Real rtol = std::min(absTol,relTol*(*
res_)[0]);
126 if ((*
res_)[0] <= rtol) {
132 V.push_back(b.
clone());
134 (V[0])->scale(one/(*
res_)[0]);
136 (*s_)(0) = (*
res_)[0];
138 for( iter=0; iter<maxit; ++iter ) {
141 itol = rtol/(maxit*(*res_)[iter]);
144 Z.push_back(x.
clone());
153 for(
int k=0; k<=iter; ++k ) {
154 (*H_)(k,iter) =
w_->dot(*(V[k]));
155 w_->axpy( -(*
H_)(k,iter), *(V[k]) );
158 (*H_)(iter+1,iter) =
w_->norm();
160 V.push_back( b.
clone() );
161 (V[iter+1])->set(*
w_);
162 (V[iter+1])->scale(one/((*
H_)(iter+1,iter)));
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;
172 if( (*
H_)(iter+1,iter) ==
zero ) {
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);
182 temp = (*H_)(iter+1,iter) / (*
H_)(iter,iter);
183 (*cs_)(iter) = one / std::sqrt( one + temp*temp );
184 (*sn_)(iter) = temp*(*
cs_)(iter);
188 temp = (*cs_)(iter)*(*
s_)(iter);
189 (*s_)(iter+1) = -(*
sn_)(iter)*(*
s_)(iter);
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));
196 *
outStream_ <<
"GMRES Iteration " << iter+1 <<
", Residual = " << (*res_)[iter+1] <<
"\n";
200 const char uplo =
'U';
201 const char trans =
'N';
202 const char diag =
'N';
203 const char normin =
'N';
207 lapack_.LATRS(uplo, trans, diag, normin, iter+1,
H_->values(), maxit+1,
y_->values(), &scaling,
cnorm_->values(), &info);
211 for(
int k=0; k<=iter;++k ) {
212 z_->axpy((*
y_)(k),*(Z[k]));
215 if( (*res_)[iter+1] <= rtol ) {
226 return (*
res_)[iter];
229 return (*
res_)[iter+1];
243 #endif // ROL_GMRES_H
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Real getRelativeTolerance(void) const
Real getAbsoluteTolerance(void) const
virtual void plus(const Vector &x)=0
Compute , where .
GMRES(ROL::ParameterList &parlist)
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
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
ROL::Ptr< Vector< Real > > r_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
ROL::Ptr< std::ostream > outStream_
LA::Matrix< Real > SDMatrix
Real run(Vector< Real > &x, LinearOperator< Real > &A, const Vector< Real > &b, LinearOperator< Real > &M, int &iter, int &flag)
ROL::Ptr< std::vector< Real > > res_
Preconditioned GMRES solver.
virtual void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
Provides definitions for Krylov solvers.
Provides the interface to apply a linear operator.
ROL::Ptr< SDVector > cnorm_
void enableOutput(std::ostream &outStream)
ROL::Ptr< Vector< Real > > w_
LA::Vector< Real > SDVector
ROL::Ptr< Vector< Real > > z_
ROL::LAPACK< int, Real > lapack_