17 #ifndef ROL_SIMPLEEQCONSTRAINED_HPP
18 #define ROL_SIMPLEEQCONSTRAINED_HPP
22 #include "Teuchos_SerialDenseVector.hpp"
23 #include "Teuchos_SerialDenseSolver.hpp"
31 template<
class Real,
class XPrim=StdVector<Real>,
class XDual=StdVector<Real> >
42 template<
class VectorType>
45 return dynamic_cast<const VectorType&
>(x).
getVector();
48 template<
class VectorType>
51 return dynamic_cast<VectorType&
>(x).
getVector();
60 ROL::Ptr<const vector> xp = getVector<XPrim>(x);
63 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective value): "
64 "Primal vector x must be of length 5.");
72 Real arg = x1*x2*x3*x4*x5;
73 Real val = -0.5*pow(pow(x1,3)+pow(x2,3)+1.0,2);
74 if (std::abs(arg) < 1e5) val += exp(x1*x2*x3*x4*x5);
75 else if (arg > 1e5) val += 1e10;
84 ROL::Ptr<const vector> xp = getVector<XPrim>(x);
85 ROL::Ptr<vector> gp = getVector<XDual>(g);
88 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective gradient): "
89 " Primal vector x must be of length 5.");
92 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective gradient): "
93 "Gradient vector g must be of length 5.");
101 Real expxi = exp(x1*x2*x3*x4*x5);
103 (*gp)[0] = x2*x3*x4*x5 * expxi - 3*pow(x1,2) * (pow(x1,3) + pow(x2,3) + 1);
104 (*gp)[1] = x1*x3*x4*x5 * expxi - 3*pow(x2,2) * (pow(x1,3) + pow(x2,3) + 1);
105 (*gp)[2] = x1*x2*x4*x5 * expxi;
106 (*gp)[3] = x1*x2*x3*x5 * expxi;
107 (*gp)[4] = x1*x2*x3*x4 * expxi;
113 ROL::Ptr<const vector> xp = getVector<XPrim>(x);
114 ROL::Ptr<const vector> vp = getVector<XPrim>(v);
115 ROL::Ptr<vector> hvp = getVector<XDual>(hv);
118 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective hessVec): "
119 "Primal vector x must be of length 5.");
122 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective hessVec): "
123 "Input vector v must be of length 5.");
126 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective hessVec): "
127 "Output vector hv must be of length 5.");
141 Real expxi = exp(x1*x2*x3*x4*x5);
143 (*hvp)[0] = ( pow(x2,2)*pow(x3,2)*pow(x4,2)*pow(x5,2)*expxi-9.0*pow(x1,4)-6.0*(pow(x1,3)+pow(x2,3)+1.0)*x1 ) * v1 +
144 ( x3*x4*x5*expxi+x2*pow(x3,2)*pow(x4,2)*pow(x5,2)*x1*expxi-9.0*pow(x2,2)*pow(x1,2) ) * v2 +
145 ( x2*x4*x5*expxi+pow(x2,2)*x3*pow(x4,2)*pow(x5,2)*x1*expxi ) * v3 +
146 ( x2*x3*x5*expxi+pow(x2,2)*pow(x3,2)*x4*pow(x5,2)*x1*expxi ) * v4 +
147 ( x2*x3*x4*expxi+pow(x2,2)*pow(x3,2)*pow(x4,2)*x5*x1*expxi ) * v5;
149 (*hvp)[1] = ( x3*x4*x5*expxi+x2*pow(x3,2)*pow(x4,2)*pow(x5,2)*x1*expxi-9.0*pow(x2,2)*pow(x1,2) ) * v1 +
150 ( pow(x1,2)*pow(x3,2)*pow(x4,2)*pow(x5,2)*expxi-9.0*pow(x2,4)-6.0*(pow(x1,3)+pow(x2,3)+1.0)*x2 ) * v2 +
151 ( x1*x4*x5*expxi+pow(x1,2)*x3*pow(x4,2)*pow(x5,2)*x2*expxi ) * v3 +
152 ( x1*x3*x5*expxi+pow(x1,2)*pow(x3,2)*x4*pow(x5,2)*x2*expxi ) * v4 +
153 ( x1*x3*x4*expxi+pow(x1,2)*pow(x3,2)*pow(x4,2)*x5*x2*expxi ) * v5;
155 (*hvp)[2] = ( x2*x4*x5*expxi+pow(x2,2)*x3*pow(x4,2)*pow(x5,2)*x1*expxi ) * v1 +
156 ( x1*x4*x5*expxi+pow(x1,2)*x3*pow(x4,2)*pow(x5,2)*x2*expxi ) * v2 +
157 ( pow(x1,2)*pow(x2,2)*pow(x4,2)*pow(x5,2)*expxi ) * v3 +
158 ( x1*x2*x5*expxi+pow(x1,2)*pow(x2,2)*x4*pow(x5,2)*x3*expxi ) * v4 +
159 ( x1*x2*x4*expxi+pow(x1,2)*pow(x2,2)*pow(x4,2)*x5*x3*expxi ) * v5;
161 (*hvp)[3] = ( x2*x3*x5*expxi+pow(x2,2)*pow(x3,2)*x4*pow(x5,2)*x1*expxi ) * v1 +
162 ( x1*x3*x5*expxi+pow(x1,2)*pow(x3,2)*x4*pow(x5,2)*x2*expxi ) * v2 +
163 ( x1*x2*x5*expxi+pow(x1,2)*pow(x2,2)*x4*pow(x5,2)*x3*expxi ) * v3 +
164 ( pow(x1,2)*pow(x2,2)*pow(x3,2)*pow(x5,2)*expxi ) * v4 +
165 ( x1*x2*x3*expxi+pow(x1,2)*pow(x2,2)*pow(x3,2)*x5*x4*expxi ) * v5;
167 (*hvp)[4] = ( x2*x3*x4*expxi+pow(x2,2)*pow(x3,2)*pow(x4,2)*x5*x1*expxi ) * v1 +
168 ( x1*x3*x4*expxi+pow(x1,2)*pow(x3,2)*pow(x4,2)*x5*x2*expxi ) * v2 +
169 ( x1*x2*x4*expxi+pow(x1,2)*pow(x2,2)*pow(x4,2)*x5*x3*expxi ) * v3 +
170 ( x1*x2*x3*expxi+pow(x1,2)*pow(x2,2)*pow(x3,2)*x5*x4*expxi ) * v4 +
171 ( pow(x1,2)*pow(x2,2)*pow(x3,2)*pow(x4,2)*expxi ) * v5;
182 template<
class Real,
class XPrim=StdVector<Real>,
class XDual=StdVector<Real>,
class CPrim=StdVector<Real>,
class CDual=StdVector<Real> >
191 template<
class VectorType>
194 return dynamic_cast<const VectorType&
>(x).
getVector();
197 template<
class VectorType>
200 return dynamic_cast<VectorType&
>(x).
getVector();
209 ROL::Ptr<const vector> xp = getVector<XPrim>(x);
210 ROL::Ptr<vector> cp = getVector<CPrim>(c);
213 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint value): "
214 "Primal vector x must be of length 5.");
217 ROL_TEST_FOR_EXCEPTION( (m != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint value): "
218 "Constraint vector c must be of length 3.");
226 (*cp)[0] = x1*x1+x2*x2+x3*x3+x4*x4+x5*x5 - 10.0;
227 (*cp)[1] = x2*x3 - 5.0*x4*x5;
228 (*cp)[2] = x1*x1*x1 + x2*x2*x2 + 1.0;
234 ROL::Ptr<const vector> xp = getVector<XPrim>(x);
235 ROL::Ptr<const vector> vp = getVector<XPrim>(v);
236 ROL::Ptr<vector> jvp = getVector<CPrim>(jv);
239 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyJacobian): "
240 "Primal vector x must be of length 5.");
243 ROL_TEST_FOR_EXCEPTION( (d != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyJacobian): "
244 "Input vector v must be of length 5.");
246 ROL_TEST_FOR_EXCEPTION( (d != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyJacobian): "
247 "Output vector jv must be of length 3.");
261 (*jvp)[0] = 2.0*(x1*v1+x2*v2+x3*v3+x4*v4+x5*v5);
262 (*jvp)[1] = x3*v2+x2*v3-5.0*x5*v4-5.0*x4*v5;
263 (*jvp)[2] = 3.0*x1*x1*v1+3.0*x2*x2*v2;
270 ROL::Ptr<const vector> xp = getVector<XPrim>(x);
271 ROL::Ptr<const vector> vp = getVector<CDual>(v);
272 ROL::Ptr<vector> ajvp = getVector<XDual>(ajv);
275 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointJacobian): "
276 "Primal vector x must be of length 5.");
279 ROL_TEST_FOR_EXCEPTION( (d != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointJacobian): "
280 "Input vector v must be of length 3.");
283 ROL_TEST_FOR_EXCEPTION( (d != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointJacobian): "
284 "Output vector ajv must be of length 5.");
296 (*ajvp)[0] = 2.0*x1*v1+3.0*x1*x1*v3;
297 (*ajvp)[1] = 2.0*x2*v1+x3*v2+3.0*x2*x2*v3;
298 (*ajvp)[2] = 2.0*x3*v1+x2*v2;
299 (*ajvp)[3] = 2.0*x4*v1-5.0*x5*v2;
300 (*ajvp)[4] = 2.0*x5*v1-5.0*x4*v2;
306 ROL::Ptr<const vector> xp = getVector<XPrim>(x);
307 ROL::Ptr<const vector> up = getVector<CDual>(u);
308 ROL::Ptr<const vector> vp = getVector<XPrim>(v);
309 ROL::Ptr<vector> ahuvp = getVector<XDual>(ahuv);
312 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): "
313 "Primal vector x must be of length 5.");
316 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): "
317 "Direction vector v must be of length 5.");
320 ROL_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): "
321 "Output vector ahuv must be of length 5.");
323 ROL_TEST_FOR_EXCEPTION( (d != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): "
324 "Dual constraint vector u must be of length 3.");
339 (*ahuvp)[0] = 2.0*u1*v1 + 6.0*u3*x1*v1;
340 (*ahuvp)[1] = 2.0*u1*v2 + u2*v3 + 6.0*u3*x2*v2;
341 (*ahuvp)[2] = 2.0*u1*v3 + u2*v2;
342 (*ahuvp)[3] = 2.0*u1*v4 - 5.0*u2*v5;
343 (*ahuvp)[4] = 2.0*u1*v5 - 5.0*u2*v4;
408 template<
class Real,
class XPrim=StdVector<Real>,
class XDual=StdVector<Real>,
class CPrim=StdVector<Real>,
class CDual=StdVector<Real> >
417 return ROL::makePtr<Objective_SimpleEqConstrained<Real,XPrim,XDual>>();
423 Ptr<vector> x0p = makePtr<vector>(n,0);
429 return makePtr<XPrim>(x0p);
435 Ptr<vector> solp = makePtr<vector>(n,0);
436 (*solp)[0] = -1.717143570394391e+00;
437 (*solp)[1] = 1.595709690183565e+00;
438 (*solp)[2] = 1.827245752927178e+00;
439 (*solp)[3] = -7.636430781841294e-01;
440 (*solp)[4] = -7.636430781841294e-01;
441 return makePtr<XPrim>(solp);
446 return ROL::makePtr<EqualityConstraint_SimpleEqConstrained<Real,XPrim,XDual,CPrim,CDual>>();
450 Ptr<vector> lp = makePtr<vector>(3,0);
451 return makePtr<CDual>(lp);
Provides the interface to evaluate objective functions.
typename PV< Real >::size_type size_type
EqualityConstraint_SimpleEqConstrained()
Equality constraints c_i(x) = 0, where: c1(x) = x1^2+x2^2+x3^2+x4^2+x5^2 - 10 c2(x) = x2*x3-5*x4*x5 c...
void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Ptr< Constraint< Real > > getEqualityConstraint(void) const
ROL::Ptr< const vector > getVector(const V &x)
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
Ptr< Vector< Real > > getEqualityMultiplier(void) const
Defines the linear algebra or vector space interface.
void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
std::vector< Real > vector
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
ROL::Ptr< vector > getVector(V &x)
ROL::Ptr< const vector > getVector(const V &x)
Contains definitions of test objective functions.
Ptr< Objective< Real > > getObjective(void) const
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
std::vector< Real > vector
ROL::Ptr< vector > getVector(V &x)
Ptr< Vector< Real > > getSolution(const int i=0) const
void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ...
getSimpleEqConstrained(void)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Objective_SimpleEqConstrained()
Ptr< Vector< Real > > getInitialGuess(void) const
Defines the general constraint operator interface.
Objective function: f(x) = exp(x1*x2*x3*x4*x5) + 0.5*(x1^3+x2^3+1)^2.
std::vector< Real > vector