ROL
function/operator/test_02.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 
15 #include "ROL_Stream.hpp"
16 #include "Teuchos_GlobalMPISession.hpp"
17 
18 typedef double RealT;
19 
20 int main(int argc, char *argv[]) {
21 
22 
23 
24  using SV = ROL::StdVector<RealT>;
27 
28  using vector = std::vector<RealT>;
29 
30 
31  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
32 
33  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
34  int iprint = argc - 1;
35  ROL::Ptr<std::ostream> outStream;
36  ROL::nullstream bhs; // outputs nothing
37  if (iprint > 0)
38  outStream = ROL::makePtrFromRef(std::cout);
39  else
40  outStream = ROL::makePtrFromRef(bhs);
41 
42  // Save the format state of the original std::cout.
43  ROL::nullstream oldFormatState;
44  oldFormatState.copyfmt(std::cout);
45 
46  int errorFlag = 0;
47 
48  RealT tol = ROL::ROL_EPSILON<RealT>();
49 
50  // *** Test body.
51 
52  try {
53 
54  int dim = 3;
55 
56  ROL::Ptr<vector> m_ptr = ROL::makePtr<vector>(
57  std::initializer_list<RealT>{ 3.0, 1.0, 0.0, -2.0, 6.0, 2.0, 0.0, -1.0, 3.0 } );
58  ROL::Ptr<vector> a_ptr = ROL::makePtr<vector>(
59  std::initializer_list<RealT>{ 3.0, 6.0, 3.0} );
60  ROL::Ptr<vector> b_ptr = ROL::makePtr<vector>(
61  std::initializer_list<RealT>{ -2.0, -1.0 } );
62  ROL::Ptr<vector> c_ptr = ROL::makePtr<vector>(
63  std::initializer_list<RealT>{ 1.0, 2.0 } );
64 
65  MAT M(m_ptr);
66  TRI T(a_ptr,b_ptr,c_ptr);
67 
68  SV xm( ROL::makePtr<vector>( std::initializer_list<RealT>{1.0, 2.0, -1.0} ) );
69  SV ym( ROL::makePtr<vector>( dim ) );
70  SV zm( ROL::makePtr<vector>( dim ) );
71 
72  SV xt( ROL::makePtr<vector>( dim ) );
73  SV yt( ROL::makePtr<vector>( dim ) );
74  SV zt( ROL::makePtr<vector>( dim ) );
75 
76  SV error( ROL::makePtr<vector>(dim) );
77  RealT nerr = 0;
78 
79  xt.set(xm);
80 
81  M.apply(ym,xm,tol);
82  M.applyInverse(zm,ym,tol);
83 
84  *outStream << "\nUsing StdLinearOperator - A is full matrix representation" << std::endl;
85  *outStream << "x = "; xm.print(*outStream);
86  *outStream << "y = Ax = "; ym.print(*outStream);
87  *outStream << "z = inv(A)y = "; zm.print(*outStream);
88 
89  *outStream << "\nUsing StdTridiagonalOperator - T is tridiagonal representation" << std::endl;
90  T.apply(yt,xt,tol);
91  *outStream << "y = Tx = "; yt.print(*outStream);
92  error.set(yt);
93  error.axpy(-1.0,ym);
94  nerr = error.norm();
95  errorFlag += static_cast<int>(nerr>tol);
96  *outStream << "apply() error = " << nerr <<std::endl;
97 
98  T.applyInverse(zt,yt,tol);
99  *outStream << "z = inv(T)y = "; yt.print(*outStream);
100  error.set(zt);
101  error.axpy(-1.0,zm);
102  nerr = error.norm();
103  errorFlag += static_cast<int>(nerr>tol);
104  *outStream << "applyInverse() error = " << nerr <<std::endl;
105 
106  M.applyAdjoint(ym,xm,tol);
107  M.applyAdjointInverse(zm,ym,tol);
108  *outStream << "\nUsing StdLinearOperator - A is full matrix representation" << std::endl;
109  *outStream << "x = "; xm.print(*outStream);
110  *outStream << "y = A'x = "; ym.print(*outStream);
111  *outStream << "z = inv(A')y = "; zm.print(*outStream);
112 
113  *outStream << "\nUsing StdTridiagonalOperator - T is tridiagonal representation" << std::endl;
114  T.applyAdjoint(yt,xt,tol);
115  *outStream << "y = T'x = "; yt.print(*outStream);
116  error.set(yt);
117  error.axpy(-1.0,ym);
118  nerr = error.norm();
119  errorFlag += static_cast<int>(nerr>tol);
120  *outStream << "applyAdjoint() error = " << nerr <<std::endl;
121 
122  T.applyAdjointInverse(zt,yt,tol);
123  *outStream << "z = inv(T')y = "; yt.print(*outStream);
124  error.set(zt);
125  error.axpy(-1.0,zm);
126  nerr = error.norm();
127  errorFlag += static_cast<int>(nerr>tol);
128  *outStream << "applyAdjointInverse() error = " << nerr <<std::endl;
129 
130 
131 
132 
133  }
134  catch (std::logic_error& err) {
135  *outStream << err.what() << "\n";
136  errorFlag = -1000;
137  }; // end try
138 
139  if (errorFlag != 0)
140  std::cout << "End Result: TEST FAILED\n";
141  else
142  std::cout << "End Result: TEST PASSED\n";
143 
144  // reset format state of std::cout
145  std::cout.copyfmt(oldFormatState);
146 
147  return 0;
148 
149 }
150 
Provides the std::vector implementation to apply a linear operator, which is a std::vector representa...
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Provides the std::vector implementation to apply a linear operator, which encapsulates a tridiagonal ...
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
int main(int argc, char *argv[])
constexpr auto dim