56 #include "Teuchos_SerialDenseMatrix.hpp" 57 #include "Teuchos_SerialDenseVector.hpp" 58 #include "Teuchos_LAPACK.hpp" 65 typedef Teuchos::SerialDenseMatrix<int, Real>
SDMatrix;
66 typedef Teuchos::SerialDenseVector<int, Real>
SDVector;
70 Teuchos::RCP<Vector<Real> >
r_;
71 Teuchos::RCP<Vector<Real> >
z_;
72 Teuchos::RCP<Vector<Real> >
w_;
74 Teuchos::RCP<SDMatrix>
H_;
75 Teuchos::RCP<SDVector>
cs_;
76 Teuchos::RCP<SDVector>
sn_;
77 Teuchos::RCP<SDVector>
s_;
78 Teuchos::RCP<SDVector>
y_;
81 Teuchos::RCP<std::vector<Real> >
res_;
100 Real zero(0), oem2(1.e-2), oem4(1.e-4);
102 Teuchos::ParameterList &gList = parlist.sublist(
"General");
103 Teuchos::ParameterList &kList = gList.sublist(
"Krylov");
105 useInexact_ = gList.get(
"Inexact Hessian-Times-A-Vector",
false);
106 maxit_ = kList.get(
"Iteration Limit",50);
107 absTol_ = kList.get(
"Absolute Tolerance", oem4);
108 relTol_ = kList.get(
"Relative Tolerance", oem2);
117 res_ = rcp(
new std::vector<Real>(
maxit_+1,zero) );
128 Real zero(0), one(1);
138 Real itol = std::sqrt(ROL_EPSILON<Real>());
155 std::vector<RCP<Vector<Real > > > V;
156 std::vector<RCP<Vector<Real > > > Z;
158 (*res_)[0] =
r_->norm();
162 V.push_back(b.
clone());
164 (V[0])->scale(one/(*
res_)[0]);
166 (*s_)(0) = (*
res_)[0];
168 for( iter=0; iter<
maxit_; ++iter ) {
173 itol = rtol/(
maxit_*(*res_)[iter]);
176 Z.push_back(x.
clone());
185 for(
int k=0; k<=iter; ++k ) {
186 (*H_)(k,iter) =
w_->dot(*(V[k]));
187 w_->axpy( -(*
H_)(k,iter), *(V[k]) );
190 (*H_)(iter+1,iter) =
w_->norm();
192 V.push_back( b.
clone() );
193 (V[iter+1])->
set(*
w_);
194 (V[iter+1])->scale(one/((*
H_)(iter+1,iter)));
197 for(
int k=0; k<=iter-1; ++k ) {
198 temp = (*cs_)(k)*(*
H_)(k,iter) + (*
sn_)(k)*(*
H_)(k+1,iter);
199 (*H_)(k+1,iter) = -(*
sn_)(k)*(*
H_)(k,iter) + (*
cs_)(k)*(*
H_)(k+1,iter);
200 (*H_)(k,iter) = temp;
204 if( (*
H_)(iter+1,iter) == zero ) {
208 else if ( std::abs((*
H_)(iter+1,iter)) > std::abs((*
H_)(iter,iter)) ) {
209 temp = (*H_)(iter,iter) / (*
H_)(iter+1,iter);
210 (*sn_)(iter) = one / std::sqrt( one + temp*temp );
211 (*cs_)(iter) = temp*(*
sn_)(iter);
214 temp = (*H_)(iter+1,iter) / (*
H_)(iter,iter);
215 (*cs_)(iter) = one / std::sqrt( one + temp*temp );
216 (*sn_)(iter) = temp*(*
cs_)(iter);
220 temp = (*cs_)(iter)*(*
s_)(iter);
221 (*s_)(iter+1) = -(*
sn_)(iter)*(*
s_)(iter);
223 (*H_)(iter,iter) = (*
cs_)(iter)*(*
H_)(iter,iter) + (*
sn_)(iter)*(*
H_)(iter+1,iter);
224 (*H_)(iter+1,iter) = zero;
225 (*res_)[iter+1] = std::abs((*
s_)(iter+1));
228 const char uplo =
'U';
229 const char trans =
'N';
230 const char diag =
'N';
231 const char normin =
'N';
235 lapack_.LATRS(uplo, trans, diag, normin, iter+1,
H_->values(),
maxit_+1,
y_->values(), &scaling,
cnorm_->values(), &info);
239 for(
int k=0; k<=iter;++k ) {
240 z_->axpy((*
y_)(k),*(Z[k]));
243 if( (*
res_)[iter+1] <= rtol ) {
261 #endif // ROL_GMRES_H
virtual void plus(const Vector &x)=0
Compute , where .
Teuchos::RCP< SDVector > cs_
Contains definitions of custom data types in ROL.
GMRES(Teuchos::ParameterList &parlist)
Teuchos::RCP< Vector< Real > > r_
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const =0
Apply linear operator.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
Preconditioned GMRES solver.
Teuchos::RCP< SDVector > y_
Teuchos::RCP< SDVector > cnorm_
Teuchos::RCP< SDVector > sn_
Teuchos::RCP< SDMatrix > H_
virtual void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
void run(Vector< Real > &x, LinearOperator< Real > &A, const Vector< Real > &b, LinearOperator< Real > &M, int &iter, int &flag)
Provides definitions for Krylov solvers.
Teuchos::SerialDenseMatrix< int, Real > SDMatrix
Provides the interface to apply a linear operator.
Teuchos::RCP< Vector< Real > > w_
Teuchos::SerialDenseVector< int, Real > SDVector
Teuchos::LAPACK< int, Real > lapack_
Teuchos::RCP< std::vector< Real > > res_
Teuchos::RCP< SDVector > s_
Teuchos::RCP< Vector< Real > > z_