ROL
ROL_PrimalDualInteriorPointReducedResidual.hpp
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
44#ifndef ROL_PRIMALDUALINTERIORPOINTREDUCEDRESIDUAL_H
45#define ROL_PRIMALDUALINTERIORPOINTREDUCEDRESIDUAL_H
46
48#include "ROL_Constraint.hpp"
49
50#include <iostream>
51
108namespace ROL {
109
110template<class Real>
112
113 typedef ROL::ParameterList PL;
114
122
123 typedef Elementwise::ValueSet<Real> ValueSet;
124
125 typedef typename PV::size_type size_type;
126
127private:
128
129 static const size_type OPT = 0;
130 static const size_type EQUAL = 1;
131 static const size_type LOWER = 2;
132 static const size_type UPPER = 3;
133
134 ROL::Ptr<const V> x_; // Optimization vector
135 ROL::Ptr<const V> l_; // constraint multiplier
136 ROL::Ptr<const V> zl_; // Lower bound multiplier
137 ROL::Ptr<const V> zu_; // Upper bound multiplier
138
139 ROL::Ptr<const V> xl_; // Lower bound
140 ROL::Ptr<const V> xu_; // Upper bound
141
142 const ROL::Ptr<const V> maskL_;
143 const ROL::Ptr<const V> maskU_;
144
145 ROL::Ptr<V> scratch_;
146
147 const ROL::Ptr<PENALTY> penalty_;
148 const ROL::Ptr<OBJ> obj_;
149 const ROL::Ptr<CON> con_;
150
151
152public:
153
154 PrimalDualInteriorPointResidual( const ROL::Ptr<PENALTY> &penalty,
155 const ROL::Ptr<CON> &con,
156 const V &x,
157 ROL::Ptr<V> &scratch ) :
158 penalty_(penalty), con_(con), scratch_(scratch) {
159
160 obj_ = penalty_->getObjective();
161 maskL_ = penalty_->getLowerMask();
162 maskU_ = penalty_->getUpperMask();
163
164 ROL::Ptr<BND> bnd = penalty_->getBoundConstraint();
165 xl_ = bnd->getLowerBound();
166 xu_ = bnd->getUpperBound();
167
168
169 // Get access to the four components
170 const PV &x_pv = dynamic_cast<const PV&>(x);
171
172 x_ = x_pv.get(OPT);
173 l_ = x_pv.get(EQUAL);
174 zl_ = x_pv.get(LOWER);
175 zu_ = x_pv.get(UPPER);
176
177 }
178
179
180 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
181
182 // Get access to the four components
183 const PV &x_pv = dynamic_cast<const PV&>(x);
184
185 x_ = x_pv.get(OPT);
186 l_ = x_pv.get(EQUAL);
187 zl_ = x_pv.get(LOWER);
188 zu_ = x_pv.get(UPPER);
189
190 obj_->update(*x_,flag,iter);
191 con_->update(*x_,flag,iter);
192 }
193
194
195 // Evaluate the gradient of the Lagrangian
196 void value( V &c, const V &x, Real &tol ) {
197
198
199 PV &c_pv = dynamic_cast<PV&>(c);
200 const PV &x_pv = dynamic_cast<const PV&>(x);
201
202 x_ = x_pv.get(OPT);
203 l_ = x_pv.get(EQUAL);
204 zl_ = x_pv.get(LOWER);
205 zu_ = x_pv.get(UPPER);
206
207 ROL::Ptr<V> cx = c_pv.get(OPT);
208 ROL::Ptr<V> cl = c_pv.get(EQUAL);
209
210 // TODO: Add check as to whether we really need to recompute these
211 penalty_->gradient(*cx,*x_,tol);
212 con_->applyAdjointJacobian(*scratch_,*l_,*x_,tol);
213
214 // \f$ c_x = \nabla \phi_mu(x) + J^\ast \lambda \f$
215 cx_->plus(*scratch_);
216
217 con_->value(*cl_,*x_,tol);
218
219 }
220
221 void applyJacobian( V &jv, const V &v, const V &x, Real &tol ) {
222
223
224
225 PV &jv_pv = dynamic_cast<PV&>(jv);
226 const PV &v_pv = dynamic_cast<const PV&>(v);
227 const PV &x_pv = dynamic_cast<const PV&>(x);
228
229 // output vector components
230 ROL::Ptr<V> jvx = jv_pv.get(OPT);
231 ROL::Ptr<V> jvl = jv_pv.get(EQUAL);
232
233 // input vector components
234 ROL::Ptr<const V> vx = v_pv.get(OPT);
235 ROL::Ptr<const V> vl = v_pv.get(EQUAL);
236
237 x_ = x_pv.get(OPT);
238 l_ = x_pv.get(EQUAL);
239
240 // \f$
241 obj_->hessVec(*jvx,*vx,*x_,tol);
242 con_->applyAdjointHessian(*scratch_,*l_,*vx,*x_,tol);
243 jvx->plus(*scratch_);
244
245 // \f$J^\ast d_\lambda \f$
246 con_->applyAdjointJacobian(*scratch_,*vl,*x_,tol);
247 jvx->plus(*scratch_);
248
249
250 Elementwise::DivideAndInvert<Real> divinv;
251 Elementwise::Multiply<Real> mult;
252
253 // Note that indices corresponding to infinite bounds should automatically lead to
254 // zero diagonal contributions to the Sigma operators
255 scratch_->set(*x_);
256 scratch_->axpy(-1.0,*xl_); // x-l
257 scratch_->applyBinary(divinv,*zl_); // zl/(x-l)
258 scratch_->applyBinary(mult,*vx); // zl*vx/(x-l) = Sigma_l*vx
259
260 jvx->plu(*scratch_);
261
262 scratch_->set(*xu_);
263 scratch_->axpy(-1.0,*x_); // u-x
264 scratch_->applyBinary(divinv,*zu_); // zu/(u-x)
265 scratch_->applyBinary(mult,*vx); // zu*vx/(u-x) = Sigma_u*vx
266
267 jvx->plus(*scratch_);
268
269
270
271
272
273
274
275// \f[ [W+(X-L)^{-1}+(U-X)^{-1}]d_x + J^\ast d_\lambda =
276// -g_x + (U-X)^{-1}g_{z_u} - (L-X)^{-1}g_{z_l} \f]
277
278
279
280 }
281
283 return nfval_;
284 }
285
287 return ngrad_;
288 }
289
291 return ncval_;
292 }
293
294
295}; // class PrimalDualInteriorPointResidual
296
297} // namespace ROL
298
299#endif // ROL_PRIMALDUALKKTOPERATOR_H
Provides the interface to apply upper and lower bound constraints.
Defines the general constraint operator interface.
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
Provides the interface to apply a linear operator.
Provides the interface to evaluate objective functions.
Defines the linear algebra of vector space on a generic partitioned vector.
std::vector< PV >::size_type size_type
ROL::Ptr< const Vector< Real > > get(size_type i) const
Symmetrized form of the KKT operator for the Type-EB problem with equality and bound multipliers.
void applyJacobian(V &jv, const V &v, const V &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable i...
void value(V &c, const V &x, Real &tol)
Evaluate the constraint operator at .
PrimalDualInteriorPointResidual(const ROL::Ptr< PENALTY > &penalty, const ROL::Ptr< CON > &con, const V &x, ROL::Ptr< V > &scratch)
Defines the linear algebra or vector space interface.