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 
44 /*! \file  test_10.cpp
45     \brief Show how to use CompositeConstraint interface.
46 */
47 
48 #include "ROL_CompositeConstraint_SimOpt.hpp"
49 #include "ROL_StdVector.hpp"
50 #include "ROL_Stream.hpp"
51 #include "Teuchos_GlobalMPISession.hpp"
52 
53 #include <iostream>
54 
55 template<class Real>
56 class valConstraint : public ROL::Constraint_SimOpt<Real> {
57 public:
valConstraint(void)58   valConstraint(void) : ROL::Constraint_SimOpt<Real>() {}
59 
value(ROL::Vector<Real> & c,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)60   void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
61     ROL::Ptr<std::vector<Real> > cp
62       = dynamic_cast<ROL::StdVector<Real>&>(c).getVector();
63     ROL::Ptr<const std::vector<Real> > up
64       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
65     ROL::Ptr<const std::vector<Real> > zp
66       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
67 
68     Real half(0.5), two(2);
69     // C(0) = U(0) - Z(0)
70     (*cp)[0] = (*up)[0]-(*zp)[0];
71     // C(1) = 0.5 * (U(0) + U(1) - Z(0))^2
72     (*cp)[1] = half*std::pow((*up)[0]+(*up)[1]-(*zp)[0],two);
73   }
74 
applyJacobian_1(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)75   void applyJacobian_1(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v,
76                        const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
77     ROL::Ptr<std::vector<Real> > jvp
78       = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
79     ROL::Ptr<const std::vector<Real> > vp
80       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
81     ROL::Ptr<const std::vector<Real> > up
82       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
83     ROL::Ptr<const std::vector<Real> > zp
84       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
85     (*jvp)[0] = (*vp)[0];
86     (*jvp)[1] = ((*up)[0] + (*up)[1] - (*zp)[0]) * ((*vp)[0] + (*vp)[1]);
87   }
88 
applyJacobian_2(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)89   void applyJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v,
90                        const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
91     ROL::Ptr<std::vector<Real> > jvp
92       = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
93     ROL::Ptr<const std::vector<Real> > vp
94       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
95     ROL::Ptr<const std::vector<Real> > up
96       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
97     ROL::Ptr<const std::vector<Real> > zp
98       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
99     (*jvp)[0] = -(*vp)[0];
100     (*jvp)[1] = ((*zp)[0] - (*up)[0] - (*up)[1]) * (*vp)[0];
101   }
102 
applyAdjointJacobian_1(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)103   void applyAdjointJacobian_1(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v,
104                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
105     ROL::Ptr<std::vector<Real> > ajvp
106       = dynamic_cast<ROL::StdVector<Real>&>(ajv).getVector();
107     ROL::Ptr<const std::vector<Real> > vp
108       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
109     ROL::Ptr<const std::vector<Real> > up
110       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
111     ROL::Ptr<const std::vector<Real> > zp
112       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
113     (*ajvp)[0] = (*vp)[0] + ((*up)[0] + (*up)[1] - (*zp)[0]) * (*vp)[1];
114     (*ajvp)[1] = ((*up)[0] + (*up)[1] - (*zp)[0]) * (*vp)[1];
115   }
116 
applyAdjointJacobian_2(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)117   void applyAdjointJacobian_2(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v,
118                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
119     ROL::Ptr<std::vector<Real> > ajvp
120       = dynamic_cast<ROL::StdVector<Real>&>(ajv).getVector();
121     ROL::Ptr<const std::vector<Real> > vp
122       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
123     ROL::Ptr<const std::vector<Real> > up
124       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
125     ROL::Ptr<const std::vector<Real> > zp
126       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
127     (*ajvp)[0] = ((*zp)[0] - (*up)[0] - (*up)[1]) * (*vp)[1] - (*vp)[0];
128   }
129 
applyAdjointHessian_11(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)130   void applyAdjointHessian_11(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
131                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
132     ROL::Ptr<std::vector<Real> > ahwvp
133       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
134     ROL::Ptr<const std::vector<Real> > wp
135       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
136     ROL::Ptr<const std::vector<Real> > vp
137       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
138     ROL::Ptr<const std::vector<Real> > up
139       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
140     ROL::Ptr<const std::vector<Real> > zp
141       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
142     (*ahwvp)[0] = (*wp)[1] * ((*vp)[0] + (*vp)[1]);
143     (*ahwvp)[1] = (*wp)[1] * ((*vp)[0] + (*vp)[1]);
144   }
145 
applyAdjointHessian_12(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)146   void applyAdjointHessian_12(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
147                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
148     ROL::Ptr<std::vector<Real> > ahwvp
149       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
150     ROL::Ptr<const std::vector<Real> > wp
151       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
152     ROL::Ptr<const std::vector<Real> > vp
153       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
154     ROL::Ptr<const std::vector<Real> > up
155       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
156     ROL::Ptr<const std::vector<Real> > zp
157       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
158     (*ahwvp)[0] = -(*wp)[1] * ((*vp)[0] + (*vp)[1]);
159   }
160 
applyAdjointHessian_21(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)161   void applyAdjointHessian_21(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
162                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
163     ROL::Ptr<std::vector<Real> > ahwvp
164       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
165     ROL::Ptr<const std::vector<Real> > wp
166       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
167     ROL::Ptr<const std::vector<Real> > vp
168       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
169     ROL::Ptr<const std::vector<Real> > up
170       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
171     ROL::Ptr<const std::vector<Real> > zp
172       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
173     (*ahwvp)[0] = -(*wp)[1] * (*vp)[0];
174     (*ahwvp)[1] = -(*wp)[1] * (*vp)[0];
175   }
176 
applyAdjointHessian_22(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)177   void applyAdjointHessian_22(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
178                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
179     ROL::Ptr<std::vector<Real> > ahwvp
180       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
181     ROL::Ptr<const std::vector<Real> > wp
182       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
183     ROL::Ptr<const std::vector<Real> > vp
184       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
185     ROL::Ptr<const std::vector<Real> > up
186       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
187     ROL::Ptr<const std::vector<Real> > zp
188       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
189     (*ahwvp)[0] = (*wp)[1] * (*vp)[0];
190   }
191 };
192 
193 template<class Real>
194 class redConstraint : public ROL::Constraint_SimOpt<Real> {
195 public:
redConstraint(void)196   redConstraint(void) : ROL::Constraint_SimOpt<Real>() {}
197 
value(ROL::Vector<Real> & c,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)198   void value(ROL::Vector<Real> &c, const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
199     ROL::Ptr<std::vector<Real> > cp
200       = dynamic_cast<ROL::StdVector<Real>&>(c).getVector();
201     ROL::Ptr<const std::vector<Real> > up
202       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
203     ROL::Ptr<const std::vector<Real> > zp
204       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
205 
206     const Real one(1), two(2);
207     // C = exp(U) - (Z^2 + 1)
208     (*cp)[0] = std::exp((*up)[0])-(std::pow((*zp)[0],two) + one);
209   }
210 
applyJacobian_1(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)211   void applyJacobian_1(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v,
212                        const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
213     ROL::Ptr<std::vector<Real> > jvp
214       = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
215     ROL::Ptr<const std::vector<Real> > vp
216       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
217     ROL::Ptr<const std::vector<Real> > up
218       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
219     ROL::Ptr<const std::vector<Real> > zp
220       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
221     (*jvp)[0] = std::exp((*up)[0]) * (*vp)[0];
222   }
223 
applyJacobian_2(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)224   void applyJacobian_2(ROL::Vector<Real> &jv, const ROL::Vector<Real> &v,
225                        const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
226     ROL::Ptr<std::vector<Real> > jvp
227       = dynamic_cast<ROL::StdVector<Real>&>(jv).getVector();
228     ROL::Ptr<const std::vector<Real> > vp
229       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
230     ROL::Ptr<const std::vector<Real> > up
231       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
232     ROL::Ptr<const std::vector<Real> > zp
233       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
234 
235     const Real two(2);
236     (*jvp)[0] = -two * (*zp)[0] * (*vp)[0];
237   }
238 
applyAdjointJacobian_1(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)239   void applyAdjointJacobian_1(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v,
240                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
241     ROL::Ptr<std::vector<Real> > ajvp
242       = dynamic_cast<ROL::StdVector<Real>&>(ajv).getVector();
243     ROL::Ptr<const std::vector<Real> > vp
244       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
245     ROL::Ptr<const std::vector<Real> > up
246       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
247     ROL::Ptr<const std::vector<Real> > zp
248       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
249     (*ajvp)[0] = std::exp((*up)[0]) * (*vp)[0];
250   }
251 
applyAdjointJacobian_2(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)252   void applyAdjointJacobian_2(ROL::Vector<Real> &ajv, const ROL::Vector<Real> &v,
253                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
254     ROL::Ptr<std::vector<Real> > ajvp
255       = dynamic_cast<ROL::StdVector<Real>&>(ajv).getVector();
256     ROL::Ptr<const std::vector<Real> > vp
257       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
258     ROL::Ptr<const std::vector<Real> > up
259       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
260     ROL::Ptr<const std::vector<Real> > zp
261       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
262 
263     const Real two(2);
264     (*ajvp)[0] = -two * (*zp)[0] * (*vp)[0];
265   }
266 
applyInverseJacobian_1(ROL::Vector<Real> & ijv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)267   void applyInverseJacobian_1(ROL::Vector<Real> &ijv, const ROL::Vector<Real> &v,
268                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
269     ROL::Ptr<std::vector<Real> > ijvp
270       = dynamic_cast<ROL::StdVector<Real>&>(ijv).getVector();
271     ROL::Ptr<const std::vector<Real> > vp
272       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
273     ROL::Ptr<const std::vector<Real> > up
274       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
275     ROL::Ptr<const std::vector<Real> > zp
276       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
277     (*ijvp)[0] = (*vp)[0] / std::exp((*up)[0]);
278   }
279 
applyInverseAdjointJacobian_1(ROL::Vector<Real> & ijv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)280   void applyInverseAdjointJacobian_1(ROL::Vector<Real> &ijv, const ROL::Vector<Real> &v,
281                                      const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
282     ROL::Ptr<std::vector<Real> > ijvp
283       = dynamic_cast<ROL::StdVector<Real>&>(ijv).getVector();
284     ROL::Ptr<const std::vector<Real> > vp
285       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
286     ROL::Ptr<const std::vector<Real> > up
287       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
288     ROL::Ptr<const std::vector<Real> > zp
289       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
290     (*ijvp)[0] = (*vp)[0] / std::exp((*up)[0]);
291   }
292 
applyAdjointHessian_11(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)293   void applyAdjointHessian_11(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
294                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
295     ROL::Ptr<std::vector<Real> > ahwvp
296       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
297     ROL::Ptr<const std::vector<Real> > wp
298       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
299     ROL::Ptr<const std::vector<Real> > vp
300       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
301     ROL::Ptr<const std::vector<Real> > up
302       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
303     ROL::Ptr<const std::vector<Real> > zp
304       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
305     (*ahwvp)[0] = std::exp((*up)[0]) * (*wp)[0] * (*vp)[0];
306   }
307 
applyAdjointHessian_12(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)308   void applyAdjointHessian_12(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
309                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
310     ROL::Ptr<std::vector<Real> > ahwvp
311       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
312     ROL::Ptr<const std::vector<Real> > wp
313       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
314     ROL::Ptr<const std::vector<Real> > vp
315       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
316     ROL::Ptr<const std::vector<Real> > up
317       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
318     ROL::Ptr<const std::vector<Real> > zp
319       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
320     (*ahwvp)[0] = static_cast<Real>(0);
321   }
322 
applyAdjointHessian_21(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)323   void applyAdjointHessian_21(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
324                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
325     ROL::Ptr<std::vector<Real> > ahwvp
326       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
327     ROL::Ptr<const std::vector<Real> > wp
328       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
329     ROL::Ptr<const std::vector<Real> > vp
330       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
331     ROL::Ptr<const std::vector<Real> > up
332       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
333     ROL::Ptr<const std::vector<Real> > zp
334       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
335     (*ahwvp)[0] = static_cast<Real>(0);
336   }
337 
applyAdjointHessian_22(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & u,const ROL::Vector<Real> & z,Real & tol)338   void applyAdjointHessian_22(ROL::Vector<Real> &ahwv, const ROL::Vector<Real> &w, const ROL::Vector<Real> &v,
339                               const ROL::Vector<Real> &u, const ROL::Vector<Real> &z, Real &tol) {
340     ROL::Ptr<std::vector<Real> > ahwvp
341       = dynamic_cast<ROL::StdVector<Real>&>(ahwv).getVector();
342     ROL::Ptr<const std::vector<Real> > wp
343       = dynamic_cast<const ROL::StdVector<Real>&>(w).getVector();
344     ROL::Ptr<const std::vector<Real> > vp
345       = dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
346     ROL::Ptr<const std::vector<Real> > up
347       = dynamic_cast<const ROL::StdVector<Real>&>(u).getVector();
348     ROL::Ptr<const std::vector<Real> > zp
349       = dynamic_cast<const ROL::StdVector<Real>&>(z).getVector();
350     (*ahwvp)[0] = static_cast<Real>(-2) * (*wp)[0] * (*vp)[0];
351   }
352 };
353 
354 typedef double RealT;
355 
main(int argc,char * argv[])356 int main(int argc, char *argv[]) {
357 
358   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
359 
360   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
361   int iprint     = argc - 1;
362   ROL::Ptr<std::ostream> outStream;
363   ROL::nullstream bhs; // outputs nothing
364   if (iprint > 0)
365     outStream = ROL::makePtrFromRef(std::cout);
366   else
367     outStream = ROL::makePtrFromRef(bhs);
368 
369   int errorFlag  = 0;
370 
371   // *** Example body.
372 
373   try {
374 
375     int dim = 2;
376     int dimz = 1;
377     ROL::Ptr<std::vector<RealT> > ustd  = ROL::makePtr<std::vector<RealT>>(dim);
378     ROL::Ptr<std::vector<RealT> > dustd = ROL::makePtr<std::vector<RealT>>(dim);
379     ROL::Ptr<std::vector<RealT> > zstd  = ROL::makePtr<std::vector<RealT>>(dimz);
380     ROL::Ptr<std::vector<RealT> > dzstd = ROL::makePtr<std::vector<RealT>>(dimz);
381     ROL::Ptr<std::vector<RealT> > cstd  = ROL::makePtr<std::vector<RealT>>(dim);
382     ROL::Ptr<std::vector<RealT> > czstd = ROL::makePtr<std::vector<RealT>>(dimz);
383     ROL::Ptr<std::vector<RealT> > sstd  = ROL::makePtr<std::vector<RealT>>(dimz);
384     ROL::Ptr<std::vector<RealT> > dsstd = ROL::makePtr<std::vector<RealT>>(dimz);
385 
386     (*ustd)[0]  = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
387     (*ustd)[1]  = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
388     (*dustd)[0] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
389     (*dustd)[1] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
390     (*zstd)[0]  = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
391     (*dzstd)[0] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
392     (*cstd)[0]  = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
393     (*cstd)[1]  = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
394     (*czstd)[0] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
395     (*sstd)[0]  = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
396     (*dsstd)[0] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
397 
398     ROL::Ptr<ROL::Vector<RealT> > u  = ROL::makePtr<ROL::StdVector<RealT>>(ustd);
399     ROL::Ptr<ROL::Vector<RealT> > du = ROL::makePtr<ROL::StdVector<RealT>>(dustd);
400     ROL::Ptr<ROL::Vector<RealT> > z  = ROL::makePtr<ROL::StdVector<RealT>>(zstd);
401     ROL::Ptr<ROL::Vector<RealT> > dz = ROL::makePtr<ROL::StdVector<RealT>>(dzstd);
402     ROL::Ptr<ROL::Vector<RealT> > c  = ROL::makePtr<ROL::StdVector<RealT>>(cstd);
403     ROL::Ptr<ROL::Vector<RealT> > cz = ROL::makePtr<ROL::StdVector<RealT>>(czstd);
404     ROL::Ptr<ROL::Vector<RealT> > s  = ROL::makePtr<ROL::StdVector<RealT>>(sstd);
405     ROL::Ptr<ROL::Vector<RealT> > ds = ROL::makePtr<ROL::StdVector<RealT>>(dsstd);
406 
407     ROL::Vector_SimOpt<RealT> x(u,s);
408     ROL::Vector_SimOpt<RealT> dx(du,ds);
409     ROL::Vector_SimOpt<RealT> y(s,z);
410     ROL::Vector_SimOpt<RealT> dy(ds,dz);
411     ROL::Vector_SimOpt<RealT> w(u,z);
412     ROL::Vector_SimOpt<RealT> dw(du,dz);
413 
414     ROL::Ptr<ROL::Constraint_SimOpt<RealT> > valCon = ROL::makePtr<valConstraint<RealT>>();
415     valCon->checkAdjointConsistencyJacobian_1(*c,*du,*u,*s,true,*outStream);
416     valCon->checkAdjointConsistencyJacobian_2(*c,*dz,*u,*s,true,*outStream);
417     valCon->checkApplyJacobian_1(*u,*s,*du,*c,true,*outStream);
418     valCon->checkApplyJacobian_2(*u,*s,*ds,*c,true,*outStream);
419     valCon->checkApplyJacobian(x,dx,*c,true,*outStream);
420     valCon->checkApplyAdjointHessian_11(*u,*s,*c,*du,*u,true,*outStream);
421     valCon->checkApplyAdjointHessian_12(*u,*s,*c,*du,*s,true,*outStream);
422     valCon->checkApplyAdjointHessian_21(*u,*s,*c,*ds,*u,true,*outStream);
423     valCon->checkApplyAdjointHessian_22(*u,*s,*c,*ds,*s,true,*outStream);
424     valCon->checkApplyAdjointHessian(x,*c,dx,x,true,*outStream);
425 
426     ROL::Ptr<ROL::Constraint_SimOpt<RealT> > redCon = ROL::makePtr<redConstraint<RealT>>();
427     redCon->checkAdjointConsistencyJacobian_1(*cz,*ds,*s,*z,true,*outStream);
428     redCon->checkAdjointConsistencyJacobian_2(*cz,*dz,*s,*z,true,*outStream);
429     redCon->checkInverseJacobian_1(*cz,*ds,*s,*z,true,*outStream);
430     redCon->checkInverseAdjointJacobian_1(*ds,*cz,*s,*z,true,*outStream);
431     redCon->checkApplyJacobian_1(*s,*z,*ds,*cz,true,*outStream);
432     redCon->checkApplyJacobian_2(*s,*z,*dz,*cz,true,*outStream);
433     redCon->checkApplyJacobian(y,dy,*cz,true,*outStream);
434     redCon->checkApplyAdjointHessian_11(*s,*z,*cz,*ds,*s,true,*outStream);
435     redCon->checkApplyAdjointHessian_12(*s,*z,*cz,*ds,*z,true,*outStream);
436     redCon->checkApplyAdjointHessian_21(*s,*z,*cz,*dz,*s,true,*outStream);
437     redCon->checkApplyAdjointHessian_22(*s,*z,*cz,*dz,*z,true,*outStream);
438     redCon->checkApplyAdjointHessian(y,*cz,dy,y,true,*outStream);
439 
440     ROL::CompositeConstraint_SimOpt<RealT> con(valCon,redCon,*c,*cz,*u,*s,*z);
441     con.checkAdjointConsistencyJacobian_1(*c,*du,*u,*z,true,*outStream);
442     con.checkAdjointConsistencyJacobian_2(*c,*dz,*u,*z,true,*outStream);
443     con.checkApplyJacobian_1(*u,*z,*du,*c,true,*outStream);
444     con.checkApplyJacobian_2(*u,*z,*dz,*c,true,*outStream);
445     con.checkApplyJacobian(w,dw,*c,true,*outStream);
446     con.checkApplyAdjointHessian_11(*u,*z,*c,*du,*u,true,*outStream);
447     con.checkApplyAdjointHessian_12(*u,*z,*c,*du,*z,true,*outStream);
448     con.checkApplyAdjointHessian_21(*u,*z,*c,*dz,*u,true,*outStream);
449     con.checkApplyAdjointHessian_22(*u,*z,*c,*dz,*z,true,*outStream);
450     con.checkApplyAdjointHessian(w,*c,dw,w,true,*outStream);
451   }
452   catch (std::logic_error& err) {
453     *outStream << err.what() << "\n";
454     errorFlag = -1000;
455   }; // end try
456 
457   if (errorFlag != 0)
458     std::cout << "End Result: TEST FAILED\n";
459   else
460     std::cout << "End Result: TEST PASSED\n";
461 
462   return 0;
463 
464 }
465 
466