ROL
ROL_PoissonInversion.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
49#ifndef USE_HESSVEC
50#define USE_HESSVEC 1
51#endif
52
53#ifndef ROL_POISSONINVERSION_HPP
54#define ROL_POISSONINVERSION_HPP
55
56#include "ROL_StdVector.hpp"
57#include "ROL_TestProblem.hpp"
59
60#include "Teuchos_LAPACK.hpp"
61
62namespace ROL {
63namespace ZOO {
64
67template<class Real>
69
70 typedef std::vector<Real> vector;
71 typedef Vector<Real> V;
73
74 typedef typename vector::size_type uint;
75
76private:
79
80 Real hu_;
81 Real hz_;
82
83 Real alpha_;
84
85 Real eps_;
87
88 ROL::Ptr<const vector> getVector( const V& x ) {
89
90 return dynamic_cast<const SV&>(x).getVector();
91 }
92
93 ROL::Ptr<vector> getVector( V& x ) {
94
95 return dynamic_cast<SV&>(x).getVector();
96 }
97
98public:
99
100 /* CONSTRUCTOR */
101 Objective_PoissonInversion(uint nz = 32, Real alpha = 1.e-4)
102 : nu_(nz-1), nz_(nz), hu_(1./((Real)nu_+1.)), hz_(hu_),
103 alpha_(alpha), eps_(1.e-4), reg_type_(2) {}
104
105 /* REGULARIZATION DEFINITIONS */
106 Real reg_value(const Vector<Real> &z) {
107
108
109 ROL::Ptr<const vector> zp = getVector(z);
110
111 Real val = 0.0;
112 for (uint i = 0; i < nz_; i++) {
113 if ( reg_type_ == 2 ) {
114 val += alpha_/2.0 * hz_ * (*zp)[i]*(*zp)[i];
115 }
116 else if ( reg_type_ == 1 ) {
117 val += alpha_ * hz_ * std::sqrt((*zp)[i]*(*zp)[i] + eps_);
118 }
119 else if ( reg_type_ == 0 ) {
120 if ( i < nz_-1 ) {
121 val += alpha_ * std::sqrt(std::pow((*zp)[i]-(*zp)[i+1],2.0)+eps_);
122 }
123 }
124 }
125 return val;
126 }
127
129
130
131 if ( reg_type_ == 2 ) {
132 g.set(z);
133 g.scale(alpha_*hz_);
134 }
135 else if ( reg_type_ == 1 ) {
136 ROL::Ptr<const vector> zp = getVector(z);
137 ROL::Ptr<vector > gp = getVector(g);
138
139 for (uint i = 0; i < nz_; i++) {
140 (*gp)[i] = alpha_ * hz_ * (*zp)[i]/std::sqrt(std::pow((*zp)[i],2.0)+eps_);
141 }
142 }
143 else if ( reg_type_ == 0 ) {
144 ROL::Ptr<const vector> zp = getVector(z);
145 ROL::Ptr<vector> gp = getVector(g);
146
147 Real diff = 0.0;
148 for (uint i = 0; i < nz_; i++) {
149 if ( i == 0 ) {
150 diff = (*zp)[i]-(*zp)[i+1];
151 (*gp)[i] = alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
152 }
153 else if ( i == nz_-1 ) {
154 diff = (*zp)[i-1]-(*zp)[i];
155 (*gp)[i] = -alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
156 }
157 else {
158 diff = (*zp)[i]-(*zp)[i+1];
159 (*gp)[i] = alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
160 diff = (*zp)[i-1]-(*zp)[i];
161 (*gp)[i] -= alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
162 }
163 }
164 }
165 }
166
167 void reg_hessVec(Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &z) {
168
169
170
171 if ( reg_type_ == 2 ) {
172 hv.set(v);
173 hv.scale(alpha_*hz_);
174 }
175 else if ( reg_type_ == 1 ) {
176 ROL::Ptr<const vector> zp = getVector(z);
177 ROL::Ptr<const vector> vp = getVector(v);
178 ROL::Ptr<vector> hvp = getVector(hv);
179
180 for (uint i = 0; i < nz_; i++) {
181 (*hvp)[i] = alpha_*hz_*(*vp)[i]*eps_/std::pow(std::pow((*zp)[i],2.0)+eps_,3.0/2.0);
182 }
183 }
184 else if ( reg_type_ == 0 ) {
185 ROL::Ptr<const vector> zp = getVector(z);
186 ROL::Ptr<const vector> vp = getVector(v);
187 ROL::Ptr<vector> hvp = getVector(hv);
188
189 Real diff1 = 0.0;
190 Real diff2 = 0.0;
191 for (uint i = 0; i < nz_; i++) {
192 if ( i == 0 ) {
193 diff1 = (*zp)[i]-(*zp)[i+1];
194 diff1 = eps_/std::pow(std::pow(diff1,2.0)+eps_,3.0/2.0);
195 (*hvp)[i] = alpha_* ((*vp)[i]*diff1 - (*vp)[i+1]*diff1);
196 }
197 else if ( i == nz_-1 ) {
198 diff2 = (*zp)[i-1]-(*zp)[i];
199 diff2 = eps_/std::pow(std::pow(diff2,2.0)+eps_,3.0/2.0);
200 (*hvp)[i] = alpha_* (-(*vp)[i-1]*diff2 + (*vp)[i]*diff2);
201 }
202 else {
203 diff1 = (*zp)[i]-(*zp)[i+1];
204 diff1 = eps_/std::pow(std::pow(diff1,2.0)+eps_,3.0/2.0);
205 diff2 = (*zp)[i-1]-(*zp)[i];
206 diff2 = eps_/std::pow(std::pow(diff2,2.0)+eps_,3.0/2.0);
207 (*hvp)[i] = alpha_* (-(*vp)[i-1]*diff2 + (*vp)[i]*(diff1 + diff2) - (*vp)[i+1]*diff1);
208 }
209 }
210 }
211 }
212
213 /* FINITE ELEMENT DEFINTIONS */
214 void apply_mass(Vector<Real> &Mf, const Vector<Real> &f ) {
215
216
217 ROL::Ptr<const vector> fp = getVector(f);
218 ROL::Ptr<vector> Mfp = getVector(Mf);
219
220 for (uint i = 0; i < nu_; i++) {
221 if ( i == 0 ) {
222 (*Mfp)[i] = hu_/6.0*(4.0*(*fp)[i] + (*fp)[i+1]);
223 }
224 else if ( i == nu_-1 ) {
225 (*Mfp)[i] = hu_/6.0*((*fp)[i-1] + 4.0*(*fp)[i]);
226 }
227 else {
228 (*Mfp)[i] = hu_/6.0*((*fp)[i-1] + 4.0*(*fp)[i] + (*fp)[i+1]);
229 }
230 }
231 }
232
234
235
236 ROL::Ptr<const vector> zp = getVector(z);
237 ROL::Ptr<vector> up = getVector(u);
238 ROL::Ptr<vector> bp = getVector(b);
239
240 // Get Diagonal and Off-Diagonal Entries of PDE Jacobian
241 vector d(nu_,1.0);
242 vector o(nu_-1,1.0);
243 for ( uint i = 0; i < nu_; i++ ) {
244 d[i] = (std::exp((*zp)[i]) + std::exp((*zp)[i+1]))/hu_;
245 if ( i < nu_-1 ) {
246 o[i] *= -std::exp((*zp)[i+1])/hu_;
247 }
248 }
249
250 // Solve Tridiagonal System Using LAPACK's SPD Tridiagonal Solver
251 Teuchos::LAPACK<int,Real> lp;
252 int info;
253 int ldb = nu_;
254 int nhrs = 1;
255 lp.PTTRF(nu_,&d[0],&o[0],&info);
256 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&(*bp)[0],ldb,&info);
257 u.set(b);
258 }
259
260 Real evaluate_target(Real x) {
261 return x*(1.0-x);
262 }
263
265 const Vector<Real> &d, const Vector<Real> &u ) {
266
267
268
269 ROL::Ptr<const vector> zp = getVector(z);
270 ROL::Ptr<const vector> up = getVector(u);
271 ROL::Ptr<const vector> dp = getVector(d);
272 ROL::Ptr<vector> Bdp = getVector(Bd);
273
274 for (uint i = 0; i < nu_; i++) {
275 if ( i == 0 ) {
276 (*Bdp)[i] = 1.0/hu_*( std::exp((*zp)[i])*(*up)[i]*(*dp)[i]
277 + std::exp((*zp)[i+1])*((*up)[i]-(*up)[i+1])*(*dp)[i+1] );
278 }
279 else if ( i == nu_-1 ) {
280 (*Bdp)[i] = 1.0/hu_*( std::exp((*zp)[i])*((*up)[i]-(*up)[i-1])*(*dp)[i]
281 + std::exp((*zp)[i+1])*(*up)[i]*(*dp)[i+1] );
282 }
283 else {
284 (*Bdp)[i] = 1.0/hu_*( std::exp((*zp)[i])*((*up)[i]-(*up)[i-1])*(*dp)[i]
285 + std::exp((*zp)[i+1])*((*up)[i]-(*up)[i+1])*(*dp)[i+1] );
286 }
287 }
288 }
289
291 const Vector<Real> &d, const Vector<Real> &u ) {
292
293
294 ROL::Ptr<const vector> zp = getVector(z);
295 ROL::Ptr<const vector> up = getVector(u);
296 ROL::Ptr<const vector> dp = getVector(d);
297 ROL::Ptr<vector> Bdp = getVector(Bd);
298
299 for (uint i = 0; i < nz_; i++) {
300 if ( i == 0 ) {
301 (*Bdp)[i] = std::exp((*zp)[i])/hu_*(*up)[i]*(*dp)[i];
302 }
303 else if ( i == nz_-1 ) {
304 (*Bdp)[i] = std::exp((*zp)[i])/hu_*(*up)[i-1]*(*dp)[i-1];
305 }
306 else {
307 (*Bdp)[i] = std::exp((*zp)[i])/hu_*( ((*up)[i]-(*up)[i-1])*((*dp)[i]-(*dp)[i-1]) );
308 }
309 }
310 }
311
313 const Vector<Real> &d, const Vector<Real> &u ) {
314
315 ROL::Ptr<const vector> zp = getVector(z);
316 ROL::Ptr<const vector> vp = getVector(v);
317 ROL::Ptr<const vector> up = getVector(u);
318 ROL::Ptr<const vector> dp = getVector(d);
319 ROL::Ptr<vector> Bdp = getVector(Bd);
320
321 for (uint i = 0; i < nz_; i++) {
322 if ( i == 0 ) {
323 (*Bdp)[i] = (*vp)[i]*std::exp((*zp)[i])/hu_*(*up)[i]*(*dp)[i];
324 }
325 else if ( i == nz_-1 ) {
326 (*Bdp)[i] = (*vp)[i]*std::exp((*zp)[i])/hu_*(*up)[i-1]*(*dp)[i-1];
327 }
328 else {
329 (*Bdp)[i] = (*vp)[i]*std::exp((*zp)[i])/hu_*( ((*up)[i]-(*up)[i-1])*((*dp)[i]-(*dp)[i-1]) );
330 }
331 }
332 }
333
334 /* STATE AND ADJOINT EQUATION DEFINTIONS */
336
337
338
339
340 Real k1 = 1.0;
341 Real k2 = 2.0;
342 // Right Hand Side
343 ROL::Ptr<vector> bp = ROL::makePtr<vector>(nu_, 0.0);
344 for ( uint i = 0; i < nu_; i++ ) {
345 if ( (Real)(i+1)*hu_ < 0.5 ) {
346 (*bp)[i] = 2.0*k1*hu_;
347 }
348 else if ( std::abs((Real)(i+1)*hu_ - 0.5) < ROL_EPSILON<Real>() ) {
349 (*bp)[i] = (k1+k2)*hu_;
350 }
351 else if ( (Real)(i+1)*hu_ > 0.5 ) {
352 (*bp)[i] = 2.0*k2*hu_;
353 }
354 }
355
356 SV b(bp);
357 // Solve Equation
358 solve_poisson(u,z,b);
359 }
360
362
363
364
365
366 ROL::Ptr<const vector> up = getVector(u);
367 ROL::Ptr<vector> rp = ROL::makePtr<vector>(nu_,0.0);
368 SV res(rp);
369
370 for ( uint i = 0; i < nu_; i++) {
371 (*rp)[i] = -((*up)[i]-evaluate_target((Real)(i+1)*hu_));
372 }
373 StdVector<Real> Mres( ROL::makePtr<std::vector<Real>>(nu_,0.0) );
374 apply_mass(Mres,res);
375 solve_poisson(p,z,Mres);
376 }
377
379 const Vector<Real> &u, const Vector<Real> &z) {
380
381
382 SV b( ROL::makePtr<vector>(nu_,0.0) );
384 solve_poisson(w,z,b);
385 }
386
388 const Vector<Real> &p, const Vector<Real> &u, const Vector<Real> &z) {
389
390
391
392 SV res( ROL::makePtr<vector>(nu_,0.0) );
393 apply_mass(res,w);
394 SV res1( ROL::makePtr<vector>(nu_,0.0) );
396 res.axpy(-1.0,res1);
397 solve_poisson(q,z,res);
398 }
399
400 /* OBJECTIVE FUNCTION DEFINITIONS */
401 Real value( const Vector<Real> &z, Real &tol ) {
402
403
404
405
406 // SOLVE STATE EQUATION
407 ROL::Ptr<vector> up = ROL::makePtr<vector>(nu_,0.0);
408 SV u( up );
409
411
412 // COMPUTE MISFIT
413 ROL::Ptr<vector> rp = ROL::makePtr<vector>(nu_,0.0);
414 SV res( rp );
415
416 for ( uint i = 0; i < nu_; i++) {
417 (*rp)[i] = ((*up)[i]-evaluate_target((Real)(i+1)*hu_));
418 }
419
420 ROL::Ptr<V> Mres = res.clone();
421 apply_mass(*Mres,res);
422 return 0.5*Mres->dot(res) + reg_value(z);
423 }
424
425 void gradient( Vector<Real> &g, const Vector<Real> &z, Real &tol ) {
426
427
428
429 // SOLVE STATE EQUATION
430 SV u( ROL::makePtr<vector>(nu_,0.0) );
432
433 // SOLVE ADJOINT EQUATION
434 SV p( ROL::makePtr<std::vector<Real>>(nu_,0.0) );
436
437 // Apply Transpose of Linearized Control Operator
439
440 // Regularization gradient
441 SV g_reg( ROL::makePtr<vector>(nz_,0.0) );
442 reg_gradient(g_reg,z);
443
444 // Build Gradient
445 g.plus(g_reg);
446 }
447#if USE_HESSVEC
448 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &z, Real &tol ) {
449
450
451
452 // SOLVE STATE EQUATION
453 SV u( ROL::makePtr<vector>(nu_,0.0) );
455
456 // SOLVE ADJOINT EQUATION
457 SV p( ROL::makePtr<vector>(nu_,0.0) );
459
460 // SOLVE STATE SENSITIVITY EQUATION
461 SV w( ROL::makePtr<vector>(nu_,0.0) );
463
464 // SOLVE ADJOINT SENSITIVITY EQUATION
465 SV q( ROL::makePtr<vector>(nu_,0.0) );
467
468 // Apply Transpose of Linearized Control Operator
470
471 // Apply Transpose of Linearized Control Operator
472 SV tmp( ROL::makePtr<vector>(nz_,0.0) );
474 hv.axpy(-1.0,tmp);
475
476 // Apply Transpose of 2nd Derivative of Control Operator
477 tmp.zero();
479 hv.plus(tmp);
480
481 // Regularization hessVec
482 SV hv_reg( ROL::makePtr<vector>(nz_,0.0) );
483 reg_hessVec(hv_reg,v,z);
484
485 // Build hessVec
486 hv.plus(hv_reg);
487 }
488#endif
489
490 void invHessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
491
492
493
494
495 // Cast hv and v vectors to std::vector.
496 ROL::Ptr<vector> hvp = getVector(hv);
497
498 std::vector<Real> vp(*getVector(v));
499
500 int dim = static_cast<int>(vp.size());
501
502
503 // Compute dense Hessian.
504 Teuchos::SerialDenseMatrix<int, Real> H(dim, dim);
506 H = computeDenseHessian<Real>(obj, x);
507
508 // Compute eigenvalues, sort real part.
509 std::vector<vector> eigenvals = computeEigenvalues<Real>(H);
510 std::sort((eigenvals[0]).begin(), (eigenvals[0]).end());
511
512 // Perform 'inertia' correction.
513 Real inertia = (eigenvals[0])[0];
514 Real correction = 0.0;
515 if ( inertia <= 0.0 ) {
516 correction = 2.0*std::abs(inertia);
517 if ( inertia == 0.0 ) {
518 int cnt = 0;
519 while ( eigenvals[0][cnt] == 0.0 ) {
520 cnt++;
521 }
522 correction = 0.5*eigenvals[0][cnt];
523 if ( cnt == dim-1 ) {
524 correction = 1.0;
525 }
526 }
527 for (int i=0; i<dim; i++) {
528 H(i,i) += correction;
529 }
530 }
531
532 // Compute dense inverse Hessian.
533 Teuchos::SerialDenseMatrix<int, Real> invH = computeInverse<Real>(H);
534
535 // Apply dense inverse Hessian.
536 Teuchos::SerialDenseVector<int, Real> hv_teuchos(Teuchos::View, &((*hvp)[0]), dim);
537 const Teuchos::SerialDenseVector<int, Real> v_teuchos(Teuchos::View, &(vp[0]), dim);
538 hv_teuchos.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, invH, v_teuchos, 0.0);
539 }
540
541};
542
543template<class Real>
544class getPoissonInversion : public TestProblem<Real> {
545public:
547
548 Ptr<Objective<Real>> getObjective(void) const {
549 // Problem dimension
550 int n = 128;
551 // Instantiate Objective Function
552 return ROL::makePtr<Objective_PoissonInversion<Real>>(n,1.e-6);
553 }
554
555 Ptr<Vector<Real>> getInitialGuess(void) const {
556 // Problem dimension
557 int n = 128;
558 // Get Initial Guess
559 ROL::Ptr<std::vector<Real> > x0p = ROL::makePtr<std::vector<Real>>(n,0.0);
560 for ( int i = 0; i < n; i++ ) {
561 (*x0p)[i] = 1.5;
562 }
563 return ROL::makePtr<StdVector<Real>>(x0p);
564 }
565
566 Ptr<Vector<Real>> getSolution(const int i = 0) const {
567 // Problem dimension
568 int n = 128;
569 // Get Solution
570 ROL::Ptr<std::vector<Real> > xp = ROL::makePtr<std::vector<Real>>(n,0.0);
571 Real h = 1.0/((Real)n+1), pt = 0.0, k1 = 1.0, k2 = 2.0;
572 for( int i = 0; i < n; i++ ) {
573 pt = (Real)(i+1)*h;
574 if ( pt >= 0.0 && pt < 0.5 ) {
575 (*xp)[i] = std::log(k1);
576 }
577 else if ( pt >= 0.5 && pt < 1.0 ) {
578 (*xp)[i] = std::log(k2);
579 }
580 }
581 return ROL::makePtr<StdVector<Real>>(xp);
582 }
583};
584
585} // End ZOO Namespace
586} // End ROL Namespace
587
588#endif
Contains definitions for helper functions in ROL.
Contains definitions of test objective functions.
Provides the interface to evaluate objective functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
void axpy(const Real alpha, const Vector< Real > &x)
Compute where .
virtual Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
virtual void scale(const Real alpha)=0
Compute where .
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
void solve_adjoint_equation(Vector< Real > &p, const Vector< Real > &u, const Vector< Real > &z)
void solve_poisson(Vector< Real > &u, const Vector< Real > &z, Vector< Real > &b)
void solve_adjoint_sensitivity_equation(Vector< Real > &q, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &u, const Vector< Real > &z)
Real reg_value(const Vector< Real > &z)
void apply_mass(Vector< Real > &Mf, const Vector< Real > &f)
void reg_hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &z)
ROL::Ptr< const vector > getVector(const V &x)
void apply_transposed_linearized_control_operator(Vector< Real > &Bd, const Vector< Real > &z, const Vector< Real > &d, const Vector< Real > &u)
void solve_state_sensitivity_equation(Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z)
void gradient(Vector< Real > &g, const Vector< Real > &z, Real &tol)
Compute gradient.
void solve_state_equation(Vector< Real > &u, const Vector< Real > &z)
void apply_linearized_control_operator(Vector< Real > &Bd, const Vector< Real > &z, const Vector< Real > &d, const Vector< Real > &u)
Objective_PoissonInversion(uint nz=32, Real alpha=1.e-4)
Real value(const Vector< Real > &z, Real &tol)
Compute value.
void reg_gradient(Vector< Real > &g, const Vector< Real > &z)
void apply_transposed_linearized_control_operator_2(Vector< Real > &Bd, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &d, const Vector< Real > &u)
void invHessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply inverse Hessian approximation to vector.
Ptr< Vector< Real > > getInitialGuess(void) const
Ptr< Objective< Real > > getObjective(void) const
Ptr< Vector< Real > > getSolution(const int i=0) const
constexpr auto dim