ROL
function/test_15.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
50#include "ROL_Bounds.hpp"
52#include "ROL_StdConstraint.hpp"
53
54#include "ROL_Stream.hpp"
55#include "Teuchos_GlobalMPISession.hpp"
56
57template<typename Real>
58class con2d : public ROL::StdConstraint<Real> {
59public:
60 void value(std::vector<Real> &c, const std::vector<Real> &x, Real &tol) {
61 c[0] = x[0]+x[1];
62 }
63 void applyJacobian(std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol) {
64 jv[0] = v[0]+v[1];
65 }
66 void applyAdjointJacobian(std::vector<Real> &ajv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol) {
67 ajv[0] = v[0];
68 ajv[1] = v[0];
69 }
70};
71
72typedef double RealT;
73
74int main(int argc, char *argv[]) {
75
76 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77
78 // This little trick lets us print to std::cout only if a
79 // (dummy) command-line argument is provided.
80 int iprint = argc - 1;
81 ROL::Ptr<std::ostream> outStream;
82 ROL::nullstream bhs; // outputs nothing
83 if (iprint > 0)
84 outStream = ROL::makePtrFromRef(std::cout);
85 else
86 outStream = ROL::makePtrFromRef(bhs);
87
88 int errorFlag = 0;
89
90 try {
91 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
92 RealT err(0);
93 ROL::Ptr<con2d<RealT>> con = ROL::makePtr<con2d<RealT>>();
94
95 ROL::Ptr<std::vector<RealT>> yptr = ROL::makePtr<std::vector<RealT>>(2);
96 (*yptr)[0] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
97 (*yptr)[1] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
99
101
102 ROL::Ptr<std::vector<RealT>> xptr = ROL::makePtr<std::vector<RealT>>(2);
103 (*xptr)[0] = (*yptr)[0];
104 (*xptr)[1] = (*yptr)[1];
105 ROL::StdVector<RealT> x(xptr);
106
107 ROL::Ptr<std::vector<RealT>> Pxptr = ROL::makePtr<std::vector<RealT>>(2,0.0);
108 ROL::StdVector<RealT> Px(Pxptr);
109
111 ns0.apply(Px,x,tol);
112
113 ROL::Ptr<std::vector<RealT>> x0ptr = ROL::makePtr<std::vector<RealT>>(2);
114 (*x0ptr)[0] = ((*yptr)[0]-(*yptr)[1])/static_cast<RealT>(2);
115 (*x0ptr)[1] = -(*x0ptr)[0];
116 ROL::StdVector<RealT> x0(x0ptr);
117
119
120 *outStream << std::setprecision(6) << std::scientific << std::endl;
121 *outStream << " x[0] = " << (*xptr)[0] << " x[1] = " << (*xptr)[1] << std::endl;
122 *outStream << " Px[0] = " << (*Pxptr)[0] << " Px[1] = " << (*Pxptr)[1] << std::endl;
123 *outStream << " x*[0] = " << (*x0ptr)[0] << " x*[1] = " << (*x0ptr)[1] << std::endl;
124
125 e0.set(x0); e0.axpy(static_cast<RealT>(-1),Px);
126 err = e0.norm();
127 *outStream << " Error in Euclidean Projection: " << err << std::endl;
128
129 e0.set(x); e0.axpy(static_cast<RealT>(-1),x0);
130 *outStream << " ||x*-x||^2 = " << e0.norm() << std::endl;
131
132 e0.set(x); e0.axpy(static_cast<RealT>(-1),Px);
133 *outStream << " ||Px-x||^2 = " << e0.norm() << std::endl << std::endl;
134
135 errorFlag += (err > tol);
136
137 ROL::Ptr<std::vector<RealT>> dptr = ROL::makePtr<std::vector<RealT>>(2);
138 (*dptr)[0] = static_cast<RealT>(1)+static_cast<RealT>(2)*static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
139 (*dptr)[1] = static_cast<RealT>(1)+static_cast<RealT>(5)*static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
140
141 ROL::Ptr<std::vector<RealT>> x1ptr = ROL::makePtr<std::vector<RealT>>(2);
142 (*x1ptr)[0] = ((*dptr)[0]*(*yptr)[0]-(*dptr)[1]*(*yptr)[1])/((*dptr)[0]+(*dptr)[1]);
143 (*x1ptr)[1] = -(*x1ptr)[0];
145
146 ROL::Ptr<std::vector<RealT>> zptr = ROL::makePtr<std::vector<RealT>>(2);
147 (*zptr)[0] = (*yptr)[0];
148 (*zptr)[1] = (*yptr)[1];
150
151 ROL::Ptr<std::vector<RealT>> Pzptr = ROL::makePtr<std::vector<RealT>>(2,0.0);
153
155 ns1.apply(Pz,z,tol);
156
157 ROL::Ptr<std::vector<RealT>> e1ptr = ROL::makePtr<std::vector<RealT>>(2);
159
160 *outStream << std::endl;
161 *outStream << " x[0] = " << (*zptr)[0] << " x[1] = " << (*zptr)[1] << std::endl;
162 *outStream << " Px[0] = " << (*Pzptr)[0] << " Px[1] = " << (*Pzptr)[1] << std::endl;
163 *outStream << " x*[0] = " << (*x1ptr)[0] << " x*[1] = " << (*x1ptr)[1] << std::endl;
164
165 e1.set(x1); e1.axpy(static_cast<RealT>(-1),Pz);
166 err = e1.norm();
167 *outStream << " Error in Euclidean Projection: " << err << std::endl;
168
169 e1.set(z); e1.axpy(static_cast<RealT>(-1),x1);
170 *outStream << " ||x*-x||^2 = " << e1.norm() << std::endl;
171
172 e1.set(z); e1.axpy(static_cast<RealT>(-1),Pz);
173 *outStream << " ||Px-x||^2 = " << e1.norm() << std::endl << std::endl;
174
175 errorFlag += (err > tol);
176 }
177
178 catch (std::logic_error& err) {
179 *outStream << err.what() << "\n";
180 errorFlag = -1000;
181 }; // end try
182
183 if (errorFlag != 0)
184 std::cout << "End Result: TEST FAILED\n";
185 else
186 std::cout << "End Result: TEST PASSED\n";
187
188 return 0;
189}
190
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Projects on to the null space of a linear constraint.
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Provides the std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
Defines the equality constraint operator interface for StdVectors.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
void axpy(const Real alpha, const Vector< Real > &x)
Compute where .
void set(const Vector< Real > &x)
Set where .
Real norm() const
Returns where .
void applyJacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
void value(std::vector< Real > &c, const std::vector< Real > &x, Real &tol)
void applyAdjointJacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
int main(int argc, char *argv[])
double RealT