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_;
94 GMRES( Teuchos::ParameterList &parlist ) : isInitialized_(false) {
100 Teuchos::ParameterList &gList = parlist.sublist(
"General");
101 Teuchos::ParameterList &kList = gList.sublist(
"Krylov");
103 useInexact_ = gList.get(
"Inexact Hessian-Times-A-Vector",
false);
104 maxit_ = kList.get(
"Iteration Limit",50);
105 absTol_ = kList.get(
"Absolute Tolerance", 1.e-4);
106 relTol_ = kList.get(
"Relative Tolerance", 1.e-2);
107 useInitialGuess_ = kList.get(
"Use Initial Guess",
false);
109 H_ = rcp(
new SDMatrix( maxit_+1, maxit_ ) );
110 cs_ = rcp(
new SDVector( maxit_ ) );
111 sn_ = rcp(
new SDVector( maxit_ ) );
112 s_ = rcp(
new SDVector( maxit_+1 ) );
113 y_ = rcp(
new SDVector( maxit_+1 ) );
114 cnorm_ = rcp(
new SDVector( maxit_ ) );
115 res_ = rcp(
new std::vector<Real>(maxit_+1,0.0) );
129 if ( !isInitialized_ ) {
134 isInitialized_ =
true;
140 if(useInitialGuess_) {
154 std::vector<RCP<Vector<Real > > > V;
155 std::vector<RCP<Vector<Real > > > Z;
157 (*res_)[0] = r_->norm();
159 Real rtol = std::min(absTol_,relTol_*(*res_)[0]);
161 V.push_back(b.
clone());
163 (V[0])->scale(one/(*res_)[0]);
165 (*s_)(0) = (*res_)[0];
167 for( iter=0; iter<
maxit_; ++iter ) {
172 itol = rtol/(maxit_*(*res_)[iter]);
175 Z.push_back(x.
clone());
178 M.
apply(*(Z[iter]),*(V[iter]),itol);
181 A.
apply(*w_,*(Z[iter]),itol);
184 for(
int k=0; k<=iter; ++k ) {
185 (*H_)(k,iter) = w_->dot(*(V[k]));
186 w_->axpy( -(*H_)(k,iter), *(V[k]) );
189 (*H_)(iter+1,iter) = w_->norm();
191 V.push_back( b.
clone() );
192 (V[iter+1])->
set(*w_);
193 (V[iter+1])->scale(one/((*H_)(iter+1,iter)));
196 for(
int k=0; k<=iter-1; ++k ) {
197 temp = (*cs_)(k)*(*H_)(k,iter) + (*sn_)(k)*(*H_)(k+1,iter);
198 (*H_)(k+1,iter) = -(*sn_)(k)*(*H_)(k,iter) + (*cs_)(k)*(*H_)(k+1,iter);
199 (*H_)(k,iter) = temp;
203 if( (*H_)(iter+1,iter) == zero ) {
207 else if ( std::abs((*H_)(iter+1,iter)) > std::abs((*H_)(iter,iter)) ) {
208 temp = (*H_)(iter,iter) / (*H_)(iter+1,iter);
209 (*sn_)(iter) = one / std::sqrt( one + temp*temp );
210 (*cs_)(iter) = temp*(*sn_)(iter);
213 temp = (*H_)(iter+1,iter) / (*H_)(iter,iter);
214 (*cs_)(iter) = one / std::sqrt( one + temp*temp );
215 (*sn_)(iter) = temp*(*cs_)(iter);
219 temp = (*cs_)(iter)*(*s_)(iter);
220 (*s_)(iter+1) = -(*sn_)(iter)*(*s_)(iter);
222 (*H_)(iter,iter) = (*cs_)(iter)*(*H_)(iter,iter) + (*sn_)(iter)*(*H_)(iter+1,iter);
223 (*H_)(iter+1,iter) = zero;
224 (*res_)[iter+1] = std::abs((*s_)(iter+1));
227 const char uplo =
'U';
228 const char trans =
'N';
229 const char diag =
'N';
230 const char normin =
'N';
234 lapack_.LATRS(uplo, trans, diag, normin, iter+1, H_->values(), maxit_+1, y_->values(), &scaling, cnorm_->values(), &info);
238 for(
int k=0; k<=iter;++k ) {
239 z_->axpy((*y_)(k),*(Z[k]));
242 if( (*res_)[iter+1] <= rtol ) {
260 #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_
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_
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Teuchos::RCP< Vector< Real > > z_