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