ROL
ROL_PrimalDualInteriorPointOperator.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_PRIMALDUALINTERIORPOINTOPERATOR_H
45#define ROL_PRIMALDUALINTERIORPOINTOPERATOR_H
46
48
49
50namespace ROL {
51
52
53
54template<class Real>
56
57 typedef Vector<Real> V;
61
62 static const size_type OPT = 0;
63 static const size_type EQUAL = 1;
64 static const size_type LOWER = 0;
65 static const size_type UPPER = 1;
66
67 ROL::Ptr<const V> x_; // Optimization vector
68 ROL::Ptr<const V> l_; // constraint multiplier
69
70 ROL::Ptr<V> scratch_;
71
72 Real delta_; // Initial correction
73
74
75public:
76
77 PrimalDualInteriorPointBlock11( ROL::Ptr<OBJ> &obj, ROL::Ptr<CON> &con,
78 const V &x, ROL::Ptr<V> & scratch,
79 Real delta=0 ) :
80 obj_(obj), con_(con), scratch_(scratch), delta_(delta) {
81
82 const PV &x_pv = dynamic_cast<const PV&>(x);
83
84 x_ = x_pv.get(OPT);
85 l_ = x_pv.get(EQUAL);
86 }
87
88 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
89
90 const PV &x_pv = dynamic_cast<const PV&>(x);
91
92 x_ = x_pv.get(OPT);
93 l_ = x_pv.get(EQUAL);
94
95 obj_->update(*x_,flag,true);
96 con_->update(*x_,flag,true);
97 }
98
99 void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
100
101
102
103 PV &Hv_pv = dynamic_cast<PV&>(Hv);
104 const PV &v_pv = dynamic_cast<const PV&>(v);
105
106 // output vector components
107 ROL::Ptr<V> Hvx = Hv_pv.get(OPT);
108 ROL::Ptr<V> Hvl = Hv_pv.get(EQUAL);
109
110 // input vector components
111 ROL::Ptr<const V> vx = v_pv.get(OPT);
112 ROL::Ptr<const V> vl = v_pv.get(EQUAL);
113
114 obj_->hessVec(*jvx,*vx,*x_,tol);
115 con_->applyAdjointHessian(*scratch_,*l_,*vx,*x_,tol);
116 jvx->plus(*scratch_);
117 con_->applyAdjointJacobian(*scratch_,*vl,*x_,tol);
118 jvx->plus(*scratch_);
119
120 // Inertial correction
121 if( delta_ != 0 ) {
122 jvx->axpy(delta_,*vx);
123 }
124
125 con_->applyJacobian(*jvl,*vx,*x_,tol);
126
127 }
128
129 void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
130 ROL_TEST_FOR_EXCEPTION( true , std::logic_error,
131 ">>> ERROR (ROL_PrimalDualInteriorPointBlock11, applyInverse): "
132 "Not implemented.");
133 }
134
135 void setInertia( Real delta ) {
136 delta_ = delta;
137 }
138
139
140}; // class PrimalDualInteriorPointBlock11
141
142
143template<class Real>
145
148
149 static const size_type OPT = 0;
150 static const size_type EQUAL = 1;
151 static const size_type LOWER = 0;
152 static const size_type UPPER = 1;
153
154public:
155
156 void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
157
158
159 PV &Hv_pv = dynamic_cast<PV&>(Hv);
160 const PV &v_pv = dynamic_cast<const PV&>(v);
161
162 // output vector components
163 ROL::Ptr<V> Hvx = Hv_pv.get(OPT);
164 ROL::Ptr<V> Hvl = Hv_pv.get(EQUAL);
165
166 // input vector components
167 ROL::Ptr<const V> vzl = v_pv.get(LOWER);
168 ROL::Ptr<const V> vzu = v_pv.get(UPPER);
169
170 Hvx->set(*vzu);
171 Hvx->axpy(-1.0,*vzl);
172 Hvl->zero();
173
174 }
175
176 void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
177 ROL_TEST_FOR_EXCEPTION( true , std::logic_error,
178 ">>> ERROR (ROL_PrimalDualInteriorPointBlock12, applyInverse): "
179 "Not implemented.");
180 }
181
182
183}; // class PrimalDualInteriorPointBlock12
184
185
186template<class Real>
188
191
192 static const size_type OPT = 0;
193 static const size_type EQUAL = 1;
194 static const size_type LOWER = 0;
195 static const size_type UPPER = 1;
196
197 ROL::Ptr<const V> zl_;
198 ROL::Ptr<const V> zu_;
199
200public:
201
203 const PV &z_pv = dynamic_cast<const PV&>(z);
204 zl_ = z_pv.get(LOWER);
205 zu_ = z_pv.get(UPPER);
206 }
207
208 void update( const Vector<Real> &z, bool flag = true, int iter = -1 ) {
209 const PV &z_pv = dynamic_cast<const PV&>(z);
210 zl_ = z_pv.get(LOWER);
211 zu_ = z_pv.get(UPPER);
212 }
213
214 virtual void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
215
216
217 PV &Hv_pv = dynamic_cast<PV&>(Hv);
218 const PV &v_pv = dynamic_cast<const PV&>(v);
219
220 // output vector components
221 ROL::Ptr<V> Hvzl = Hv_pv.get(LOWER);
222 ROL::Ptr<V> Hvzu = Hv_pv.get(UPPER);
223
224 // input vector components
225 ROL::Ptr<const V> vx = v_pv.get(OPT);
226 ROL::Ptr<const V> vl = v_pv.get(EQUAL);
227
228 Hvzl->set(*vx);
229 Hvzl->applyBinary(mult_,*zl_);
230
231 Hvzu->set(*vx);
232 Hvzu->applyBinary(mult_,*zu_);
233 Hvzu->scale(-1.0);
234 }
235
236 virtual void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
237 ROL_TEST_FOR_EXCEPTION( true , std::logic_error,
238 ">>> ERROR (ROL_PrimalDualInteriorPointBlock21, applyInverse): "
239 "Not implemented.");
240 }
241
242}; // class PrimalDualInteriorPointBlock21
243
244
245template<class Real>
247
251
252 static const size_type OPT = 0;
253 static const size_type LOWER = 0;
254 static const size_type UPPER = 1;
255
256 ROL::Ptr<const V> x_;
257 ROL::Ptr<const V> xl_;
258 ROL::Ptr<const V> xu_;
259
260 Elementwise::Multiply<Real> mult_;
261 Elementwise::Multiply<Real> divinv_;
262
263
264public:
265
266 PrimalDualInteriorPointBlock22( const ROL::Ptr<BND> &bnd, const Vector<Real> &x ) {
267
268 const PV &x_pv = dynamic_cast<const PV&>(x);
269
270 x_ = x_pv.get(OPT);
271 xl_ = bnd.getLowerBound();
272 xu_ = bnd.getUpperBound();
273
274 }
275
276 virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
277
278 const PV &x_pv = dynamic_cast<const PV&>(x);
279 x_ = x_pv.get(OPT);
280 }
281
282 virtual void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
283
284
285 PV &Hv_pv = dynamic_cast<PV&>(Hv);
286 const PV &v_pv = dynamic_cast<const PV&>(v);
287
288 // output vector components
289 ROL::Ptr<V> Hvzl = Hv_pv.get(LOWER);
290 ROL::Ptr<V> Hvzu = Hv_pv.get(UPPER);
291
292 // input vector components
293 ROL::Ptr<const V> vzl = v_pv.get(LOWER);
294 ROL::Ptr<const V> vzu = v_pv.get(UPPER);
295
296 Hvzl->set(*x_);
297 Hvzl->axpy(-1.0,*xl_);
298 Hvzl->applyBinary(mult_,*vzl);
299
300 Hvzu->set(*xu_);
301 Hvzu->axpy(-1.0,*x_);
302 Hvzu->applyBinary(mult_,*vzu);
303
304 }
305
306 virtual void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
307
308
309
310 PV &Hv_pv = dynamic_cast<PV&>(Hv);
311 const PV &v_pv = dynamic_cast<const PV&>(v);
312
313 // output vector components
314 ROL::Ptr<V> Hvzl = Hv_pv.get(LOWER);
315 ROL::Ptr<V> Hvzu = Hv_pv.get(UPPER);
316
317 // input vector components
318 ROL::Ptr<const V> vzl = v_pv.get(LOWER);
319 ROL::Ptr<const V> vzu = v_pv.get(UPPER);
320
321 Hvzl->set(*x_);
322 Hvzl->axpy(-1.0,*xl_);
323 Hvzl->applyBinary(divinv_,*vzl); // Hvzl = vzl/(x-xl)
324
325 Hvzu->set(*xu_);
326 Hvzu->axpy(-1.0,*x_);
327 Hvzu->applyBinary(divinv_,*vzu); // Hvzu = vzu/(xu-x)
328 }
329
330}; // class PrimalDualInteriorPointBlock22
331
332
333} // namespace ROL
334
335
336
337#endif // ROL_PRIMALDUALINTERIORPOINTOPERATOR_H
338
const Ptr< Obj > obj_
typename PV< Real >::size_type size_type
Provides the interface to apply upper and lower bound constraints.
Defines the general constraint operator interface.
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.
ROL::Ptr< const Vector< Real > > get(size_type i) const
PrimalDualInteriorPointBlock11(ROL::Ptr< OBJ > &obj, ROL::Ptr< CON > &con, const V &x, ROL::Ptr< V > &scratch, Real delta=0)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update linear operator.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
virtual void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
void update(const Vector< Real > &z, bool flag=true, int iter=-1)
Update linear operator.
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
PrimalDualInteriorPointBlock22(const ROL::Ptr< BND > &bnd, const Vector< Real > &x)
virtual void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update linear operator.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84