ROL
step/krylov/test_01.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
14 #include "ROL_StdVector.hpp"
15 #include "ROL_GMRES.hpp"
16 #include "ROL_KrylovFactory.hpp"
17 #include "ROL_RandomVector.hpp"
18 #include "ROL_ParameterList.hpp"
19 
20 #include "ROL_Stream.hpp"
21 #include "Teuchos_GlobalMPISession.hpp"
22 
23 #include<iomanip>
24 
25 // Identity operator for preconditioner
26 template<class Real>
27 class Identity : public ROL::LinearOperator<Real> {
29 public:
30  void apply( V& Hv, const V& v, Real &tol ) const {
31  Hv.set(v);
32  }
33 }; // class Identity
34 
35 
36 // Apply a tridiagonal Toeplitz matrix to a ROL::StdVector to test Krylov solvers
37 template<class Real>
39 
40  typedef std::vector<Real> vector;
43 
44  typedef typename vector::size_type uint;
45 
46 private:
47 
48  Real a_; // subdiagonal
49  Real b_; // diagonal
50  Real c_; // superdiagonal
51 
52  ROL::LAPACK<int,Real> lapack_;
53 
54 public:
55 
56  TridiagonalToeplitzOperator( Real &a, Real &b, Real &c ) : a_(a), b_(b), c_(c) {}
57 
58  // Tridiagonal multiplication
59  void apply( V &Hv, const V &v, Real &tol ) const {
60 
61 
62 
63  SV &Hvs = dynamic_cast<SV&>(Hv);
64  ROL::Ptr<vector> Hvp = Hvs.getVector();
65 
66  const SV &vs = dynamic_cast<const SV&>(v);
67  ROL::Ptr<const vector> vp = vs.getVector();
68 
69  uint n = vp->size();
70 
71  (*Hvp)[0] = b_*(*vp)[0] + c_*(*vp)[1];
72 
73  for(uint k=1; k<n-1; ++k) {
74  (*Hvp)[k] = a_*(*vp)[k-1] + b_*(*vp)[k] + c_*(*vp)[k+1];
75  }
76 
77  (*Hvp)[n-1] = a_*(*vp)[n-2] + b_*(*vp)[n-1];
78 
79  }
80 
81  // Tridiagonal solve - compare against GMRES
82  void applyInverse( V &Hv, const V &v, Real &tol ) const {
83 
84 
85 
86  SV &Hvs = dynamic_cast<SV&>(Hv);
87  ROL::Ptr<vector> Hvp = Hvs.getVector();
88 
89  const SV &vs = dynamic_cast<const SV&>(v);
90  ROL::Ptr<const vector> vp = vs.getVector();
91 
92  uint n = vp->size();
93 
94  const char TRANS = 'N';
95  const int NRHS = 1;
96 
97  vector dl(n-1,a_);
98  vector d(n,b_);
99  vector du(n-1,c_);
100  vector du2(n-2,0.0);
101 
102  std::vector<int> ipiv(n);
103  int info;
104 
105  Hv.set(v); // LAPACK will modify this in place
106 
107  // Do Tridiagonal LU factorization
108  lapack_.GTTRF(n,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
109 
110  // Solve the system with the LU factors
111  lapack_.GTTRS(TRANS,n,NRHS,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&(*Hvp)[0],n,&info);
112 
113  }
114 
115 }; // class TridiagonalToeplitzOperator
116 
117 
118 
119 typedef double RealT;
120 
121 int main(int argc, char *argv[]) {
122 
123 
124 
125 
126  typedef std::vector<RealT> vector;
127  typedef ROL::StdVector<RealT> SV;
128 
129  typedef typename vector::size_type uint;
130 
131  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
132 
133  int iprint = argc - 1;
134  ROL::Ptr<std::ostream> outStream;
135  ROL::nullstream bhs; // outputs nothing
136  if (iprint > 0)
137  outStream = ROL::makePtrFromRef(std::cout);
138  else
139  outStream = ROL::makePtrFromRef(bhs);
140 
141  int errorFlag = 0;
142 
143  try {
144 
145  ROL::ParameterList parlist;
146  ROL::ParameterList &gList = parlist.sublist("General");
147  ROL::ParameterList &kList = gList.sublist("Krylov");
148 
149  kList.set("Type","GMRES");
150  kList.set("Iteration Limit",20);
151  kList.set("Absolute Tolerance",1.e-8);
152  kList.set("Relative Tolerance",1.e-6);
153  kList.set("Use Initial Guess",false);
154 
155  uint dim = 10;
156 
157  ROL::Ptr<vector> xp = ROL::makePtr<vector>(dim,0.0);
158  ROL::Ptr<vector> yp = ROL::makePtr<vector>(dim,0.0);
159  ROL::Ptr<vector> zp = ROL::makePtr<vector>(dim,0.0);
160  ROL::Ptr<vector> bp = ROL::makePtr<vector>(dim,0.0);
161 
162  SV x(xp); // Exact solution
163  SV y(yp); // Solution using direct solve
164  SV z(zp); // Solution using GMRES
165 
166  SV b(bp); // Right-hand-side
167 
168  RealT left = -1.0;
169  RealT right = 1.0;
170 
171  ROL::RandomizeVector(x,left,right);
172 
173  RealT sub = -1.0;
174  RealT diag = 2.0;
175  RealT super = -1.0;
176 
177  TridiagonalToeplitzOperator<RealT> T(sub,diag,super);
178  Identity<RealT> I;
179 
180  RealT tol = 0.0;
181 
182  T.apply(b,x,tol);
183 
184  T.applyInverse(y,b,tol);
185 
186  ROL::Ptr<ROL::Krylov<RealT> > krylov = ROL::KrylovFactory<RealT>( parlist );
187 
188  int iter;
189  int flag;
190 
191  krylov->run(z,T,b,I,iter,flag);
192 
193  *outStream << std::setw(10) << "Exact"
194  << std::setw(10) << "LAPACK"
195  << std::setw(10) << "GMRES " << std::endl;
196  *outStream << "---------------------------------" << std::endl;
197 
198  for(uint k=0;k<dim;++k) {
199  *outStream << std::setw(10) << (*xp)[k] << " "
200  << std::setw(10) << (*yp)[k] << " "
201  << std::setw(10) << (*zp)[k] << " " << std::endl;
202  }
203 
204  *outStream << "GMRES performed " << iter << " iterations." << std::endl;
205 
206  z.axpy(-1.0,x);
207 
208  if( z.norm() > std::sqrt(ROL::ROL_EPSILON<RealT>()) ) {
209  ++errorFlag;
210  }
211 
212  }
213  catch (std::logic_error& err) {
214  *outStream << err.what() << "\n";
215  errorFlag = -1000;
216  }; // end try
217 
218  if (errorFlag != 0)
219  std::cout << "End Result: TEST FAILED\n";
220  else
221  std::cout << "End Result: TEST PASSED\n";
222 
223  return 0;
224 }
typename PV< Real >::size_type size_type
Ptr< const std::vector< Element > > getVector() const
void RandomizeVector(Vector< Real > &x, const Real &lower=0.0, const Real &upper=1.0)
Fill a ROL::Vector with uniformly-distributed random numbers in the interval [lower,upper].
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
TridiagonalToeplitzOperator(Real &a, Real &b, Real &c)
ROL::LAPACK< int, Real > lapack_
Provides the interface to apply a linear operator.
void applyInverse(V &Hv, const V &v, Real &tol) const
Apply inverse of linear operator.
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
int main(int argc, char *argv[])
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
constexpr auto dim
ROL::Vector< Real > V