51 #ifndef ROL_SIMPLEEQCONSTRAINED_HPP 52 #define ROL_SIMPLEEQCONSTRAINED_HPP 57 #include "Teuchos_SerialDenseVector.hpp" 58 #include "Teuchos_SerialDenseSolver.hpp" 66 template<
class Real,
class XPrim=StdVector<Real>,
class XDual=StdVector<Real> >
72 typedef typename vector::size_type
uint;
77 template<
class VectorType>
79 using Teuchos::dyn_cast;
80 return dyn_cast<
const VectorType>(x).
getVector();
83 template<
class VectorType>
85 using Teuchos::dyn_cast;
86 return dyn_cast<VectorType>(x).
getVector();
95 RCP<const vector> xp = getVector<XPrim>(x);
98 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective value): " 99 "Primal vector x must be of length 5.");
107 Real val = exp(x1*x2*x3*x4*x5) - 0.5 * pow( (pow(x1,3)+pow(x2,3)+1.0), 2);
115 RCP<const vector> xp = getVector<XPrim>(x);
116 RCP<vector> gp = getVector<XDual>(g);
119 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective gradient): " 120 " Primal vector x must be of length 5.");
123 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective gradient): " 124 "Gradient vector g must be of length 5.");
132 Real expxi = exp(x1*x2*x3*x4*x5);
134 (*gp)[0] = x2*x3*x4*x5 * expxi - 3*pow(x1,2) * (pow(x1,3) + pow(x2,3) + 1);
135 (*gp)[1] = x1*x3*x4*x5 * expxi - 3*pow(x2,2) * (pow(x1,3) + pow(x2,3) + 1);
136 (*gp)[2] = x1*x2*x4*x5 * expxi;
137 (*gp)[3] = x1*x2*x3*x5 * expxi;
138 (*gp)[4] = x1*x2*x3*x4 * expxi;
144 RCP<const vector> xp = getVector<XPrim>(x);
145 RCP<const vector> vp = getVector<XPrim>(v);
146 RCP<vector> hvp = getVector<XDual>(hv);
149 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective hessVec): " 150 "Primal vector x must be of length 5.");
153 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective hessVec): " 154 "Input vector v must be of length 5.");
157 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, objective hessVec): " 158 "Output vector hv must be of length 5.");
172 Real expxi = exp(x1*x2*x3*x4*x5);
174 (*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 +
175 ( 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 +
176 ( x2*x4*x5*expxi+pow(x2,2)*x3*pow(x4,2)*pow(x5,2)*x1*expxi ) * v3 +
177 ( x2*x3*x5*expxi+pow(x2,2)*pow(x3,2)*x4*pow(x5,2)*x1*expxi ) * v4 +
178 ( x2*x3*x4*expxi+pow(x2,2)*pow(x3,2)*pow(x4,2)*x5*x1*expxi ) * v5;
180 (*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 +
181 ( 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 +
182 ( x1*x4*x5*expxi+pow(x1,2)*x3*pow(x4,2)*pow(x5,2)*x2*expxi ) * v3 +
183 ( x1*x3*x5*expxi+pow(x1,2)*pow(x3,2)*x4*pow(x5,2)*x2*expxi ) * v4 +
184 ( x1*x3*x4*expxi+pow(x1,2)*pow(x3,2)*pow(x4,2)*x5*x2*expxi ) * v5;
186 (*hvp)[2] = ( x2*x4*x5*expxi+pow(x2,2)*x3*pow(x4,2)*pow(x5,2)*x1*expxi ) * v1 +
187 ( x1*x4*x5*expxi+pow(x1,2)*x3*pow(x4,2)*pow(x5,2)*x2*expxi ) * v2 +
188 ( pow(x1,2)*pow(x2,2)*pow(x4,2)*pow(x5,2)*expxi ) * v3 +
189 ( x1*x2*x5*expxi+pow(x1,2)*pow(x2,2)*x4*pow(x5,2)*x3*expxi ) * v4 +
190 ( x1*x2*x4*expxi+pow(x1,2)*pow(x2,2)*pow(x4,2)*x5*x3*expxi ) * v5;
192 (*hvp)[3] = ( x2*x3*x5*expxi+pow(x2,2)*pow(x3,2)*x4*pow(x5,2)*x1*expxi ) * v1 +
193 ( x1*x3*x5*expxi+pow(x1,2)*pow(x3,2)*x4*pow(x5,2)*x2*expxi ) * v2 +
194 ( x1*x2*x5*expxi+pow(x1,2)*pow(x2,2)*x4*pow(x5,2)*x3*expxi ) * v3 +
195 ( pow(x1,2)*pow(x2,2)*pow(x3,2)*pow(x5,2)*expxi ) * v4 +
196 ( x1*x2*x3*expxi+pow(x1,2)*pow(x2,2)*pow(x3,2)*x5*x4*expxi ) * v5;
198 (*hvp)[4] = ( x2*x3*x4*expxi+pow(x2,2)*pow(x3,2)*pow(x4,2)*x5*x1*expxi ) * v1 +
199 ( x1*x3*x4*expxi+pow(x1,2)*pow(x3,2)*pow(x4,2)*x5*x2*expxi ) * v2 +
200 ( x1*x2*x4*expxi+pow(x1,2)*pow(x2,2)*pow(x4,2)*x5*x3*expxi ) * v3 +
201 ( x1*x2*x3*expxi+pow(x1,2)*pow(x2,2)*pow(x3,2)*x5*x4*expxi ) * v4 +
202 ( pow(x1,2)*pow(x2,2)*pow(x3,2)*pow(x4,2)*expxi ) * v5;
213 template<
class Real,
class XPrim=StdVector<Real>,
class XDual=StdVector<Real>,
class CPrim=StdVector<Real>,
class CDual=StdVector<Real> >
219 typedef typename vector::size_type
uint;
222 template<
class VectorType>
224 using Teuchos::dyn_cast;
225 return dyn_cast<
const VectorType>(x).
getVector();
228 template<
class VectorType>
230 using Teuchos::dyn_cast;
231 return dyn_cast<VectorType>(x).
getVector();
240 RCP<const vector> xp = getVector<XPrim>(x);
241 RCP<vector> cp = getVector<CPrim>(c);
244 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint value): " 245 "Primal vector x must be of length 5.");
248 TEUCHOS_TEST_FOR_EXCEPTION( (m != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint value): " 249 "Constraint vector c must be of length 3.");
257 (*cp)[0] = x1*x1+x2*x2+x3*x3+x4*x4+x5*x5 - 10.0;
258 (*cp)[1] = x2*x3 - 5.0*x4*x5;
259 (*cp)[2] = x1*x1*x1 + x2*x2*x2 + 1.0;
265 RCP<const vector> xp = getVector<XPrim>(x);
266 RCP<const vector> vp = getVector<XPrim>(v);
267 RCP<vector> jvp = getVector<CPrim>(jv);
270 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyJacobian): " 271 "Primal vector x must be of length 5.");
274 TEUCHOS_TEST_FOR_EXCEPTION( (d != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyJacobian): " 275 "Input vector v must be of length 5.");
277 TEUCHOS_TEST_FOR_EXCEPTION( (d != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyJacobian): " 278 "Output vector jv must be of length 3.");
292 (*jvp)[0] = 2.0*(x1*v1+x2*v2+x3*v3+x4*v4+x5*v5);
293 (*jvp)[1] = x3*v2+x2*v3-5.0*x5*v4-5.0*x4*v5;
294 (*jvp)[2] = 3.0*x1*x1*v1+3.0*x2*x2*v2;
301 RCP<const vector> xp = getVector<XPrim>(x);
302 RCP<const vector> vp = getVector<CDual>(v);
303 RCP<vector> ajvp = getVector<XDual>(ajv);
306 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointJacobian): " 307 "Primal vector x must be of length 5.");
310 TEUCHOS_TEST_FOR_EXCEPTION( (d != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointJacobian): " 311 "Input vector v must be of length 3.");
314 TEUCHOS_TEST_FOR_EXCEPTION( (d != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointJacobian): " 315 "Output vector ajv must be of length 5.");
327 (*ajvp)[0] = 2.0*x1*v1+3.0*x1*x1*v3;
328 (*ajvp)[1] = 2.0*x2*v1+x3*v2+3.0*x2*x2*v3;
329 (*ajvp)[2] = 2.0*x3*v1+x2*v2;
330 (*ajvp)[3] = 2.0*x4*v1-5.0*x5*v2;
331 (*ajvp)[4] = 2.0*x5*v1-5.0*x4*v2;
337 RCP<const vector> xp = getVector<XPrim>(x);
338 RCP<const vector> up = getVector<CDual>(u);
339 RCP<const vector> vp = getVector<XPrim>(v);
340 RCP<vector> ahuvp = getVector<XDual>(ahuv);
343 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): " 344 "Primal vector x must be of length 5.");
347 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): " 348 "Direction vector v must be of length 5.");
351 TEUCHOS_TEST_FOR_EXCEPTION( (n != 5), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): " 352 "Output vector ahuv must be of length 5.");
354 TEUCHOS_TEST_FOR_EXCEPTION( (d != 3), std::invalid_argument,
">>> ERROR (ROL_SimpleEqConstrained, constraint applyAdjointHessian): " 355 "Dual constraint vector u must be of length 3.");
370 (*ahuvp)[0] = 2.0*u1*v1 + 6.0*u3*x1*v1;
371 (*ahuvp)[1] = 2.0*u1*v2 + u2*v3 + 6.0*u3*x2*v2;
372 (*ahuvp)[2] = 2.0*u1*v3 + u2*v2;
373 (*ahuvp)[3] = 2.0*u1*v4 - 5.0*u2*v5;
374 (*ahuvp)[4] = 2.0*u1*v5 - 5.0*u2*v4;
439 template<
class Real,
class XPrim,
class XDual,
class CPrim,
class CDual>
445 typedef std::vector<Real> vector;
447 typedef typename vector::size_type uint;
449 using Teuchos::RCP;
using Teuchos::rcp;
450 using Teuchos::dyn_cast;
453 RCP<vector> x0p = dyn_cast<XPrim>(x0).getVector();
454 RCP<vector> solp = dyn_cast<XPrim>(sol).getVector();
476 (*solp)[0] = -1.717143570394391e+00;
477 (*solp)[1] = 1.595709690183565e+00;
478 (*solp)[2] = 1.827245752927178e+00;
479 (*solp)[3] = -7.636430781841294e-01;
480 (*solp)[4] = -7.636430781841294e-01;
Provides the interface to evaluate objective functions.
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 .
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
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 .
Defines the equality constraint operator interface.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::RCP< vector > getVector(V &x)
Teuchos::RCP< vector > getVector(V &x)
Teuchos::RCP< const vector > getVector(const V &x)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
std::vector< Real > vector
void getSimpleEqConstrained(Teuchos::RCP< Objective< Real > > &obj, Teuchos::RCP< EqualityConstraint< Real > > &constr, Vector< Real > &x0, Vector< Real > &sol)
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 ...
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Objective_SimpleEqConstrained()
Teuchos::RCP< const vector > getVector(const V &x)
Objective function: f(x) = exp(x1*x2*x3*x4*x5) + 0.5*(x1^3+x2^3+1)^2.
std::vector< Real > vector