ROL
step/krylov/test_01.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
48 #include "ROL_StdVector.hpp"
49 #include "ROL_GMRES.hpp"
50 #include "ROL_KrylovFactory.hpp"
51 #include "ROL_RandomVector.hpp"
52 
53 #include "Teuchos_oblackholestream.hpp"
54 #include "Teuchos_GlobalMPISession.hpp"
55 
56 #include<iomanip>
57 
58 // Identity operator for preconditioner
59 template<class Real>
60 class Identity : public ROL::LinearOperator<Real> {
62 public:
63  void apply( V& Hv, const V& v, Real &tol ) const {
64  Hv.set(v);
65  }
66 }; // class Identity
67 
68 
69 // Apply a tridiagonal Toeplitz matrix to a ROL::StdVector to test Krylov solvers
70 template<class Real>
72 
73  typedef std::vector<Real> vector;
76 
77  typedef typename vector::size_type uint;
78 
79 private:
80 
81  Real a_; // subdiagonal
82  Real b_; // diagonal
83  Real c_; // superdiagonal
84 
85  Teuchos::LAPACK<int,Real> lapack_;
86 
87 public:
88 
89  TridiagonalToeplitzOperator( Real &a, Real &b, Real &c ) : a_(a), b_(b), c_(c) {}
90 
91  // Tridiagonal multiplication
92  void apply( V &Hv, const V &v, Real &tol ) const {
93 
94  using Teuchos::RCP; using Teuchos::dyn_cast;
95 
96  SV &Hvs = dyn_cast<SV>(Hv);
97  RCP<vector> Hvp = Hvs.getVector();
98 
99  const SV &vs = dyn_cast<const SV>(v);
100  RCP<const vector> vp = vs.getVector();
101 
102  uint n = vp->size();
103 
104  (*Hvp)[0] = b_*(*vp)[0] + c_*(*vp)[1];
105 
106  for(uint k=1; k<n-1; ++k) {
107  (*Hvp)[k] = a_*(*vp)[k-1] + b_*(*vp)[k] + c_*(*vp)[k+1];
108  }
109 
110  (*Hvp)[n-1] = a_*(*vp)[n-2] + b_*(*vp)[n-1];
111 
112  }
113 
114  // Tridiagonal solve - compare against GMRES
115  void applyInverse( V &Hv, const V &v, Real &tol ) const {
116 
117  using Teuchos::RCP; using Teuchos::dyn_cast;
118 
119  SV &Hvs = dyn_cast<SV>(Hv);
120  RCP<vector> Hvp = Hvs.getVector();
121 
122  const SV &vs = dyn_cast<const SV>(v);
123  RCP<const vector> vp = vs.getVector();
124 
125  uint n = vp->size();
126 
127  const char TRANS = 'N';
128  const int NRHS = 1;
129 
130  vector dl(n-1,a_);
131  vector d(n,b_);
132  vector du(n-1,c_);
133  vector du2(n-2,0.0);
134 
135  std::vector<int> ipiv(n);
136  int info;
137 
138  Hv.set(v); // LAPACK will modify this in place
139 
140  // Do Tridiagonal LU factorization
141  lapack_.GTTRF(n,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
142 
143  // Solve the system with the LU factors
144  lapack_.GTTRS(TRANS,n,NRHS,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&(*Hvp)[0],n,&info);
145 
146  }
147 
148 }; // class TridiagonalToeplitzOperator
149 
150 
151 
152 typedef double RealT;
153 
154 int main(int argc, char *argv[]) {
155 
156  using Teuchos::RCP;
157  using Teuchos::rcp;
158 
159  typedef std::vector<RealT> vector;
160  typedef ROL::StdVector<RealT> SV;
161 
162  typedef typename vector::size_type uint;
163 
164  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
165 
166  int iprint = argc - 1;
167  RCP<std::ostream> outStream;
168  Teuchos::oblackholestream bhs; // outputs nothing
169  if (iprint > 0)
170  outStream = rcp(&std::cout, false);
171  else
172  outStream = rcp(&bhs, false);
173 
174  int errorFlag = 0;
175 
176  try {
177 
178  Teuchos::ParameterList parlist;
179  Teuchos::ParameterList &gList = parlist.sublist("General");
180  Teuchos::ParameterList &kList = gList.sublist("Krylov");
181 
182  kList.set("Type","GMRES");
183  kList.set("Iteration Limit",20);
184  kList.set("Absolute Tolerance",1.e-8);
185  kList.set("Relative Tolerance",1.e-6);
186  kList.set("Use Initial Guess",false);
187 
188  uint dim = 10;
189 
190  RCP<vector> xp = rcp( new vector(dim,0.0) );
191  RCP<vector> yp = rcp( new vector(dim,0.0) );
192  RCP<vector> zp = rcp( new vector(dim,0.0) );
193  RCP<vector> bp = rcp( new vector(dim,0.0) );
194 
195  SV x(xp); // Exact solution
196  SV y(yp); // Solution using direct solve
197  SV z(zp); // Solution using GMRES
198 
199  SV b(bp); // Right-hand-side
200 
201  RealT left = -1.0;
202  RealT right = 1.0;
203 
204  ROL::RandomizeVector(x,left,right);
205 
206  RealT sub = -1.0;
207  RealT diag = 2.0;
208  RealT super = -1.0;
209 
210  TridiagonalToeplitzOperator<RealT> T(sub,diag,super);
211  Identity<RealT> I;
212 
213  RealT tol = 0.0;
214 
215  T.apply(b,x,tol);
216 
217  T.applyInverse(y,b,tol);
218 
219  RCP<ROL::Krylov<RealT> > krylov = ROL::KrylovFactory<RealT>( parlist );
220 
221  int iter;
222  int flag;
223 
224  krylov->run(z,T,b,I,iter,flag);
225 
226  *outStream << std::setw(10) << "Exact"
227  << std::setw(10) << "LAPACK"
228  << std::setw(10) << "GMRES " << std::endl;
229  *outStream << "---------------------------------" << std::endl;
230 
231  for(uint k=0;k<dim;++k) {
232  *outStream << std::setw(10) << (*xp)[k] << " "
233  << std::setw(10) << (*yp)[k] << " "
234  << std::setw(10) << (*zp)[k] << " " << std::endl;
235  }
236 
237  *outStream << "GMRES performed " << iter << " iterations." << std::endl;
238 
239  z.axpy(-1.0,x);
240 
241  if( z.norm() > std::sqrt(ROL::ROL_EPSILON<RealT>()) ) {
242  ++errorFlag;
243  }
244 
245  }
246  catch (std::logic_error err) {
247  *outStream << err.what() << "\n";
248  errorFlag = -1000;
249  }; // end try
250 
251  if (errorFlag != 0)
252  std::cout << "End Result: TEST FAILED\n";
253  else
254  std::cout << "End Result: TEST PASSED\n";
255 
256  return 0;
257 }
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
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].
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void applyInverse(V &Hv, const V &v, Real &tol) const
Apply inverse of linear operator.
TridiagonalToeplitzOperator(Real &a, Real &b, Real &c)
Teuchos::RCP< const std::vector< Element > > getVector() const
int main(int argc, char *argv[])
Provides the interface to apply a linear operator.
Teuchos::LAPACK< int, Real > lapack_
double RealT
double RealT
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
ROL::Vector< Real > V
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.