ROL
function/test_03.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 
52 #include "ROL_BlockOperator.hpp"
53 #include "ROL_DiagonalOperator.hpp"
54 #include "ROL_Elementwise_Function.hpp"
56 #include "ROL_StdVector.hpp"
57 #include "ROL_Types.hpp"
58 
59 #include "Teuchos_oblackholestream.hpp"
60 #include "Teuchos_GlobalMPISession.hpp"
61 
62 template<class Real>
63 void print_vector( const ROL::Vector<Real> &x ) {
64 
65  typedef ROL::Vector<Real> V;
66  typedef ROL::StdVector<Real> SV;
68  typedef typename PV::size_type size_type;
69 
70  const PV eb = Teuchos::dyn_cast<const PV>(x);
71  size_type n = eb.numVectors();
72 
73  for(size_type k=0; k<n; ++k) {
74  std::cout << "[subvector " << k << "]" << std::endl;
75  Teuchos::RCP<const V> vec = eb.get(k);
76  Teuchos::RCP<const std::vector<Real> > vp =
77  Teuchos::dyn_cast<SV>(const_cast<V&>(*vec)).getVector();
78  for(size_type i=0;i<vp->size();++i) {
79  std::cout << (*vp)[i] << std::endl;
80  }
81  }
82 }
83 
84 // Implemantation of a Dyadic operator x*y'
85 template<class Real>
86 class DyadicOperator : public ROL::LinearOperator<Real> {
87 
89 
90 private:
91 
92  const Teuchos::RCP<const V> x_;
93  const Teuchos::RCP<const V> y_;
94 
95 public:
96 
97  DyadicOperator( const Teuchos::RCP<const V> &x,
98  const Teuchos::RCP<const V> &y ) : x_(x), y_(y) {}
99 
100  void apply( V &Hv, const V &v, Real &tol ) const {
101  Hv.set(*x_);
102  Hv.scale(v.dot(*y_));
103  }
104 
105 };
106 
107 // Implementation of a Null operator
108 template<class Real>
109 class NullOperator : public ROL::LinearOperator<Real> {
110 
112 
113 public:
114 
115  void apply( V &Hv, const V &v, Real &tol ) const {
116  Hv.zero();
117  }
118 };
119 
120 
121 
122 
123 typedef double RealT;
124 
125 int main(int argc, char *argv[]) {
126 
127  // Define vectors
128  typedef std::vector<RealT> vector;
129  typedef ROL::Vector<RealT> V;
130  typedef ROL::StdVector<RealT> SV;
131 
132  // Define operators
133  typedef ROL::LinearOperator<RealT> LinOp;
134  typedef ROL::DiagonalOperator<RealT> DiagOp;
135  typedef DyadicOperator<RealT> DyadOp;
136  typedef NullOperator<RealT> NullOp;
137 
138  typedef typename vector::size_type uint;
139 
140  using namespace Teuchos;
141 
142  GlobalMPISession mpiSession(&argc, &argv);
143 
144  int iprint = argc - 1;
145 
146  RCP<std::ostream> outStream;
147  oblackholestream bhs; // no output
148 
149  if( iprint>0 )
150  outStream = rcp(&std::cout,false);
151  else
152  outStream = rcp(&bhs,false);
153 
154  int errorFlag = 0;
155 
156  double errtol = ROL::ROL_THRESHOLD;
157 
158  try {
159 
160  uint dim = 3; // Number of elements in each subvector (could be different)
161 
162  RCP<vector> x1_rcp = rcp( new vector(dim,1.0) );
163  RCP<vector> x2_rcp = rcp( new vector(dim,2.0) );
164 
165  RCP<vector> y1_rcp = rcp( new vector(dim,0.0) );
166  RCP<vector> y2_rcp = rcp( new vector(dim,0.0) );
167 
168  RCP<vector> z1_rcp = rcp( new vector(dim,0.0) );
169  RCP<vector> z2_rcp = rcp( new vector(dim,0.0) );
170 
171  RCP<V> x1 = rcp( new SV( x1_rcp) );
172  RCP<V> x2 = rcp( new SV( x2_rcp) );
173 
174  RCP<V> y1 = rcp( new SV( y1_rcp) );
175  RCP<V> y2 = rcp( new SV( y2_rcp) );
176 
177  RCP<V> z1 = rcp( new SV( z1_rcp) );
178  RCP<V> z2 = rcp( new SV( z2_rcp) );
179 
180  RCP<V> x = ROL::CreatePartitionedVector( x1, x2 );
181  RCP<V> y = ROL::CreatePartitionedVector( y1, y2 );
182  RCP<V> z = ROL::CreatePartitionedVector( z1, z2 );
183 
184  // Operator diagonals
185  RCP<vector> d1_rcp = rcp( new vector(dim,0.0) );
186  RCP<vector> d2_rcp = rcp( new vector(dim,0.0) );
187 
188  // Dyadic components
189  RCP<vector> u_rcp = rcp( new vector(dim,0.0) );
190  RCP<vector> v_rcp = rcp( new vector(dim,1.0) );
191 
192  (*d1_rcp)[0] = 6.0; (*d2_rcp)[0] = 3.0;
193  (*d1_rcp)[1] = 5.0; (*d2_rcp)[1] = 2.0;
194  (*d1_rcp)[2] = 4.0; (*d2_rcp)[2] = 1.0;
195 
196  (*z1_rcp)[0] = 6.0; (*z2_rcp)[0] = 6.0;
197  (*z1_rcp)[1] = 11.0; (*z2_rcp)[1] = 4.0;
198  (*z1_rcp)[2] = 4.0; (*z2_rcp)[2] = 2.0;
199 
200  (*u_rcp)[1] = 1.0;
201 
202  RCP<V> d1 = rcp( new SV(d1_rcp) );
203  RCP<V> d2 = rcp( new SV(d2_rcp) );
204  RCP<V> u = rcp( new SV(u_rcp) );
205  RCP<V> v = rcp( new SV(v_rcp) );
206 
207  RCP<LinOp> D1 = rcp( new DiagOp(*d1) );
208  RCP<LinOp> NO = rcp( new NullOp() );
209  RCP<LinOp> UV = rcp( new DyadOp(u,v) );
210  RCP<LinOp> D2 = rcp( new DiagOp(*d2) );
211 
212 
213  RealT tol = 0.0;
214 
215  D1->apply(*x1,*x1,tol);
216  D1->applyInverse(*x1,*x1,tol);
217 
218  ROL::BlockOperator2<RealT> bkop(D1,UV,NO,D2);
219 
220 
221  bkop.apply(*y,*x,tol);
222 
223  z->axpy(-1.0,*y);
224 
225  errorFlag += static_cast<int>(z->norm()>errtol);
226 
227  }
228 
229  catch (std::logic_error err) {
230 
231 
232  }; // end try
233 
234  if (errorFlag != 0)
235  std::cout << "End Result: TEST FAILED\n";
236  else
237  std::cout << "End Result: TEST PASSED\n";
238 
239  return 0;
240 }
241 
int main(int argc, char *argv[])
virtual void scale(const Real alpha)=0
Compute where .
Defines the linear algebra of vector space on a generic partitioned vector.
double RealT
Contains definitions of custom data types in ROL.
static const double ROL_THRESHOLD
Tolerance for various equality tests.
Definition: ROL_Types.hpp:122
Teuchos::RCP< Vector< Real > > CreatePartitionedVector(const Teuchos::RCP< Vector< Real > > &a)
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void print_vector(const ROL::Vector< Real > &x)
virtual Real dot(const Vector &x) const =0
Compute where .
Provides the interface to apply a diagonal operator which acts like elementwise multiplication when a...
ROL::Vector< Real > V
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
ROL::Vector< Real > V
Provides the interface to apply a linear operator.
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
const Teuchos::RCP< const V > y_
double RealT
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
DyadicOperator(const Teuchos::RCP< const V > &x, const Teuchos::RCP< const V > &y)
const Teuchos::RCP< const V > x_