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