ROL
poisson-inversion/example_02.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
49#include "ROL_StdVector.hpp"
50#include "ROL_Objective.hpp"
51#include "ROL_Bounds.hpp"
52#include "ROL_Algorithm.hpp"
55#include "ROL_StatusTest.hpp"
56#include "ROL_Types.hpp"
58#include "ROL_Stream.hpp"
59
60
61#include "Teuchos_GlobalMPISession.hpp"
62
63#include <fstream>
64#include <iostream>
65#include <algorithm>
66
67template<class Real>
69
70 typedef std::vector<Real> vector;
73
74 typedef typename vector::size_type uint;
75
76
77private:
80
81 Real hu_;
82 Real hz_;
83
84 Real u0_;
85 Real u1_;
86
87 Real alpha_;
88
90 Teuchos::SerialDenseMatrix<int, Real> H_;
91
92 ROL::Ptr<const vector> getVector( const V& x ) {
93
94 return dynamic_cast<const SV&>(x).getVector();
95 }
96
97 ROL::Ptr<vector> getVector( V& x ) {
98
99 return dynamic_cast<SV&>(x).getVector();
100 }
101
102public:
103
104 /* CONSTRUCTOR */
105 Objective_PoissonInversion(int nz = 32, Real alpha = 1.e-4)
106 : nz_(nz), u0_(1.0), u1_(0.0), alpha_(alpha), useCorrection_(false) {
107 nu_ = nz_-1;
108 hu_ = 1.0/((Real)nu_+1.0);
109 hz_ = hu_;
110 }
111
112 void apply_mass(std::vector<Real> &Mz, const std::vector<Real> &z ) {
113 Mz.resize(nu_,0.0);
114 for (uint i=0; i<nu_; i++) {
115 if ( i == 0 ) {
116 Mz[i] = hu_/6.0*(2.0*z[i] + z[i+1]);
117 }
118 else if ( i == nu_-1 ) {
119 Mz[i] = hu_/6.0*(z[i-1] + 2.0*z[i]);
120 }
121 else {
122 Mz[i] = hu_/6.0*(z[i-1] + 4.0*z[i] + z[i+1]);
123 }
124 }
125 }
126
127 Real evaluate_target(Real x) {
128 return (x <= 0.5) ? 1.0 : 0.0;
129 }
130
131 void apply_linearized_control_operator( std::vector<Real> &Bd, const std::vector<Real> &z,
132 const std::vector<Real> &d, const std::vector<Real> &u,
133 bool addBC = true ) {
134 Bd.clear();
135 Bd.resize(nu_,0.0);
136 for (uint i = 0; i < nu_; i++) {
137 if ( i == 0 ) {
138 Bd[i] = 1.0/hu_*( u[i]*d[i] + (u[i]-u[i+1])*d[i+1] );
139 }
140 else if ( i == nu_-1 ) {
141 Bd[i] = 1.0/hu_*( (u[i]-u[i-1])*d[i] + u[i]*d[i+1] );
142 }
143 else {
144 Bd[i] = 1.0/hu_*( (u[i]-u[i-1])*d[i] + (u[i]-u[i+1])*d[i+1] );
145 }
146 }
147 if ( addBC ) {
148 Bd[ 0] -= u0_*d[ 0]/hu_;
149 Bd[nu_-1] -= u1_*d[nz_-1]/hu_;
150 }
151 }
152
153 void apply_transposed_linearized_control_operator( std::vector<Real> &Bd, const std::vector<Real> &z,
154 const std::vector<Real> &d, const std::vector<Real> &u,
155 bool addBC = true ) {
156 Bd.clear();
157 Bd.resize(nz_,0.0);
158 for (uint i = 0; i < nz_; i++) {
159 if ( i == 0 ) {
160 Bd[i] = 1.0/hu_*u[i]*d[i];
161 }
162 else if ( i == nz_-1 ) {
163 Bd[i] = 1.0/hu_*u[i-1]*d[i-1];
164 }
165 else {
166 Bd[i] = 1.0/hu_*( (u[i]-u[i-1])*(d[i]-d[i-1]) );
167 }
168 }
169 if ( addBC ) {
170 Bd[ 0] -= u0_*d[ 0]/hu_;
171 Bd[nz_-1] -= u1_*d[nu_-1]/hu_;
172 }
173 }
174
175 /* STATE AND ADJOINT EQUATION DEFINTIONS */
176 void solve_state_equation(std::vector<Real> &u, const std::vector<Real> &z) {
177 // Get Diagonal and Off-Diagonal Entries of PDE Jacobian
178 std::vector<Real> d(nu_,1.0);
179 std::vector<Real> o(nu_-1,1.0);
180 for ( uint i = 0; i < nu_; i++ ) {
181 d[i] = (z[i] + z[i+1])/hu_;
182 if ( i < nu_-1 ) {
183 o[i] *= -z[i+1]/hu_;
184 }
185 }
186 // Set right hand side
187 u.clear();
188 u.resize(nu_,0.0);
189 u[ 0] = z[ 0]/hu_ * u0_;
190 u[nu_-1] = z[nz_-1]/hu_ * u1_;
191 // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
192 Teuchos::LAPACK<int,Real> lp;
193 int info;
194 int ldb = nu_;
195 int nhrs = 1;
196 lp.PTTRF(nu_,&d[0],&o[0],&info);
197 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&u[0],ldb,&info);
198 }
199
200 void solve_adjoint_equation(std::vector<Real> &p, const std::vector<Real> &u, const std::vector<Real> &z) {
201 // Get Diagonal and Off-Diagonal Entries of PDE Jacobian
202 vector d(nu_,1.0);
203 vector o(nu_-1,1.0);
204 for ( uint i = 0; i < nu_; i++ ) {
205 d[i] = (z[i] + z[i+1])/hu_;
206 if ( i < nu_-1 ) {
207 o[i] *= -z[i+1]/hu_;
208 }
209 }
210 // Set right hand side
211 vector r(nu_,0.0);
212 for (uint i = 0; i < nu_; i++) {
213 r[i] = -(u[i]-evaluate_target((Real)(i+1)*hu_));
214 }
215 p.clear();
216 p.resize(nu_,0.0);
217 apply_mass(p,r);
218 // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
219 Teuchos::LAPACK<int,Real> lp;
220 int info;
221 int ldb = nu_;
222 int nhrs = 1;
223 lp.PTTRF(nu_,&d[0],&o[0],&info);
224 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&p[0],ldb,&info);
225 }
226
227 void solve_state_sensitivity_equation(std::vector<Real> &w, const std::vector<Real> &v,
228 const std::vector<Real> &u, const std::vector<Real> &z) {
229 // Get Diagonal and Off-Diagonal Entries of PDE Jacobian
230 vector d(nu_,1.0);
231 vector o(nu_-1,1.0);
232 for ( uint i = 0; i < nu_; i++ ) {
233 d[i] = (z[i] + z[i+1])/hu_;
234 if ( i < nu_-1 ) {
235 o[i] *= -z[i+1]/hu_;
236 }
237 }
238 // Set right hand side
239 w.clear();
240 w.resize(nu_,0.0);
242 // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
243 Teuchos::LAPACK<int,Real> lp;
244 int info;
245 int ldb = nu_;
246 int nhrs = 1;
247 lp.PTTRF(nu_,&d[0],&o[0],&info);
248 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&w[0],ldb,&info);
249 }
250
251 void solve_adjoint_sensitivity_equation(std::vector<Real> &q, const std::vector<Real> &w,
252 const std::vector<Real> &v, const std::vector<Real> &p,
253 const std::vector<Real> &u, const std::vector<Real> &z) {
254 // Get Diagonal and Off-Diagonal Entries of PDE Jacobian
255 vector d(nu_,1.0);
256 vector o(nu_-1,1.0);
257 for ( uint i = 0; i < nu_; i++ ) {
258 d[i] = (z[i] + z[i+1])/hu_;
259 if ( i < nu_-1 ) {
260 o[i] *= -z[i+1]/hu_;
261 }
262 }
263 // Set right hand side
264 q.clear();
265 q.resize(nu_,0.0);
266 apply_mass(q,w);
267 std::vector<Real> res(nu_,0.0);
268 apply_linearized_control_operator(res,z,v,p,false);
269 for (uint i = 0; i < nu_; i++) {
270 q[i] -= res[i];
271 }
272 // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
273 Teuchos::LAPACK<int,Real> lp;
274 int info;
275 int ldb = nu_;
276 int nhrs = 1;
277 lp.PTTRF(nu_,&d[0],&o[0],&info);
278 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&q[0],ldb,&info);
279 }
280
281 void update(const ROL::Vector<Real> &z, bool flag, int iter) {
282
283
284
285 if ( flag && useCorrection_ ) {
286 Real tol = std::sqrt(ROL::ROL_EPSILON<Real>());
287 H_.shape(nz_,nz_);
288 ROL::Ptr<V> e = z.clone();
289 ROL::Ptr<V> h = z.clone();
290 for ( uint i = 0; i < nz_; i++ ) {
291 e = z.basis(i);
292 hessVec_true(*h,*e,z,tol);
293 for ( uint j = 0; j < nz_; j++ ) {
294 e = z.basis(j);
295 (H_)(j,i) = e->dot(*h);
296 }
297 }
298 std::vector<vector> eigenvals = ROL::computeEigenvalues<Real>(H_);
299 std::sort((eigenvals[0]).begin(), (eigenvals[0]).end());
300 Real inertia = (eigenvals[0])[0];
301 Real correction = 0.0;
302 if ( inertia <= 0.0 ) {
303 correction = (1.0+std::sqrt(ROL::ROL_EPSILON<Real>()))*std::abs(inertia);
304 if ( inertia == 0.0 ) {
305 uint cnt = 0;
306 while ( eigenvals[0][cnt] == 0.0 ) {
307 cnt++;
308 }
309 correction = std::sqrt(ROL::ROL_EPSILON<Real>())*eigenvals[0][cnt];
310 if ( cnt == nz_-1 ) {
311 correction = 1.0;
312 }
313 }
314 for ( uint i = 0; i < nz_; i++ ) {
315 (H_)(i,i) += correction;
316 }
317 }
318 }
319 }
320
321 /* OBJECTIVE FUNCTION DEFINITIONS */
322 Real value( const ROL::Vector<Real> &z, Real &tol ) {
323
324
325 ROL::Ptr<const vector> zp = getVector(z);
326
327 // SOLVE STATE EQUATION
328 vector u(nu_,0.0);
330
331 // EVALUATE OBJECTIVE
332 Real val = 0.0;
333 for (uint i=0; i<nz_;i++) {
334 val += hz_*alpha_*0.5*(*zp)[i]*(*zp)[i];
335 }
336 Real res = 0.0;
337 Real res1 = 0.0;
338 Real res2 = 0.0;
339 Real res3 = 0.0;
340 for (uint i=0; i<nu_; i++) {
341 if ( i == 0 ) {
342 res1 = u[i]-evaluate_target((Real)(i+1)*hu_);
343 res2 = u[i+1]-evaluate_target((Real)(i+2)*hu_);
344 res = hu_/6.0*(2.0*res1 + res2)*res1;
345 }
346 else if ( i == nu_-1 ) {
347 res1 = u[i-1]-evaluate_target((Real)i*hu_);
348 res2 = u[i]-evaluate_target((Real)(i+1)*hu_);
349 res = hu_/6.0*(res1 + 2.0*res2)*res2;
350 }
351 else {
352 res1 = u[i-1]-evaluate_target((Real)i*hu_);
353 res2 = u[i]-evaluate_target((Real)(i+1)*hu_);
354 res3 = u[i+1]-evaluate_target((Real)(i+2)*hu_);
355 res = hu_/6.0*(res1 + 4.0*res2 + res3)*res2;
356 }
357 val += 0.5*res;
358 }
359 return val;
360 }
361
362 void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &z, Real &tol ) {
363
364
365
366 ROL::Ptr<const vector> zp = getVector(z);
367 ROL::Ptr<vector> gp = getVector(g);
368
369 // SOLVE STATE EQUATION
370 vector u(nu_,0.0);
372
373 // SOLVE ADJOINT EQUATION
374 vector p(nu_,0.0);
375 solve_adjoint_equation(p,u,*zp);
376
377 // Apply Transpose of Linearized Control Operator
379 // Build Gradient
380 for ( uint i = 0; i < nz_; i++ ) {
381 (*gp)[i] += hz_*alpha_*(*zp)[i];
382 }
383 }
384
385 void hessVec( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &z, Real &tol ) {
386 if ( useCorrection_ ) {
387 hessVec_inertia(hv,v,z,tol);
388 }
389 else {
390 hessVec_true(hv,v,z,tol);
391 }
392 }
393
394 void activateInertia(void) {
395 useCorrection_ = true;
396 }
397
398 void deactivateInertia(void) {
399 useCorrection_ = false;
400 }
401
402 void hessVec_true( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &z, Real &tol ) {
403
404
405
406 ROL::Ptr<const vector> vp = getVector(v);
407 ROL::Ptr<const vector> zp = getVector(z);
408 ROL::Ptr<vector> hvp = getVector(hv);
409
410 // SOLVE STATE EQUATION
411 vector u(nu_,0.0);
413
414 // SOLVE ADJOINT EQUATION
415 vector p(nu_,0.0);
416 solve_adjoint_equation(p,u,*zp);
417
418 // SOLVE STATE SENSITIVITY EQUATION
419 vector w(nu_,0.0);
421 // SOLVE ADJOINT SENSITIVITY EQUATION
422 vector q(nu_,0.0);
423 solve_adjoint_sensitivity_equation(q,w,*vp,p,u,*zp);
424
425 // Apply Transpose of Linearized Control Operator
427
428 // Apply Transpose of Linearized Control Operator
429 std::vector<Real> tmp(nz_,0.0);
431 for (uint i=0; i < nz_; i++) {
432 (*hvp)[i] -= tmp[i];
433 }
434 // Regularization hessVec
435 for (uint i=0; i < nz_; i++) {
436 (*hvp)[i] += hz_*alpha_*(*vp)[i];
437 }
438 }
439
440 void hessVec_inertia( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &z, Real &tol ) {
441
442
443 using ROL::constPtrCast;
444
445 ROL::Ptr<vector> hvp = getVector(hv);
446
447
448 ROL::Ptr<vector> vp = ROL::constPtrCast<vector>(getVector(v));
449
450 Teuchos::SerialDenseVector<int, Real> hv_teuchos(Teuchos::View, &((*hvp)[0]), static_cast<int>(nz_));
451 Teuchos::SerialDenseVector<int, Real> v_teuchos(Teuchos::View, &(( *vp)[0]), static_cast<int>(nz_));
452 hv_teuchos.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, H_, v_teuchos, 0.0);
453 }
454
455};
456
457
458
459typedef double RealT;
460
461int main(int argc, char *argv[]) {
462
463 typedef std::vector<RealT> vector;
464 typedef ROL::Vector<RealT> V;
465 typedef ROL::StdVector<RealT> SV;
466
467 typedef typename vector::size_type uint;
468
469
470
471 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
472
473 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
474 int iprint = argc - 1;
475 ROL::Ptr<std::ostream> outStream;
476 ROL::nullstream bhs; // outputs nothing
477 if (iprint > 0)
478 outStream = ROL::makePtrFromRef(std::cout);
479 else
480 outStream = ROL::makePtrFromRef(bhs);
481
482 int errorFlag = 0;
483
484 // *** Example body.
485
486 try {
487
488 uint dim = 128; // Set problem dimension.
489 RealT alpha = 1.e-6;
491
492 // Iteration vector.
493 ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(dim, 0.0);
494 ROL::Ptr<vector> y_ptr = ROL::makePtr<vector>(dim, 0.0);
495
496 // Set initial guess.
497 for (uint i=0; i<dim; i++) {
498 (*x_ptr)[i] = (RealT)rand()/(RealT)RAND_MAX + 1.e2;
499 (*y_ptr)[i] = (RealT)rand()/(RealT)RAND_MAX + 1.e2;
500 }
501
502 SV x(x_ptr);
503 SV y(y_ptr);
504
505 obj.checkGradient(x,y,true);
506 obj.checkHessVec(x,y,true);
507
508 ROL::Ptr<vector> l_ptr = ROL::makePtr<vector>(dim,1.0);
509 ROL::Ptr<vector> u_ptr = ROL::makePtr<vector>(dim,10.0);
510
511 ROL::Ptr<V> lo = ROL::makePtr<SV>(l_ptr);
512 ROL::Ptr<V> up = ROL::makePtr<SV>(u_ptr);
513
514 ROL::Bounds<RealT> icon(lo,up);
515
516 ROL::ParameterList parlist;
517
518 // Krylov parameters.
519 parlist.sublist("General").sublist("Krylov").set("Absolute Tolerance",1.e-8);
520 parlist.sublist("General").sublist("Krylov").set("Relative Tolerance",1.e-4);
521 parlist.sublist("General").sublist("Krylov").set("Iteration Limit",static_cast<int>(dim));
522 // PDAS parameters.
523 parlist.sublist("Step").sublist("Primal Dual Active Set").set("Relative Step Tolerance",1.e-8);
524 parlist.sublist("Step").sublist("Primal Dual Active Set").set("Relative Gradient Tolerance",1.e-6);
525 parlist.sublist("Step").sublist("Primal Dual Active Set").set("Iteration Limit", 10);
526 parlist.sublist("Step").sublist("Primal Dual Active Set").set("Dual Scaling",(alpha>0.0)?alpha:1.e-4);
527 parlist.sublist("General").sublist("Secant").set("Use as Hessian",true);
528 // Status test parameters.
529 parlist.sublist("Status Test").set("Gradient Tolerance",1.e-12);
530 parlist.sublist("Status Test").set("Step Tolerance",1.e-14);
531 parlist.sublist("Status Test").set("Iteration Limit",100);
532
533 ROL::Ptr<ROL::Step<RealT>> step;
534 ROL::Ptr<ROL::StatusTest<RealT>> status;
535
536 // Define algorithm.
537 step = ROL::makePtr<ROL::PrimalDualActiveSetStep<RealT>>(parlist);
538 status = ROL::makePtr<ROL::StatusTest<RealT>>(parlist);
539 ROL::Algorithm<RealT> algo(step,status,false);
540
541 x.zero();
542 obj.deactivateInertia();
543 algo.run(x,obj,icon,true,*outStream);
544
545 // Output control to file.
546 std::ofstream file;
547 file.open("control_PDAS.txt");
548 for ( uint i = 0; i < dim; i++ ) {
549 file << (*x_ptr)[i] << "\n";
550 }
551 file.close();
552
553 // Projected Newtion.
554 // Define step.
555 parlist.sublist("General").sublist("Secant").set("Use as Hessian",false);
556 parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver", "Truncated CG");
557 parlist.sublist("Step").sublist("Trust Region").set("Initial Radius", 1e3);
558 parlist.sublist("Step").sublist("Trust Region").set("Maximum Radius", 1e8);
559 step = ROL::makePtr<ROL::TrustRegionStep<RealT>>(parlist);
560 status = ROL::makePtr<ROL::StatusTest<RealT>>(parlist);
561 ROL::Algorithm<RealT> algo_tr(step,status,false);
562 // Run Algorithm
563 y.zero();
564 obj.deactivateInertia();
565 algo_tr.run(y,obj,icon,true,*outStream);
566
567 std::ofstream file_tr;
568 file_tr.open("control_TR.txt");
569 for ( uint i = 0; i < dim; i++ ) {
570 file_tr << (*y_ptr)[i] << "\n";
571 }
572 file_tr.close();
573
574 std::vector<RealT> u;
575 obj.solve_state_equation(u,*y_ptr);
576 std::ofstream file_u;
577 file_u.open("state.txt");
578 for ( uint i = 0; i < (dim-1); i++ ) {
579 file_u << u[i] << "\n";
580 }
581 file_u.close();
582
583 ROL::Ptr<V> diff = x.clone();
584 diff->set(x);
585 diff->axpy(-1.0,y);
586 RealT error = diff->norm()/std::sqrt((RealT)dim-1.0);
587 std::cout << "\nError between PDAS solution and TR solution is " << error << "\n";
588 errorFlag = ((error > 1.e2*std::sqrt(ROL::ROL_EPSILON<RealT>())) ? 1 : 0);
589
590 }
591 catch (std::logic_error& err) {
592 *outStream << err.what() << "\n";
593 errorFlag = -1000;
594 }; // end try
595
596 if (errorFlag != 0)
597 std::cout << "End Result: TEST FAILED\n";
598 else
599 std::cout << "End Result: TEST PASSED\n";
600
601 return 0;
602
603}
604
Vector< Real > V
Contains definitions for helper functions in ROL.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Contains definitions of custom data types in ROL.
void solve_state_sensitivity_equation(std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
Teuchos::SerialDenseMatrix< int, Real > H_
void solve_state_equation(std::vector< Real > &u, const std::vector< Real > &z)
void solve_adjoint_equation(std::vector< Real > &p, const std::vector< Real > &u, const std::vector< Real > &z)
void update(const ROL::Vector< Real > &z, bool flag, int iter)
Update objective function.
void hessVec(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
void hessVec_inertia(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
ROL::Ptr< const vector > getVector(const V &x)
void apply_transposed_linearized_control_operator(std::vector< Real > &Bd, const std::vector< Real > &z, const std::vector< Real > &d, const std::vector< Real > &u, bool addBC=true)
void hessVec_true(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
Real value(const ROL::Vector< Real > &z, Real &tol)
Compute value.
Objective_PoissonInversion(int nz=32, Real alpha=1.e-4)
void gradient(ROL::Vector< Real > &g, const ROL::Vector< Real > &z, Real &tol)
Compute gradient.
void apply_linearized_control_operator(std::vector< Real > &Bd, const std::vector< Real > &z, const std::vector< Real > &d, const std::vector< Real > &u, bool addBC=true)
void apply_mass(std::vector< Real > &Mz, const std::vector< Real > &z)
void solve_adjoint_sensitivity_equation(std::vector< Real > &q, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &p, const std::vector< Real > &u, const std::vector< Real > &z)
Provides an interface to run optimization algorithms.
Provides the elementwise interface to apply upper and lower bound constraints.
Provides the interface to evaluate objective functions.
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
int main(int argc, char *argv[])
constexpr auto dim