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_CONSTRAINT_SIMOPT_H 45 #define ROL_CONSTRAINT_SIMOPT_H 46 47 #include "ROL_Constraint.hpp" 48 #include "ROL_Vector_SimOpt.hpp" 49 #include "ROL_Types.hpp" 50 #include <iostream> 51 52 #include "ROL_NonlinearLeastSquaresObjective_SimOpt.hpp" 53 #include "ROL_Constraint_State.hpp" 54 #include "ROL_Objective_FSsolver.hpp" 55 #include "ROL_Algorithm.hpp" 56 #include "ROL_TrustRegionStep.hpp" 57 #include "ROL_CompositeStep.hpp" 58 #include "ROL_ConstraintStatusTest.hpp" 59 60 /** @ingroup func_group 61 \class ROL::Constraint_SimOpt 62 \brief Defines the constraint operator interface for simulation-based optimization. 63 64 This constraint interface inherits from ROL_Constraint, for the 65 use case when \f$\mathcal{X}=\mathcal{U}\times\mathcal{Z}\f$ where \f$\mathcal{U}\f$ and 66 \f$\mathcal{Z}\f$ are Banach spaces. \f$\mathcal{U}\f$ denotes the "simulation space" 67 and \f$\mathcal{Z}\f$ denotes the "optimization space" (of designs, controls, parameters). 68 The simulation-based constraints are of the form 69 \f[ 70 c(u,z) = 0 \,. 71 \f] 72 The basic operator interface, to be implemented by the user, requires: 73 \li #value -- constraint evaluation. 74 \li #applyJacobian_1 -- action of the partial constraint Jacobian --derivatives are 75 with respect to the first component \f$\mathcal{U}\f$; 76 \li #applyJacobian_2 -- action of the partial constraint Jacobian --derivatives are 77 with respect to the second component \f$\mathcal{Z}\f$; 78 \li #applyAdjointJacobian_1 -- action of the adjoint of the partial constraint Jacobian --derivatives are 79 with respect to the first component \f$\mathcal{U}\f$; 80 \li #applyAdjointJacobian_2 -- action of the adjoint of the partial constraint Jacobian --derivatives are 81 with respect to the second component \f$\mathcal{Z}\f$; 82 83 The user may also overload: 84 \li #applyAdjointHessian_11 -- action of the adjoint of the partial constraint Hessian --derivatives 85 are with respect to the first component only; 86 \li #applyAdjointHessian_12 -- action of the adjoint of the partial constraint Hessian --derivatives 87 are with respect to the first and second components; 88 \li #applyAdjointHessian_21 -- action of the adjoint of the partial constraint Hessian --derivatives 89 are with respect to the second and first components; 90 \li #applyAdjointHessian_22 -- action of the adjoint of the partial constraint Hessian --derivatives 91 are with respect to the second component only; 92 \li #solveAugmentedSystem -- solution of the augmented system --the default is an iterative 93 scheme based on the action of the Jacobian and its adjoint. 94 \li #applyPreconditioner -- action of a constraint preconditioner --the default is null-op. 95 96 --- 97 */ 98 99 100 namespace ROL { 101 102 template <class Real> 103 class Constraint_SimOpt : public Constraint<Real> { 104 private: 105 // Additional vector storage for solve 106 Ptr<Vector<Real>> unew_; 107 Ptr<Vector<Real>> jv_; 108 109 // Default parameters for solve (backtracking Newton) 110 const Real DEFAULT_atol_; 111 const Real DEFAULT_rtol_; 112 const Real DEFAULT_stol_; 113 const Real DEFAULT_factor_; 114 const Real DEFAULT_decr_; 115 const int DEFAULT_maxit_; 116 const bool DEFAULT_print_; 117 const bool DEFAULT_zero_; 118 const int DEFAULT_solverType_; 119 120 // User-set parameters for solve (backtracking Newton) 121 122 protected: 123 Real atol_; 124 Real rtol_; 125 Real stol_; 126 Real factor_; 127 Real decr_; 128 int maxit_; 129 bool print_; 130 bool zero_; 131 int solverType_; 132 133 // Flag to initialize vector storage in solve 134 bool firstSolve_; 135 136 public: Constraint_SimOpt()137 Constraint_SimOpt() 138 : Constraint<Real>(), 139 unew_(nullPtr), jv_(nullPtr), 140 DEFAULT_atol_(1.e-4*std::sqrt(ROL_EPSILON<Real>())), 141 DEFAULT_rtol_(1.e0), 142 DEFAULT_stol_(std::sqrt(ROL_EPSILON<Real>())), 143 DEFAULT_factor_(0.5), 144 DEFAULT_decr_(1.e-4), 145 DEFAULT_maxit_(500), 146 DEFAULT_print_(false), 147 DEFAULT_zero_(false), 148 DEFAULT_solverType_(0), 149 atol_(DEFAULT_atol_), rtol_(DEFAULT_rtol_), stol_(DEFAULT_stol_), factor_(DEFAULT_factor_), 150 decr_(DEFAULT_decr_), maxit_(DEFAULT_maxit_), print_(DEFAULT_print_), zero_(DEFAULT_zero_), 151 solverType_(DEFAULT_solverType_), firstSolve_(true) {} 152 153 /** \brief Update constraint functions. 154 x is the optimization variable, 155 flag = true if optimization variable is changed, 156 iter is the outer algorithm iterations count. 157 */ update(const Vector<Real> & u,const Vector<Real> & z,bool flag=true,int iter=-1)158 virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) { 159 update_1(u,flag,iter); 160 update_2(z,flag,iter); 161 } 162 163 /** \brief Update constraint functions with respect to Sim variable. 164 x is the optimization variable, 165 flag = true if optimization variable is changed, 166 iter is the outer algorithm iterations count. 167 */ update_1(const Vector<Real> & u,bool flag=true,int iter=-1)168 virtual void update_1( const Vector<Real> &u, bool flag = true, int iter = -1 ) {} 169 170 /** \brief Update constraint functions with respect to Opt variable. 171 x is the optimization variable, 172 flag = true if optimization variable is changed, 173 iter is the outer algorithm iterations count. 174 */ update_2(const Vector<Real> & z,bool flag=true,int iter=-1)175 virtual void update_2( const Vector<Real> &z, bool flag = true, int iter = -1 ) {} 176 177 178 /** \brief Evaluate the constraint operator \f$c:\mathcal{U}\times\mathcal{Z} \rightarrow \mathcal{C}\f$ 179 at \f$(u,z)\f$. 180 181 @param[out] c is the result of evaluating the constraint operator at @b \f$(u,z)\f$; a constraint-space vector 182 @param[in] u is the constraint argument; a simulation-space vector 183 @param[in] z is the constraint argument; an optimization-space vector 184 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 185 186 On return, \f$\mathsf{c} = c(u,z)\f$, 187 where \f$\mathsf{c} \in \mathcal{C}\f$, \f$\mathsf{u} \in \mathcal{U}\f$, and $\f$\mathsf{z} \in\mathcal{Z}\f$. 188 189 --- 190 */ 191 virtual void value(Vector<Real> &c, 192 const Vector<Real> &u, 193 const Vector<Real> &z, 194 Real &tol) = 0; 195 196 /** \brief Given \f$z\f$, solve \f$c(u,z)=0\f$ for \f$u\f$. 197 198 @param[out] c is the result of evaluating the constraint operator at @b \f$(u,z)\f$; a constraint-space vector 199 @param[in,out] u is the solution vector; a simulation-space vector 200 @param[in] z is the constraint argument; an optimization-space vector 201 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 202 203 The defualt implementation is Newton's method globalized with a backtracking line search. 204 205 --- 206 */ solve(Vector<Real> & c,Vector<Real> & u,const Vector<Real> & z,Real & tol)207 virtual void solve(Vector<Real> &c, 208 Vector<Real> &u, 209 const Vector<Real> &z, 210 Real &tol) { 211 if ( zero_ ) { 212 u.zero(); 213 } 214 update(u,z,true); 215 value(c,u,z,tol); 216 Real cnorm = c.norm(); 217 const Real ctol = std::min(atol_, rtol_*cnorm); 218 if (solverType_==0 || solverType_==3 || solverType_==4) { 219 if ( firstSolve_ ) { 220 unew_ = u.clone(); 221 jv_ = u.clone(); 222 firstSolve_ = false; 223 } 224 Real alpha(1), tmp(0); 225 int cnt = 0; 226 if ( print_ ) { 227 std::cout << "\n Default Constraint_SimOpt::solve\n"; 228 std::cout << " "; 229 std::cout << std::setw(6) << std::left << "iter"; 230 std::cout << std::setw(15) << std::left << "rnorm"; 231 std::cout << std::setw(15) << std::left << "alpha"; 232 std::cout << "\n"; 233 } 234 for (cnt = 0; cnt < maxit_; ++cnt) { 235 // Compute Newton step 236 applyInverseJacobian_1(*jv_,c,u,z,tol); 237 unew_->set(u); 238 unew_->axpy(-alpha, *jv_); 239 update_1(*unew_); 240 //update(*unew_,z); 241 value(c,*unew_,z,tol); 242 tmp = c.norm(); 243 // Perform backtracking line search 244 while ( tmp > (1.0-decr_*alpha)*cnorm && 245 alpha > stol_ ) { 246 alpha *= factor_; 247 unew_->set(u); 248 unew_->axpy(-alpha,*jv_); 249 update_1(*unew_); 250 //update(*unew_,z); 251 value(c,*unew_,z,tol); 252 tmp = c.norm(); 253 } 254 if ( print_ ) { 255 std::cout << " "; 256 std::cout << std::setw(6) << std::left << cnt; 257 std::cout << std::scientific << std::setprecision(6); 258 std::cout << std::setw(15) << std::left << tmp; 259 std::cout << std::scientific << std::setprecision(6); 260 std::cout << std::setw(15) << std::left << alpha; 261 std::cout << "\n"; 262 } 263 // Update iterate 264 cnorm = tmp; 265 u.set(*unew_); 266 if (cnorm <= ctol) { // = covers the case of identically zero residual 267 break; 268 } 269 update(u,z,true); 270 alpha = 1.0; 271 } 272 } 273 if (solverType_==1 || (solverType_==3 && cnorm > ctol)) { 274 Ptr<Constraint_SimOpt<Real>> con = makePtrFromRef(*this); 275 Ptr<Objective<Real>> obj = makePtr<NonlinearLeastSquaresObjective_SimOpt<Real>>(con,u,z,c,true); 276 ParameterList parlist; 277 parlist.sublist("Status Test").set("Gradient Tolerance",ctol); 278 parlist.sublist("Status Test").set("Step Tolerance",stol_); 279 parlist.sublist("Status Test").set("Iteration Limit",maxit_); 280 parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG"); 281 parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100); 282 Ptr<Step<Real>> step = makePtr<TrustRegionStep<Real>>(parlist); 283 Ptr<StatusTest<Real>> status = makePtr<StatusTest<Real>>(parlist); 284 Ptr<Algorithm<Real>> algo = makePtr<Algorithm<Real>>(step,status,false); 285 algo->run(u,*obj,print_); 286 value(c,u,z,tol); 287 } 288 if (solverType_==2 || (solverType_==4 && cnorm > ctol)) { 289 Ptr<Constraint_SimOpt<Real>> con = makePtrFromRef(*this); 290 Ptr<const Vector<Real>> zVec = makePtrFromRef(z); 291 Ptr<Constraint<Real>> conU 292 = makePtr<Constraint_State<Real>>(con,zVec); 293 Ptr<Objective<Real>> objU 294 = makePtr<Objective_FSsolver<Real>>(); 295 ParameterList parlist; 296 parlist.sublist("Status Test").set("Constraint Tolerance",ctol); 297 parlist.sublist("Status Test").set("Step Tolerance",stol_); 298 parlist.sublist("Status Test").set("Iteration Limit",maxit_); 299 Ptr<Step<Real>> step = makePtr<CompositeStep<Real>>(parlist); 300 Ptr<StatusTest<Real>> status = makePtr<ConstraintStatusTest<Real>>(parlist); 301 Ptr<Algorithm<Real>> algo = makePtr<Algorithm<Real>>(step,status,false); 302 Ptr<Vector<Real>> l = c.dual().clone(); 303 algo->run(u,*l,*objU,*conU,print_); 304 value(c,u,z,tol); 305 } 306 if (solverType_ > 4 || solverType_ < 0) { 307 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument, 308 ">>> ERROR (ROL:Constraint_SimOpt:solve): Invalid solver type!"); 309 } 310 } 311 312 /** \brief Set solve parameters. 313 314 @param[in] parlist ParameterList containing solve parameters 315 316 For the default implementation, parlist has two sublist ("SimOpt" 317 and "Solve") and the "Solve" sublist has six input parameters. 318 319 - "Residual Tolerance": Absolute tolerance for the norm of the residual (Real) 320 - "Iteration Limit": Maximum number of Newton iterations (int) 321 - "Sufficient Decrease Tolerance": Tolerance signifying sufficient decrease in the residual norm, between 0 and 1 (Real) 322 - "Step Tolerance": Absolute tolerance for the step size parameter (Real) 323 - "Backtracking Factor": Rate for decreasing step size during backtracking, between 0 and 1 (Real) 324 - "Output Iteration History": Set to true in order to print solve iteration history (bool) 325 - "Zero Initial Guess": Use a vector of zeros as an initial guess for the solve (bool) 326 - "Solver Type": Determine which solver to use (0: Newton with line search, 1: Levenberg-Marquardt, 2: SQP) (int) 327 328 These parameters are accessed as parlist.sublist("SimOpt").sublist("Solve").get(...). 329 330 --- 331 */ setSolveParameters(ParameterList & parlist)332 virtual void setSolveParameters(ParameterList &parlist) { 333 ParameterList & list = parlist.sublist("SimOpt").sublist("Solve"); 334 atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_); 335 rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_); 336 maxit_ = list.get("Iteration Limit", DEFAULT_maxit_); 337 decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_); 338 stol_ = list.get("Step Tolerance", DEFAULT_stol_); 339 factor_ = list.get("Backtracking Factor", DEFAULT_factor_); 340 print_ = list.get("Output Iteration History", DEFAULT_print_); 341 zero_ = list.get("Zero Initial Guess", DEFAULT_zero_); 342 solverType_ = list.get("Solver Type", DEFAULT_solverType_); 343 } 344 345 /** \brief Apply the partial constraint Jacobian at \f$(u,z)\f$, 346 \f$c_u(u,z) \in L(\mathcal{U}, \mathcal{C})\f$, 347 to the vector \f$v\f$. 348 349 @param[out] jv is the result of applying the constraint Jacobian to @b v at @b \f$(u,z)\f$; a constraint-space vector 350 @param[in] v is a simulation-space vector 351 @param[in] u is the constraint argument; an simulation-space vector 352 @param[in] z is the constraint argument; an optimization-space vector 353 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 354 355 On return, \f$\mathsf{jv} = c_u(u,z)v\f$, where 356 \f$v \in \mathcal{U}\f$, \f$\mathsf{jv} \in \mathcal{C}\f$. 357 358 --- 359 */ applyJacobian_1(Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)360 virtual void applyJacobian_1(Vector<Real> &jv, 361 const Vector<Real> &v, 362 const Vector<Real> &u, 363 const Vector<Real> &z, 364 Real &tol) { 365 Real ctol = std::sqrt(ROL_EPSILON<Real>()); 366 // Compute step length 367 Real h = tol; 368 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 369 h = std::max(1.0,u.norm()/v.norm())*tol; 370 } 371 // Update state vector to u + hv 372 Ptr<Vector<Real>> unew = u.clone(); 373 unew->set(u); 374 unew->axpy(h,v); 375 // Compute new constraint value 376 update(*unew,z); 377 value(jv,*unew,z,ctol); 378 // Compute current constraint value 379 Ptr<Vector<Real>> cold = jv.clone(); 380 update(u,z); 381 value(*cold,u,z,ctol); 382 // Compute Newton quotient 383 jv.axpy(-1.0,*cold); 384 jv.scale(1.0/h); 385 } 386 387 388 /** \brief Apply the partial constraint Jacobian at \f$(u,z)\f$, 389 \f$c_z(u,z) \in L(\mathcal{Z}, \mathcal{C})\f$, 390 to the vector \f$v\f$. 391 392 @param[out] jv is the result of applying the constraint Jacobian to @b v at @b \f$(u,z)\f$; a constraint-space vector 393 @param[in] v is an optimization-space vector 394 @param[in] u is the constraint argument; a simulation-space vector 395 @param[in] z is the constraint argument; an optimization-space vector 396 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 397 398 On return, \f$\mathsf{jv} = c_z(u,z)v\f$, where 399 \f$v \in \mathcal{Z}\f$, \f$\mathsf{jv} \in \mathcal{C}\f$. 400 401 --- 402 */ applyJacobian_2(Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)403 virtual void applyJacobian_2(Vector<Real> &jv, 404 const Vector<Real> &v, 405 const Vector<Real> &u, 406 const Vector<Real> &z, 407 Real &tol) { 408 Real ctol = std::sqrt(ROL_EPSILON<Real>()); 409 // Compute step length 410 Real h = tol; 411 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 412 h = std::max(1.0,u.norm()/v.norm())*tol; 413 } 414 // Update state vector to u + hv 415 Ptr<Vector<Real>> znew = z.clone(); 416 znew->set(z); 417 znew->axpy(h,v); 418 // Compute new constraint value 419 update(u,*znew); 420 value(jv,u,*znew,ctol); 421 // Compute current constraint value 422 Ptr<Vector<Real>> cold = jv.clone(); 423 update(u,z); 424 value(*cold,u,z,ctol); 425 // Compute Newton quotient 426 jv.axpy(-1.0,*cold); 427 jv.scale(1.0/h); 428 } 429 430 /** \brief Apply the inverse partial constraint Jacobian at \f$(u,z)\f$, 431 \f$c_u(u,z)^{-1} \in L(\mathcal{C}, \mathcal{U})\f$, 432 to the vector \f$v\f$. 433 434 @param[out] ijv is the result of applying the inverse constraint Jacobian to @b v at @b \f$(u,z)\f$; a simulation-space vector 435 @param[in] v is a constraint-space vector 436 @param[in] u is the constraint argument; a simulation-space vector 437 @param[in] z is the constraint argument; an optimization-space vector 438 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 439 440 On return, \f$\mathsf{ijv} = c_u(u,z)^{-1}v\f$, where 441 \f$v \in \mathcal{C}\f$, \f$\mathsf{ijv} \in \mathcal{U}\f$. 442 443 --- 444 */ applyInverseJacobian_1(Vector<Real> & ijv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)445 virtual void applyInverseJacobian_1(Vector<Real> &ijv, 446 const Vector<Real> &v, 447 const Vector<Real> &u, 448 const Vector<Real> &z, 449 Real &tol) { 450 ROL_TEST_FOR_EXCEPTION(true, std::logic_error, 451 "The method applyInverseJacobian_1 is used but not implemented!\n"); 452 } 453 454 /** \brief Apply the adjoint of the partial constraint Jacobian at \f$(u,z)\f$, 455 \f$c_u(u,z)^* \in L(\mathcal{C}^*, \mathcal{U}^*)\f$, 456 to the vector \f$v\f$. This is the primary interface. 457 458 @param[out] ajv is the result of applying the adjoint of the constraint Jacobian to @b v at @b (u,z); a dual simulation-space vector 459 @param[in] v is a dual constraint-space vector 460 @param[in] u is the constraint argument; a simulation-space vector 461 @param[in] z is the constraint argument; an optimization-space vector 462 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 463 464 On return, \f$\mathsf{ajv} = c_u(u,z)^*v\f$, where 465 \f$v \in \mathcal{C}^*\f$, \f$\mathsf{ajv} \in \mathcal{U}^*\f$. 466 467 --- 468 */ applyAdjointJacobian_1(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)469 virtual void applyAdjointJacobian_1(Vector<Real> &ajv, 470 const Vector<Real> &v, 471 const Vector<Real> &u, 472 const Vector<Real> &z, 473 Real &tol) { 474 applyAdjointJacobian_1(ajv, v, u, z, v.dual(), tol); 475 } 476 477 478 /** \brief Apply the adjoint of the partial constraint Jacobian at \f$(u,z)\f$, 479 \f$c_u(u,z)^* \in L(\mathcal{C}^*, \mathcal{U}^*)\f$, 480 to the vector \f$v\f$. This is the secondary interface, for use 481 with dual spaces where the user does not define the dual() operation. 482 483 @param[out] ajv is the result of applying the adjoint of the constraint Jacobian to @b v at @b (u,z); a dual simulation-space vector 484 @param[in] v is a dual constraint-space vector 485 @param[in] u is the constraint argument; a simulation-space vector 486 @param[in] z is the constraint argument; an optimization-space vector 487 @param[in] dualv is a vector used for temporary variables; a constraint-space vector 488 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 489 490 On return, \f$\mathsf{ajv} = c_u(u,z)^*v\f$, where 491 \f$v \in \mathcal{C}^*\f$, \f$\mathsf{ajv} \in \mathcal{U}^*\f$. 492 493 --- 494 */ applyAdjointJacobian_1(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & dualv,Real & tol)495 virtual void applyAdjointJacobian_1(Vector<Real> &ajv, 496 const Vector<Real> &v, 497 const Vector<Real> &u, 498 const Vector<Real> &z, 499 const Vector<Real> &dualv, 500 Real &tol) { 501 Real ctol = std::sqrt(ROL_EPSILON<Real>()); 502 Real h = tol; 503 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 504 h = std::max(1.0,u.norm()/v.norm())*tol; 505 } 506 Ptr<Vector<Real>> cold = dualv.clone(); 507 Ptr<Vector<Real>> cnew = dualv.clone(); 508 update(u,z); 509 value(*cold,u,z,ctol); 510 Ptr<Vector<Real>> unew = u.clone(); 511 ajv.zero(); 512 for (int i = 0; i < u.dimension(); i++) { 513 unew->set(u); 514 unew->axpy(h,*(u.basis(i))); 515 update(*unew,z); 516 value(*cnew,*unew,z,ctol); 517 cnew->axpy(-1.0,*cold); 518 cnew->scale(1.0/h); 519 ajv.axpy(cnew->dot(v),*((u.dual()).basis(i))); 520 } 521 update(u,z); 522 } 523 524 525 /** \brief Apply the adjoint of the partial constraint Jacobian at \f$(u,z)\f$, 526 \f$c_z(u,z)^* \in L(\mathcal{C}^*, \mathcal{Z}^*)\f$, 527 to vector \f$v\f$. This is the primary interface. 528 529 @param[out] ajv is the result of applying the adjoint of the constraint Jacobian to @b v at @b \f$(u,z)\f$; a dual optimization-space vector 530 @param[in] v is a dual constraint-space vector 531 @param[in] u is the constraint argument; a simulation-space vector 532 @param[in] z is the constraint argument; an optimization-space vector 533 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 534 535 On return, \f$\mathsf{ajv} = c_z(u,z)^*v\f$, where 536 \f$v \in \mathcal{C}^*\f$, \f$\mathsf{ajv} \in \mathcal{Z}^*\f$. 537 538 --- 539 */ applyAdjointJacobian_2(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)540 virtual void applyAdjointJacobian_2(Vector<Real> &ajv, 541 const Vector<Real> &v, 542 const Vector<Real> &u, 543 const Vector<Real> &z, 544 Real &tol) { 545 applyAdjointJacobian_2(ajv, v, u, z, v.dual(), tol); 546 } 547 548 549 /** \brief Apply the adjoint of the partial constraint Jacobian at \f$(u,z)\f$, 550 \f$c_z(u,z)^* \in L(\mathcal{C}^*, \mathcal{Z}^*)\f$, 551 to vector \f$v\f$. This is the secondary interface, for use 552 with dual spaces where the user does not define the dual() operation. 553 554 @param[out] ajv is the result of applying the adjoint of the constraint Jacobian to @b v at @b \f$(u,z)\f$; a dual optimization-space vector 555 @param[in] v is a dual constraint-space vector 556 @param[in] u is the constraint argument; a simulation-space vector 557 @param[in] z is the constraint argument; an optimization-space vector 558 @param[in] dualv is a vector used for temporary variables; a constraint-space vector 559 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 560 561 On return, \f$\mathsf{ajv} = c_z(u,z)^*v\f$, where 562 \f$v \in \mathcal{C}^*\f$, \f$\mathsf{ajv} \in \mathcal{Z}^*\f$. 563 564 --- 565 */ applyAdjointJacobian_2(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & dualv,Real & tol)566 virtual void applyAdjointJacobian_2(Vector<Real> &ajv, 567 const Vector<Real> &v, 568 const Vector<Real> &u, 569 const Vector<Real> &z, 570 const Vector<Real> &dualv, 571 Real &tol) { 572 Real ctol = std::sqrt(ROL_EPSILON<Real>()); 573 Real h = tol; 574 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 575 h = std::max(1.0,u.norm()/v.norm())*tol; 576 } 577 Ptr<Vector<Real>> cold = dualv.clone(); 578 Ptr<Vector<Real>> cnew = dualv.clone(); 579 update(u,z); 580 value(*cold,u,z,ctol); 581 Ptr<Vector<Real>> znew = z.clone(); 582 ajv.zero(); 583 for (int i = 0; i < z.dimension(); i++) { 584 znew->set(z); 585 znew->axpy(h,*(z.basis(i))); 586 update(u,*znew); 587 value(*cnew,u,*znew,ctol); 588 cnew->axpy(-1.0,*cold); 589 cnew->scale(1.0/h); 590 ajv.axpy(cnew->dot(v),*((z.dual()).basis(i))); 591 } 592 update(u,z); 593 } 594 595 /** \brief Apply the inverse of the adjoint of the partial constraint Jacobian at \f$(u,z)\f$, 596 \f$c_u(u,z)^{-*} \in L(\mathcal{U}^*, \mathcal{C}^*)\f$, 597 to the vector \f$v\f$. 598 599 @param[out] iajv is the result of applying the inverse adjoint of the constraint Jacobian to @b v at @b (u,z); a dual constraint-space vector 600 @param[in] v is a dual simulation-space vector 601 @param[in] u is the constraint argument; a simulation-space vector 602 @param[in] z is the constraint argument; an optimization-space vector 603 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 604 605 On return, \f$\mathsf{iajv} = c_u(u,z)^{-*}v\f$, where 606 \f$v \in \mathcal{U}^*\f$, \f$\mathsf{iajv} \in \mathcal{C}^*\f$. 607 608 --- 609 */ applyInverseAdjointJacobian_1(Vector<Real> & iajv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)610 virtual void applyInverseAdjointJacobian_1(Vector<Real> &iajv, 611 const Vector<Real> &v, 612 const Vector<Real> &u, 613 const Vector<Real> &z, 614 Real &tol) { 615 ROL_TEST_FOR_EXCEPTION(true, std::logic_error, 616 "The method applyInverseAdjointJacobian_1 is used but not implemented!\n"); 617 }; 618 619 /** \brief Apply the simulation-space derivative of the adjoint of the constraint 620 simulation-space Jacobian at \f$(u,z)\f$ to the vector \f$w\f$ in the 621 direction \f$v\f$, according to \f$v\mapsto c_{uu}(u,z)(v,\cdot)^*w\f$. 622 623 @param[out] ahwv is the result of applying the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at @b \f$(u,z)\f$ to the vector @b \f$w\f$ in direction @b \f$w\f$; a dual simulation-space vector 624 @param[in] w is the direction vector; a dual constraint-space vector 625 @param[in] v is a simulation-space vector 626 @param[in] u is the constraint argument; a simulation-space vector 627 @param[in] z is the constraint argument; an optimization-space vector 628 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 629 630 On return, \f$\mathsf{ahwv} = c_{uu}(u,z)(v,\cdot)^*w\f$, where 631 \f$w \in \mathcal{C}^*\f$, \f$v \in \mathcal{U}\f$, and 632 \f$\mathsf{ahwv} \in \mathcal{U}^*\f$. 633 634 --- 635 */ applyAdjointHessian_11(Vector<Real> & ahwv,const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)636 virtual void applyAdjointHessian_11(Vector<Real> &ahwv, 637 const Vector<Real> &w, 638 const Vector<Real> &v, 639 const Vector<Real> &u, 640 const Vector<Real> &z, 641 Real &tol) { 642 Real jtol = std::sqrt(ROL_EPSILON<Real>()); 643 // Compute step size 644 Real h = tol; 645 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 646 h = std::max(1.0,u.norm()/v.norm())*tol; 647 } 648 // Evaluate Jacobian at new state 649 Ptr<Vector<Real>> unew = u.clone(); 650 unew->set(u); 651 unew->axpy(h,v); 652 update(*unew,z); 653 applyAdjointJacobian_1(ahwv,w,*unew,z,jtol); 654 // Evaluate Jacobian at old state 655 Ptr<Vector<Real>> jv = ahwv.clone(); 656 update(u,z); 657 applyAdjointJacobian_1(*jv,w,u,z,jtol); 658 // Compute Newton quotient 659 ahwv.axpy(-1.0,*jv); 660 ahwv.scale(1.0/h); 661 } 662 663 664 /** \brief Apply the optimization-space derivative of the adjoint of the constraint 665 simulation-space Jacobian at \f$(u,z)\f$ to the vector \f$w\f$ in the 666 direction \f$v\f$, according to \f$v\mapsto c_{uz}(u,z)(v,\cdot)^*w\f$. 667 668 @param[out] ahwv is the result of applying the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at @b \f$(u,z)\f$ to the vector @b \f$w\f$ in direction @b \f$w\f$; a dual optimization-space vector 669 @param[in] w is the direction vector; a dual constraint-space vector 670 @param[in] v is a simulation-space vector 671 @param[in] u is the constraint argument; a simulation-space vector 672 @param[in] z is the constraint argument; an optimization-space vector 673 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 674 675 On return, \f$\mathsf{ahwv} = c_{uz}(u,z)(v,\cdot)^*w\f$, where 676 \f$w \in \mathcal{C}^*\f$, \f$v \in \mathcal{U}\f$, and 677 \f$\mathsf{ahwv} \in \mathcal{Z}^*\f$. 678 679 --- 680 */ applyAdjointHessian_12(Vector<Real> & ahwv,const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)681 virtual void applyAdjointHessian_12(Vector<Real> &ahwv, 682 const Vector<Real> &w, 683 const Vector<Real> &v, 684 const Vector<Real> &u, 685 const Vector<Real> &z, 686 Real &tol) { 687 Real jtol = std::sqrt(ROL_EPSILON<Real>()); 688 // Compute step size 689 Real h = tol; 690 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 691 h = std::max(1.0,u.norm()/v.norm())*tol; 692 } 693 // Evaluate Jacobian at new state 694 Ptr<Vector<Real>> unew = u.clone(); 695 unew->set(u); 696 unew->axpy(h,v); 697 update(*unew,z); 698 applyAdjointJacobian_2(ahwv,w,*unew,z,jtol); 699 // Evaluate Jacobian at old state 700 Ptr<Vector<Real>> jv = ahwv.clone(); 701 update(u,z); 702 applyAdjointJacobian_2(*jv,w,u,z,jtol); 703 // Compute Newton quotient 704 ahwv.axpy(-1.0,*jv); 705 ahwv.scale(1.0/h); 706 } 707 708 709 /** \brief Apply the simulation-space derivative of the adjoint of the constraint 710 optimization-space Jacobian at \f$(u,z)\f$ to the vector \f$w\f$ in the 711 direction \f$v\f$, according to \f$v\mapsto c_{zu}(u,z)(v,\cdot)^*w\f$. 712 713 @param[out] ahwv is the result of applying the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at @b \f$(u,z)\f$ to the vector @b \f$w\f$ in direction @b \f$w\f$; a dual simulation-space vector 714 @param[in] w is the direction vector; a dual constraint-space vector 715 @param[in] v is a optimization-space vector 716 @param[in] u is the constraint argument; a simulation-space vector 717 @param[in] z is the constraint argument; an optimization-space vector 718 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 719 720 On return, \f$\mathsf{ahwv} = c_{zu}(u,z)(v,\cdot)^*w\f$, where 721 \f$w \in \mathcal{C}^*\f$, \f$v \in \mathcal{Z}\f$, and 722 \f$\mathsf{ahwv} \in \mathcal{U}^*\f$. 723 724 --- 725 */ applyAdjointHessian_21(Vector<Real> & ahwv,const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)726 virtual void applyAdjointHessian_21(Vector<Real> &ahwv, 727 const Vector<Real> &w, 728 const Vector<Real> &v, 729 const Vector<Real> &u, 730 const Vector<Real> &z, 731 Real &tol) { 732 Real jtol = std::sqrt(ROL_EPSILON<Real>()); 733 // Compute step size 734 Real h = tol; 735 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 736 h = std::max(1.0,u.norm()/v.norm())*tol; 737 } 738 // Evaluate Jacobian at new control 739 Ptr<Vector<Real>> znew = z.clone(); 740 znew->set(z); 741 znew->axpy(h,v); 742 update(u,*znew); 743 applyAdjointJacobian_1(ahwv,w,u,*znew,jtol); 744 // Evaluate Jacobian at old control 745 Ptr<Vector<Real>> jv = ahwv.clone(); 746 update(u,z); 747 applyAdjointJacobian_1(*jv,w,u,z,jtol); 748 // Compute Newton quotient 749 ahwv.axpy(-1.0,*jv); 750 ahwv.scale(1.0/h); 751 } 752 753 /** \brief Apply the optimization-space derivative of the adjoint of the constraint 754 optimization-space Jacobian at \f$(u,z)\f$ to the vector \f$w\f$ in the 755 direction \f$v\f$, according to \f$v\mapsto c_{zz}(u,z)(v,\cdot)^*w\f$. 756 757 @param[out] ahwv is the result of applying the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian at @b \f$(u,z)\f$ to the vector @b \f$w\f$ in direction @b \f$w\f$; a dual optimization-space vector 758 @param[in] w is the direction vector; a dual constraint-space vector 759 @param[in] v is a optimization-space vector 760 @param[in] u is the constraint argument; a simulation-space vector 761 @param[in] z is the constraint argument; an optimization-space vector 762 @param[in,out] tol is a tolerance for inexact evaluations; currently unused 763 764 On return, \f$\mathsf{ahwv} = c_{zz}(u,z)(v,\cdot)^*w\f$, where 765 \f$w \in \mathcal{C}^*\f$, \f$v \in \mathcal{Z}\f$, and 766 \f$\mathsf{ahwv} \in \mathcal{Z}^*\f$. 767 768 --- 769 */ applyAdjointHessian_22(Vector<Real> & ahwv,const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,Real & tol)770 virtual void applyAdjointHessian_22(Vector<Real> &ahwv, 771 const Vector<Real> &w, 772 const Vector<Real> &v, 773 const Vector<Real> &u, 774 const Vector<Real> &z, 775 Real &tol) { 776 Real jtol = std::sqrt(ROL_EPSILON<Real>()); 777 // Compute step size 778 Real h = tol; 779 if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) { 780 h = std::max(1.0,u.norm()/v.norm())*tol; 781 } 782 // Evaluate Jacobian at new control 783 Ptr<Vector<Real>> znew = z.clone(); 784 znew->set(z); 785 znew->axpy(h,v); 786 update(u,*znew); 787 applyAdjointJacobian_2(ahwv,w,u,*znew,jtol); 788 // Evaluate Jacobian at old control 789 Ptr<Vector<Real>> jv = ahwv.clone(); 790 update(u,z); 791 applyAdjointJacobian_2(*jv,w,u,z,jtol); 792 // Compute Newton quotient 793 ahwv.axpy(-1.0,*jv); 794 ahwv.scale(1.0/h); 795 } 796 797 /** \brief Approximately solves the <em> augmented system </em> 798 \f[ 799 \begin{pmatrix} 800 I & c'(x)^* \\ 801 c'(x) & 0 802 \end{pmatrix} 803 \begin{pmatrix} 804 v_{1} \\ 805 v_{2} 806 \end{pmatrix} 807 = 808 \begin{pmatrix} 809 b_{1} \\ 810 b_{2} 811 \end{pmatrix} 812 \f] 813 where \f$v_{1} \in \mathcal{X}\f$, \f$v_{2} \in \mathcal{C}^*\f$, 814 \f$b_{1} \in \mathcal{X}^*\f$, \f$b_{2} \in \mathcal{C}\f$, 815 \f$I : \mathcal{X} \rightarrow \mathcal{X}^*\f$ is an identity 816 operator, and \f$0 : \mathcal{C}^* \rightarrow \mathcal{C}\f$ 817 is a zero operator. 818 819 @param[out] v1 is the optimization-space component of the result 820 @param[out] v2 is the dual constraint-space component of the result 821 @param[in] b1 is the dual optimization-space component of the right-hand side 822 @param[in] b2 is the constraint-space component of the right-hand side 823 @param[in] x is the constraint argument; an optimization-space vector 824 @param[in,out] tol is the nominal relative residual tolerance 825 826 On return, \f$ [\mathsf{v1} \,\, \mathsf{v2}] \f$ approximately 827 solves the augmented system, where the size of the residual is 828 governed by special stopping conditions. \n\n 829 The default implementation is the preconditioned generalized 830 minimal residual (GMRES) method, which enables the use of 831 nonsymmetric preconditioners. 832 833 --- 834 */ solveAugmentedSystem(Vector<Real> & v1,Vector<Real> & v2,const Vector<Real> & b1,const Vector<Real> & b2,const Vector<Real> & x,Real & tol)835 virtual std::vector<Real> solveAugmentedSystem(Vector<Real> &v1, 836 Vector<Real> &v2, 837 const Vector<Real> &b1, 838 const Vector<Real> &b2, 839 const Vector<Real> &x, 840 Real &tol) { 841 return Constraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,x,tol); 842 } 843 844 845 /** \brief Apply a constraint preconditioner at \f$x\f$, \f$P(x) \in L(\mathcal{C}, \mathcal{C})\f$, 846 to vector \f$v\f$. In general, this preconditioner satisfies the following relationship: 847 \f[ 848 c'(x) c'(x)^* P(x) v \approx v \,. 849 \f] 850 It is used by the #solveAugmentedSystem method. 851 852 @param[out] pv is the result of applying the constraint preconditioner to @b v at @b x; a constraint-space vector 853 @param[in] v is a constraint-space vector 854 @param[in] x is the preconditioner argument; an optimization-space vector 855 @param[in,out] tol is a tolerance for inexact evaluations 856 857 On return, \f$\mathsf{pv} = P(x)v\f$, where 858 \f$v \in \mathcal{C}\f$, \f$\mathsf{pv} \in \mathcal{C}\f$. \n\n 859 The default implementation is a null-op. 860 861 --- 862 */ applyPreconditioner(Vector<Real> & pv,const Vector<Real> & v,const Vector<Real> & x,const Vector<Real> & g,Real & tol)863 virtual void applyPreconditioner(Vector<Real> &pv, 864 const Vector<Real> &v, 865 const Vector<Real> &x, 866 const Vector<Real> &g, 867 Real &tol) { 868 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>(x); 869 Ptr<Vector<Real>> ijv = (xs.get_1())->clone(); 870 871 try { 872 applyInverseJacobian_1(*ijv, v, *(xs.get_1()), *(xs.get_2()), tol); 873 } 874 catch (const std::logic_error &e) { 875 Constraint<Real>::applyPreconditioner(pv, v, x, g, tol); 876 return; 877 } 878 879 const Vector_SimOpt<Real> &gs = dynamic_cast<const Vector_SimOpt<Real>&>(g); 880 Ptr<Vector<Real>> ijv_dual = (gs.get_1())->clone(); 881 ijv_dual->set(ijv->dual()); 882 883 try { 884 applyInverseAdjointJacobian_1(pv, *ijv_dual, *(xs.get_1()), *(xs.get_2()), tol); 885 } 886 catch (const std::logic_error &e) { 887 Constraint<Real>::applyPreconditioner(pv, v, x, g, tol); 888 return; 889 } 890 891 } 892 893 /** \brief Update constraint functions. 894 x is the optimization variable, 895 flag = true if optimization variable is changed, 896 iter is the outer algorithm iterations count. 897 */ update(const Vector<Real> & x,bool flag=true,int iter=-1)898 virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) { 899 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>( 900 dynamic_cast<const Vector<Real>&>(x)); 901 update(*(xs.get_1()),*(xs.get_2()),flag,iter); 902 } 903 value(Vector<Real> & c,const Vector<Real> & x,Real & tol)904 virtual void value(Vector<Real> &c, 905 const Vector<Real> &x, 906 Real &tol) { 907 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>( 908 dynamic_cast<const Vector<Real>&>(x)); 909 value(c,*(xs.get_1()),*(xs.get_2()),tol); 910 } 911 912 applyJacobian(Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)913 virtual void applyJacobian(Vector<Real> &jv, 914 const Vector<Real> &v, 915 const Vector<Real> &x, 916 Real &tol) { 917 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>( 918 dynamic_cast<const Vector<Real>&>(x)); 919 const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>( 920 dynamic_cast<const Vector<Real>&>(v)); 921 applyJacobian_1(jv,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol); 922 Ptr<Vector<Real>> jv2 = jv.clone(); 923 applyJacobian_2(*jv2,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol); 924 jv.plus(*jv2); 925 } 926 927 928 using Constraint<Real>::applyAdjointJacobian; applyAdjointJacobian(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)929 virtual void applyAdjointJacobian(Vector<Real> &ajv, 930 const Vector<Real> &v, 931 const Vector<Real> &x, 932 Real &tol) { 933 Vector_SimOpt<Real> &ajvs = dynamic_cast<Vector_SimOpt<Real>&>( 934 dynamic_cast<Vector<Real>&>(ajv)); 935 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>( 936 dynamic_cast<const Vector<Real>&>(x)); 937 Ptr<Vector<Real>> ajv1 = (ajvs.get_1())->clone(); 938 applyAdjointJacobian_1(*ajv1,v,*(xs.get_1()),*(xs.get_2()),tol); 939 ajvs.set_1(*ajv1); 940 Ptr<Vector<Real>> ajv2 = (ajvs.get_2())->clone(); 941 applyAdjointJacobian_2(*ajv2,v,*(xs.get_1()),*(xs.get_2()),tol); 942 ajvs.set_2(*ajv2); 943 } 944 945 applyAdjointHessian(Vector<Real> & ahwv,const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & x,Real & tol)946 virtual void applyAdjointHessian(Vector<Real> &ahwv, 947 const Vector<Real> &w, 948 const Vector<Real> &v, 949 const Vector<Real> &x, 950 Real &tol) { 951 Vector_SimOpt<Real> &ahwvs = dynamic_cast<Vector_SimOpt<Real>&>( 952 dynamic_cast<Vector<Real>&>(ahwv)); 953 const Vector_SimOpt<Real> &xs = dynamic_cast<const Vector_SimOpt<Real>&>( 954 dynamic_cast<const Vector<Real>&>(x)); 955 const Vector_SimOpt<Real> &vs = dynamic_cast<const Vector_SimOpt<Real>&>( 956 dynamic_cast<const Vector<Real>&>(v)); 957 // Block-row 1 958 Ptr<Vector<Real>> C11 = (ahwvs.get_1())->clone(); 959 Ptr<Vector<Real>> C21 = (ahwvs.get_1())->clone(); 960 applyAdjointHessian_11(*C11,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol); 961 applyAdjointHessian_21(*C21,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol); 962 C11->plus(*C21); 963 ahwvs.set_1(*C11); 964 // Block-row 2 965 Ptr<Vector<Real>> C12 = (ahwvs.get_2())->clone(); 966 Ptr<Vector<Real>> C22 = (ahwvs.get_2())->clone(); 967 applyAdjointHessian_12(*C12,w,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol); 968 applyAdjointHessian_22(*C22,w,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol); 969 C22->plus(*C12); 970 ahwvs.set_2(*C22); 971 } 972 973 974 checkSolve(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & c,const bool printToStream=true,std::ostream & outStream=std::cout)975 virtual Real checkSolve(const Vector<Real> &u, 976 const Vector<Real> &z, 977 const Vector<Real> &c, 978 const bool printToStream = true, 979 std::ostream & outStream = std::cout) { 980 // Solve constraint for u. 981 Real tol = ROL_EPSILON<Real>(); 982 Ptr<Vector<Real>> r = c.clone(); 983 Ptr<Vector<Real>> s = u.clone(); 984 solve(*r,*s,z,tol); 985 // Evaluate constraint residual at (u,z). 986 Ptr<Vector<Real>> cs = c.clone(); 987 update(*s,z); 988 value(*cs,*s,z,tol); 989 // Output norm of residual. 990 Real rnorm = r->norm(); 991 Real cnorm = cs->norm(); 992 if ( printToStream ) { 993 std::stringstream hist; 994 hist << std::scientific << std::setprecision(8); 995 hist << "\nTest SimOpt solve at feasible (u,z):\n"; 996 hist << " Solver Residual = " << rnorm << "\n"; 997 hist << " ||c(u,z)|| = " << cnorm << "\n"; 998 outStream << hist.str(); 999 } 1000 return cnorm; 1001 } 1002 1003 1004 /** \brief Check the consistency of the Jacobian and its adjoint. 1005 This is the primary interface. 1006 1007 @param[out] w is a dual constraint-space vector 1008 @param[in] v is a simulation-space vector 1009 @param[in] u is the constraint argument; a simulation-space vector 1010 @param[in] z is the constraint argument; an optimization-space vector 1011 @param[in] printToStream is is a flag that turns on/off output 1012 @param[in] outStream is the output stream 1013 1014 --- 1015 */ checkAdjointConsistencyJacobian_1(const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const bool printToStream=true,std::ostream & outStream=std::cout)1016 virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w, 1017 const Vector<Real> &v, 1018 const Vector<Real> &u, 1019 const Vector<Real> &z, 1020 const bool printToStream = true, 1021 std::ostream & outStream = std::cout) { 1022 return checkAdjointConsistencyJacobian_1(w, v, u, z, w.dual(), v.dual(), printToStream, outStream); 1023 } 1024 1025 1026 /** \brief Check the consistency of the Jacobian and its adjoint. 1027 This is the secondary interface, for use with dual spaces where 1028 the user does not define the dual() operation. 1029 1030 @param[out] w is a dual constraint-space vector 1031 @param[in] v is a simulation-space vector u_lo->zero(); 1032 u_up->setScalar( height ); 1033 1034 @param[in] u is the constraint argument; a simulation-space vector 1035 @param[in] z is the constraint argument; an optimization-space vector 1036 @param[in] dualw is a constraint-space vector 1037 @param[in] dualv is a dual simulation-space vector 1038 @param[in] printToStream is is a flag that turns on/off output 1039 @param[in] outStream is the output stream 1040 1041 --- 1042 */ checkAdjointConsistencyJacobian_1(const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & dualw,const Vector<Real> & dualv,const bool printToStream=true,std::ostream & outStream=std::cout)1043 virtual Real checkAdjointConsistencyJacobian_1(const Vector<Real> &w, 1044 const Vector<Real> &v, 1045 const Vector<Real> &u, 1046 const Vector<Real> &z, 1047 const Vector<Real> &dualw, 1048 const Vector<Real> &dualv, 1049 const bool printToStream = true, 1050 std::ostream & outStream = std::cout) { 1051 Real tol = ROL_EPSILON<Real>(); 1052 Ptr<Vector<Real>> Jv = dualw.clone(); 1053 update(u,z); 1054 applyJacobian_1(*Jv,v,u,z,tol); 1055 Real wJv = w.dot(Jv->dual()); 1056 Ptr<Vector<Real>> Jw = dualv.clone(); 1057 update(u,z); 1058 applyAdjointJacobian_1(*Jw,w,u,z,tol); 1059 Real vJw = v.dot(Jw->dual()); 1060 Real diff = std::abs(wJv-vJw); 1061 if ( printToStream ) { 1062 std::stringstream hist; 1063 hist << std::scientific << std::setprecision(8); 1064 hist << "\nTest SimOpt consistency of Jacobian_1 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = " 1065 << diff << "\n"; 1066 hist << " |<w,Jv>| = " << std::abs(wJv) << "\n"; 1067 hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n"; 1068 outStream << hist.str(); 1069 } 1070 return diff; 1071 } 1072 1073 1074 /** \brief Check the consistency of the Jacobian and its adjoint. 1075 This is the primary interface. 1076 1077 @param[out] w is a dual constraint-space vector 1078 @param[in] v is an optimization-space vector 1079 @param[in] u is the constraint argument; a simulation-space vector 1080 @param[in] z is the constraint argument; an optimization-space vector 1081 @param[in] printToStream is is a flag that turns on/off output 1082 @param[in] outStream is the output stream 1083 1084 --- 1085 */ checkAdjointConsistencyJacobian_2(const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const bool printToStream=true,std::ostream & outStream=std::cout)1086 virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w, 1087 const Vector<Real> &v, 1088 const Vector<Real> &u, 1089 const Vector<Real> &z, 1090 const bool printToStream = true, 1091 std::ostream & outStream = std::cout) { 1092 return checkAdjointConsistencyJacobian_2(w, v, u, z, w.dual(), v.dual(), printToStream, outStream); 1093 } 1094 1095 /** \brief Check the consistency of the Jacobian and its adjoint. 1096 This is the secondary interface, for use with dual spaces where 1097 the user does not define the dual() operation. 1098 1099 @param[out] w is a dual constraint-space vector 1100 @param[in] v is an optimization-space vector 1101 @param[in] u is the constraint argument; a simulation-space vector 1102 @param[in] z is the constraint argument; an optimization-space vector 1103 @param[in] dualw is a constraint-space vector 1104 @param[in] dualv is a dual optimization-space vector 1105 @param[in] printToStream is is a flag that turns on/off output 1106 @param[in] outStream is the output stream 1107 1108 --- 1109 */ checkAdjointConsistencyJacobian_2(const Vector<Real> & w,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & dualw,const Vector<Real> & dualv,const bool printToStream=true,std::ostream & outStream=std::cout)1110 virtual Real checkAdjointConsistencyJacobian_2(const Vector<Real> &w, 1111 const Vector<Real> &v, 1112 const Vector<Real> &u, 1113 const Vector<Real> &z, 1114 const Vector<Real> &dualw, 1115 const Vector<Real> &dualv, 1116 const bool printToStream = true, 1117 std::ostream & outStream = std::cout) { 1118 Real tol = ROL_EPSILON<Real>(); 1119 Ptr<Vector<Real>> Jv = dualw.clone(); 1120 update(u,z); 1121 applyJacobian_2(*Jv,v,u,z,tol); 1122 Real wJv = w.dot(Jv->dual()); 1123 Ptr<Vector<Real>> Jw = dualv.clone(); 1124 update(u,z); 1125 applyAdjointJacobian_2(*Jw,w,u,z,tol); 1126 Real vJw = v.dot(Jw->dual()); 1127 Real diff = std::abs(wJv-vJw); 1128 if ( printToStream ) { 1129 std::stringstream hist; 1130 hist << std::scientific << std::setprecision(8); 1131 hist << "\nTest SimOpt consistency of Jacobian_2 and its adjoint: \n |<w,Jv> - <adj(J)w,v>| = " 1132 << diff << "\n"; 1133 hist << " |<w,Jv>| = " << std::abs(wJv) << "\n"; 1134 hist << " Relative Error = " << diff / (std::abs(wJv)+ROL_UNDERFLOW<Real>()) << "\n"; 1135 outStream << hist.str(); 1136 } 1137 return diff; 1138 } 1139 checkInverseJacobian_1(const Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const bool printToStream=true,std::ostream & outStream=std::cout)1140 virtual Real checkInverseJacobian_1(const Vector<Real> &jv, 1141 const Vector<Real> &v, 1142 const Vector<Real> &u, 1143 const Vector<Real> &z, 1144 const bool printToStream = true, 1145 std::ostream & outStream = std::cout) { 1146 Real tol = ROL_EPSILON<Real>(); 1147 Ptr<Vector<Real>> Jv = jv.clone(); 1148 update(u,z); 1149 applyJacobian_1(*Jv,v,u,z,tol); 1150 Ptr<Vector<Real>> iJJv = u.clone(); 1151 update(u,z); // Does this update do anything? 1152 applyInverseJacobian_1(*iJJv,*Jv,u,z,tol); 1153 Ptr<Vector<Real>> diff = v.clone(); 1154 diff->set(v); 1155 diff->axpy(-1.0,*iJJv); 1156 Real dnorm = diff->norm(); 1157 Real vnorm = v.norm(); 1158 if ( printToStream ) { 1159 std::stringstream hist; 1160 hist << std::scientific << std::setprecision(8); 1161 hist << "\nTest SimOpt consistency of inverse Jacobian_1: \n ||v-inv(J)Jv|| = " 1162 << dnorm << "\n"; 1163 hist << " ||v|| = " << vnorm << "\n"; 1164 hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n"; 1165 outStream << hist.str(); 1166 } 1167 return dnorm; 1168 } 1169 checkInverseAdjointJacobian_1(const Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & u,const Vector<Real> & z,const bool printToStream=true,std::ostream & outStream=std::cout)1170 virtual Real checkInverseAdjointJacobian_1(const Vector<Real> &jv, 1171 const Vector<Real> &v, 1172 const Vector<Real> &u, 1173 const Vector<Real> &z, 1174 const bool printToStream = true, 1175 std::ostream & outStream = std::cout) { 1176 Real tol = ROL_EPSILON<Real>(); 1177 Ptr<Vector<Real>> Jv = jv.clone(); 1178 update(u,z); 1179 applyAdjointJacobian_1(*Jv,v,u,z,tol); 1180 Ptr<Vector<Real>> iJJv = v.clone(); 1181 update(u,z); 1182 applyInverseAdjointJacobian_1(*iJJv,*Jv,u,z,tol); 1183 Ptr<Vector<Real>> diff = v.clone(); 1184 diff->set(v); 1185 diff->axpy(-1.0,*iJJv); 1186 Real dnorm = diff->norm(); 1187 Real vnorm = v.norm(); 1188 if ( printToStream ) { 1189 std::stringstream hist; 1190 hist << std::scientific << std::setprecision(8); 1191 hist << "\nTest SimOpt consistency of inverse adjoint Jacobian_1: \n ||v-inv(adj(J))adj(J)v|| = " 1192 << dnorm << "\n"; 1193 hist << " ||v|| = " << vnorm << "\n"; 1194 hist << " Relative Error = " << dnorm / (vnorm+ROL_UNDERFLOW<Real>()) << "\n"; 1195 outStream << hist.str(); 1196 } 1197 return dnorm; 1198 } 1199 1200 1201 checkApplyJacobian_1(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & v,const Vector<Real> & jv,const bool printToStream=true,std::ostream & outStream=std::cout,const int numSteps=ROL_NUM_CHECKDERIV_STEPS,const int order=1)1202 std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u, 1203 const Vector<Real> &z, 1204 const Vector<Real> &v, 1205 const Vector<Real> &jv, 1206 const bool printToStream = true, 1207 std::ostream & outStream = std::cout, 1208 const int numSteps = ROL_NUM_CHECKDERIV_STEPS, 1209 const int order = 1) { 1210 std::vector<Real> steps(numSteps); 1211 for(int i=0;i<numSteps;++i) { 1212 steps[i] = pow(10,-i); 1213 } 1214 1215 return checkApplyJacobian_1(u,z,v,jv,steps,printToStream,outStream,order); 1216 } 1217 1218 1219 1220 checkApplyJacobian_1(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & v,const Vector<Real> & jv,const std::vector<Real> & steps,const bool printToStream=true,std::ostream & outStream=std::cout,const int order=1)1221 std::vector<std::vector<Real>> checkApplyJacobian_1(const Vector<Real> &u, 1222 const Vector<Real> &z, 1223 const Vector<Real> &v, 1224 const Vector<Real> &jv, 1225 const std::vector<Real> &steps, 1226 const bool printToStream = true, 1227 std::ostream & outStream = std::cout, 1228 const int order = 1) { 1229 1230 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument, 1231 "Error: finite difference order must be 1,2,3, or 4" ); 1232 1233 Real one(1.0); 1234 1235 using Finite_Difference_Arrays::shifts; 1236 using Finite_Difference_Arrays::weights; 1237 1238 Real tol = std::sqrt(ROL_EPSILON<Real>()); 1239 1240 int numSteps = steps.size(); 1241 int numVals = 4; 1242 std::vector<Real> tmp(numVals); 1243 std::vector<std::vector<Real>> jvCheck(numSteps, tmp); 1244 1245 // Save the format state of the original outStream. 1246 nullstream oldFormatState; 1247 oldFormatState.copyfmt(outStream); 1248 1249 // Compute constraint value at x. 1250 Ptr<Vector<Real>> c = jv.clone(); 1251 this->update(u,z); 1252 this->value(*c, u, z, tol); 1253 1254 // Compute (Jacobian at x) times (vector v). 1255 Ptr<Vector<Real>> Jv = jv.clone(); 1256 this->applyJacobian_1(*Jv, v, u, z, tol); 1257 Real normJv = Jv->norm(); 1258 1259 // Temporary vectors. 1260 Ptr<Vector<Real>> cdif = jv.clone(); 1261 Ptr<Vector<Real>> cnew = jv.clone(); 1262 Ptr<Vector<Real>> unew = u.clone(); 1263 1264 for (int i=0; i<numSteps; i++) { 1265 1266 Real eta = steps[i]; 1267 1268 unew->set(u); 1269 1270 cdif->set(*c); 1271 cdif->scale(weights[order-1][0]); 1272 1273 for(int j=0; j<order; ++j) { 1274 1275 unew->axpy(eta*shifts[order-1][j], v); 1276 1277 if( weights[order-1][j+1] != 0 ) { 1278 this->update(*unew,z); 1279 this->value(*cnew,*unew,z,tol); 1280 cdif->axpy(weights[order-1][j+1],*cnew); 1281 } 1282 1283 } 1284 1285 cdif->scale(one/eta); 1286 1287 // Compute norms of Jacobian-vector products, finite-difference approximations, and error. 1288 jvCheck[i][0] = eta; 1289 jvCheck[i][1] = normJv; 1290 jvCheck[i][2] = cdif->norm(); 1291 cdif->axpy(-one, *Jv); 1292 jvCheck[i][3] = cdif->norm(); 1293 1294 if (printToStream) { 1295 std::stringstream hist; 1296 if (i==0) { 1297 hist << std::right 1298 << std::setw(20) << "Step size" 1299 << std::setw(20) << "norm(Jac*vec)" 1300 << std::setw(20) << "norm(FD approx)" 1301 << std::setw(20) << "norm(abs error)" 1302 << "\n" 1303 << std::setw(20) << "---------" 1304 << std::setw(20) << "-------------" 1305 << std::setw(20) << "---------------" 1306 << std::setw(20) << "---------------" 1307 << "\n"; 1308 } 1309 hist << std::scientific << std::setprecision(11) << std::right 1310 << std::setw(20) << jvCheck[i][0] 1311 << std::setw(20) << jvCheck[i][1] 1312 << std::setw(20) << jvCheck[i][2] 1313 << std::setw(20) << jvCheck[i][3] 1314 << "\n"; 1315 outStream << hist.str(); 1316 } 1317 1318 } 1319 1320 // Reset format state of outStream. 1321 outStream.copyfmt(oldFormatState); 1322 1323 return jvCheck; 1324 } // checkApplyJacobian 1325 1326 checkApplyJacobian_2(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & v,const Vector<Real> & jv,const bool printToStream=true,std::ostream & outStream=std::cout,const int numSteps=ROL_NUM_CHECKDERIV_STEPS,const int order=1)1327 std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u, 1328 const Vector<Real> &z, 1329 const Vector<Real> &v, 1330 const Vector<Real> &jv, 1331 const bool printToStream = true, 1332 std::ostream & outStream = std::cout, 1333 const int numSteps = ROL_NUM_CHECKDERIV_STEPS, 1334 const int order = 1) { 1335 std::vector<Real> steps(numSteps); 1336 for(int i=0;i<numSteps;++i) { 1337 steps[i] = pow(10,-i); 1338 } 1339 1340 return checkApplyJacobian_2(u,z,v,jv,steps,printToStream,outStream,order); 1341 } 1342 1343 1344 1345 checkApplyJacobian_2(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & v,const Vector<Real> & jv,const std::vector<Real> & steps,const bool printToStream=true,std::ostream & outStream=std::cout,const int order=1)1346 std::vector<std::vector<Real>> checkApplyJacobian_2(const Vector<Real> &u, 1347 const Vector<Real> &z, 1348 const Vector<Real> &v, 1349 const Vector<Real> &jv, 1350 const std::vector<Real> &steps, 1351 const bool printToStream = true, 1352 std::ostream & outStream = std::cout, 1353 const int order = 1) { 1354 1355 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument, 1356 "Error: finite difference order must be 1,2,3, or 4" ); 1357 1358 Real one(1.0); 1359 1360 using Finite_Difference_Arrays::shifts; 1361 using Finite_Difference_Arrays::weights; 1362 1363 Real tol = std::sqrt(ROL_EPSILON<Real>()); 1364 1365 int numSteps = steps.size(); 1366 int numVals = 4; 1367 std::vector<Real> tmp(numVals); 1368 std::vector<std::vector<Real>> jvCheck(numSteps, tmp); 1369 1370 // Save the format state of the original outStream. 1371 nullstream oldFormatState; 1372 oldFormatState.copyfmt(outStream); 1373 1374 // Compute constraint value at x. 1375 Ptr<Vector<Real>> c = jv.clone(); 1376 this->update(u,z); 1377 this->value(*c, u, z, tol); 1378 1379 // Compute (Jacobian at x) times (vector v). 1380 Ptr<Vector<Real>> Jv = jv.clone(); 1381 this->applyJacobian_2(*Jv, v, u, z, tol); 1382 Real normJv = Jv->norm(); 1383 1384 // Temporary vectors. 1385 Ptr<Vector<Real>> cdif = jv.clone(); 1386 Ptr<Vector<Real>> cnew = jv.clone(); 1387 Ptr<Vector<Real>> znew = z.clone(); 1388 1389 for (int i=0; i<numSteps; i++) { 1390 1391 Real eta = steps[i]; 1392 1393 znew->set(z); 1394 1395 cdif->set(*c); 1396 cdif->scale(weights[order-1][0]); 1397 1398 for(int j=0; j<order; ++j) { 1399 1400 znew->axpy(eta*shifts[order-1][j], v); 1401 1402 if( weights[order-1][j+1] != 0 ) { 1403 this->update(u,*znew); 1404 this->value(*cnew,u,*znew,tol); 1405 cdif->axpy(weights[order-1][j+1],*cnew); 1406 } 1407 1408 } 1409 1410 cdif->scale(one/eta); 1411 1412 // Compute norms of Jacobian-vector products, finite-difference approximations, and error. 1413 jvCheck[i][0] = eta; 1414 jvCheck[i][1] = normJv; 1415 jvCheck[i][2] = cdif->norm(); 1416 cdif->axpy(-one, *Jv); 1417 jvCheck[i][3] = cdif->norm(); 1418 1419 if (printToStream) { 1420 std::stringstream hist; 1421 if (i==0) { 1422 hist << std::right 1423 << std::setw(20) << "Step size" 1424 << std::setw(20) << "norm(Jac*vec)" 1425 << std::setw(20) << "norm(FD approx)" 1426 << std::setw(20) << "norm(abs error)" 1427 << "\n" 1428 << std::setw(20) << "---------" 1429 << std::setw(20) << "-------------" 1430 << std::setw(20) << "---------------" 1431 << std::setw(20) << "---------------" 1432 << "\n"; 1433 } 1434 hist << std::scientific << std::setprecision(11) << std::right 1435 << std::setw(20) << jvCheck[i][0] 1436 << std::setw(20) << jvCheck[i][1] 1437 << std::setw(20) << jvCheck[i][2] 1438 << std::setw(20) << jvCheck[i][3] 1439 << "\n"; 1440 outStream << hist.str(); 1441 } 1442 1443 } 1444 1445 // Reset format state of outStream. 1446 outStream.copyfmt(oldFormatState); 1447 1448 return jvCheck; 1449 } // checkApplyJacobian 1450 1451 1452 checkApplyAdjointHessian_11(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const bool printToStream=true,std::ostream & outStream=std::cout,const int numSteps=ROL_NUM_CHECKDERIV_STEPS,const int order=1)1453 std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u, 1454 const Vector<Real> &z, 1455 const Vector<Real> &p, 1456 const Vector<Real> &v, 1457 const Vector<Real> &hv, 1458 const bool printToStream = true, 1459 std::ostream & outStream = std::cout, 1460 const int numSteps = ROL_NUM_CHECKDERIV_STEPS, 1461 const int order = 1 ) { 1462 std::vector<Real> steps(numSteps); 1463 for(int i=0;i<numSteps;++i) { 1464 steps[i] = pow(10,-i); 1465 } 1466 1467 return checkApplyAdjointHessian_11(u,z,p,v,hv,steps,printToStream,outStream,order); 1468 1469 } 1470 checkApplyAdjointHessian_11(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const std::vector<Real> & steps,const bool printToStream=true,std::ostream & outStream=std::cout,const int order=1)1471 std::vector<std::vector<Real>> checkApplyAdjointHessian_11(const Vector<Real> &u, 1472 const Vector<Real> &z, 1473 const Vector<Real> &p, 1474 const Vector<Real> &v, 1475 const Vector<Real> &hv, 1476 const std::vector<Real> &steps, 1477 const bool printToStream = true, 1478 std::ostream & outStream = std::cout, 1479 const int order = 1 ) { 1480 using Finite_Difference_Arrays::shifts; 1481 using Finite_Difference_Arrays::weights; 1482 1483 Real one(1.0); 1484 1485 Real tol = std::sqrt(ROL_EPSILON<Real>()); 1486 1487 int numSteps = steps.size(); 1488 int numVals = 4; 1489 std::vector<Real> tmp(numVals); 1490 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp); 1491 1492 // Temporary vectors. 1493 Ptr<Vector<Real>> AJdif = hv.clone(); 1494 Ptr<Vector<Real>> AJp = hv.clone(); 1495 Ptr<Vector<Real>> AHpv = hv.clone(); 1496 Ptr<Vector<Real>> AJnew = hv.clone(); 1497 Ptr<Vector<Real>> unew = u.clone(); 1498 1499 // Save the format state of the original outStream. 1500 nullstream oldFormatState; 1501 oldFormatState.copyfmt(outStream); 1502 1503 // Apply adjoint Jacobian to p. 1504 this->update(u,z); 1505 this->applyAdjointJacobian_1(*AJp, p, u, z, tol); 1506 1507 // Apply adjoint Hessian at (u,z), in direction v, to p. 1508 this->applyAdjointHessian_11(*AHpv, p, v, u, z, tol); 1509 Real normAHpv = AHpv->norm(); 1510 1511 for (int i=0; i<numSteps; i++) { 1512 1513 Real eta = steps[i]; 1514 1515 // Apply adjoint Jacobian to p at (u+eta*v,z). 1516 unew->set(u); 1517 1518 AJdif->set(*AJp); 1519 AJdif->scale(weights[order-1][0]); 1520 1521 for(int j=0; j<order; ++j) { 1522 1523 unew->axpy(eta*shifts[order-1][j],v); 1524 1525 if( weights[order-1][j+1] != 0 ) { 1526 this->update(*unew,z); 1527 this->applyAdjointJacobian_1(*AJnew, p, *unew, z, tol); 1528 AJdif->axpy(weights[order-1][j+1],*AJnew); 1529 } 1530 } 1531 1532 AJdif->scale(one/eta); 1533 1534 // Compute norms of Jacobian-vector products, finite-difference approximations, and error. 1535 ahpvCheck[i][0] = eta; 1536 ahpvCheck[i][1] = normAHpv; 1537 ahpvCheck[i][2] = AJdif->norm(); 1538 AJdif->axpy(-one, *AHpv); 1539 ahpvCheck[i][3] = AJdif->norm(); 1540 1541 if (printToStream) { 1542 std::stringstream hist; 1543 if (i==0) { 1544 hist << std::right 1545 << std::setw(20) << "Step size" 1546 << std::setw(20) << "norm(adj(H)(u,v))" 1547 << std::setw(20) << "norm(FD approx)" 1548 << std::setw(20) << "norm(abs error)" 1549 << "\n" 1550 << std::setw(20) << "---------" 1551 << std::setw(20) << "-----------------" 1552 << std::setw(20) << "---------------" 1553 << std::setw(20) << "---------------" 1554 << "\n"; 1555 } 1556 hist << std::scientific << std::setprecision(11) << std::right 1557 << std::setw(20) << ahpvCheck[i][0] 1558 << std::setw(20) << ahpvCheck[i][1] 1559 << std::setw(20) << ahpvCheck[i][2] 1560 << std::setw(20) << ahpvCheck[i][3] 1561 << "\n"; 1562 outStream << hist.str(); 1563 } 1564 1565 } 1566 1567 // Reset format state of outStream. 1568 outStream.copyfmt(oldFormatState); 1569 1570 return ahpvCheck; 1571 } // checkApplyAdjointHessian_11 1572 1573 /** 1574 \brief \f$ u\in U \f$, \f$ z\in Z \f$, \f$ p\in C^\ast \f$, \f$ v \in U \f$, \f$ hv \in U^\ast \f$ 1575 */ checkApplyAdjointHessian_21(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const bool printToStream=true,std::ostream & outStream=std::cout,const int numSteps=ROL_NUM_CHECKDERIV_STEPS,const int order=1)1576 std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u, 1577 const Vector<Real> &z, 1578 const Vector<Real> &p, 1579 const Vector<Real> &v, 1580 const Vector<Real> &hv, 1581 const bool printToStream = true, 1582 std::ostream & outStream = std::cout, 1583 const int numSteps = ROL_NUM_CHECKDERIV_STEPS, 1584 const int order = 1 ) { 1585 std::vector<Real> steps(numSteps); 1586 for(int i=0;i<numSteps;++i) { 1587 steps[i] = pow(10,-i); 1588 } 1589 1590 return checkApplyAdjointHessian_21(u,z,p,v,hv,steps,printToStream,outStream,order); 1591 1592 } 1593 1594 /** 1595 \brief \f$ u\in U \f$, \f$ z\in Z \f$, \f$ p\in C^\ast \f$, \f$ v \in U \f$, \f$ hv \in U^\ast \f$ 1596 */ checkApplyAdjointHessian_21(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const std::vector<Real> & steps,const bool printToStream=true,std::ostream & outStream=std::cout,const int order=1)1597 std::vector<std::vector<Real>> checkApplyAdjointHessian_21(const Vector<Real> &u, 1598 const Vector<Real> &z, 1599 const Vector<Real> &p, 1600 const Vector<Real> &v, 1601 const Vector<Real> &hv, 1602 const std::vector<Real> &steps, 1603 const bool printToStream = true, 1604 std::ostream & outStream = std::cout, 1605 const int order = 1 ) { 1606 using Finite_Difference_Arrays::shifts; 1607 using Finite_Difference_Arrays::weights; 1608 1609 Real one(1.0); 1610 1611 Real tol = std::sqrt(ROL_EPSILON<Real>()); 1612 1613 int numSteps = steps.size(); 1614 int numVals = 4; 1615 std::vector<Real> tmp(numVals); 1616 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp); 1617 1618 // Temporary vectors. 1619 Ptr<Vector<Real>> AJdif = hv.clone(); 1620 Ptr<Vector<Real>> AJp = hv.clone(); 1621 Ptr<Vector<Real>> AHpv = hv.clone(); 1622 Ptr<Vector<Real>> AJnew = hv.clone(); 1623 Ptr<Vector<Real>> znew = z.clone(); 1624 1625 // Save the format state of the original outStream. 1626 nullstream oldFormatState; 1627 oldFormatState.copyfmt(outStream); 1628 1629 // Apply adjoint Jacobian to p. 1630 this->update(u,z); 1631 this->applyAdjointJacobian_1(*AJp, p, u, z, tol); 1632 1633 // Apply adjoint Hessian at (u,z), in direction v, to p. 1634 this->applyAdjointHessian_21(*AHpv, p, v, u, z, tol); 1635 Real normAHpv = AHpv->norm(); 1636 1637 for (int i=0; i<numSteps; i++) { 1638 1639 Real eta = steps[i]; 1640 1641 // Apply adjoint Jacobian to p at (u,z+eta*v). 1642 znew->set(z); 1643 1644 AJdif->set(*AJp); 1645 AJdif->scale(weights[order-1][0]); 1646 1647 for(int j=0; j<order; ++j) { 1648 1649 znew->axpy(eta*shifts[order-1][j],v); 1650 1651 if( weights[order-1][j+1] != 0 ) { 1652 this->update(u,*znew); 1653 this->applyAdjointJacobian_1(*AJnew, p, u, *znew, tol); 1654 AJdif->axpy(weights[order-1][j+1],*AJnew); 1655 } 1656 } 1657 1658 AJdif->scale(one/eta); 1659 1660 // Compute norms of Jacobian-vector products, finite-difference approximations, and error. 1661 ahpvCheck[i][0] = eta; 1662 ahpvCheck[i][1] = normAHpv; 1663 ahpvCheck[i][2] = AJdif->norm(); 1664 AJdif->axpy(-one, *AHpv); 1665 ahpvCheck[i][3] = AJdif->norm(); 1666 1667 if (printToStream) { 1668 std::stringstream hist; 1669 if (i==0) { 1670 hist << std::right 1671 << std::setw(20) << "Step size" 1672 << std::setw(20) << "norm(adj(H)(u,v))" 1673 << std::setw(20) << "norm(FD approx)" 1674 << std::setw(20) << "norm(abs error)" 1675 << "\n" 1676 << std::setw(20) << "---------" 1677 << std::setw(20) << "-----------------" 1678 << std::setw(20) << "---------------" 1679 << std::setw(20) << "---------------" 1680 << "\n"; 1681 } 1682 hist << std::scientific << std::setprecision(11) << std::right 1683 << std::setw(20) << ahpvCheck[i][0] 1684 << std::setw(20) << ahpvCheck[i][1] 1685 << std::setw(20) << ahpvCheck[i][2] 1686 << std::setw(20) << ahpvCheck[i][3] 1687 << "\n"; 1688 outStream << hist.str(); 1689 } 1690 1691 } 1692 1693 // Reset format state of outStream. 1694 outStream.copyfmt(oldFormatState); 1695 1696 return ahpvCheck; 1697 } // checkApplyAdjointHessian_21 1698 1699 /** 1700 \brief \f$ u\in U \f$, \f$ z\in Z \f$, \f$ p\in C^\ast \f$, \f$ v \in U \f$, \f$ hv \in U^\ast \f$ 1701 */ checkApplyAdjointHessian_12(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const bool printToStream=true,std::ostream & outStream=std::cout,const int numSteps=ROL_NUM_CHECKDERIV_STEPS,const int order=1)1702 std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u, 1703 const Vector<Real> &z, 1704 const Vector<Real> &p, 1705 const Vector<Real> &v, 1706 const Vector<Real> &hv, 1707 const bool printToStream = true, 1708 std::ostream & outStream = std::cout, 1709 const int numSteps = ROL_NUM_CHECKDERIV_STEPS, 1710 const int order = 1 ) { 1711 std::vector<Real> steps(numSteps); 1712 for(int i=0;i<numSteps;++i) { 1713 steps[i] = pow(10,-i); 1714 } 1715 1716 return checkApplyAdjointHessian_12(u,z,p,v,hv,steps,printToStream,outStream,order); 1717 1718 } 1719 1720 checkApplyAdjointHessian_12(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const std::vector<Real> & steps,const bool printToStream=true,std::ostream & outStream=std::cout,const int order=1)1721 std::vector<std::vector<Real>> checkApplyAdjointHessian_12(const Vector<Real> &u, 1722 const Vector<Real> &z, 1723 const Vector<Real> &p, 1724 const Vector<Real> &v, 1725 const Vector<Real> &hv, 1726 const std::vector<Real> &steps, 1727 const bool printToStream = true, 1728 std::ostream & outStream = std::cout, 1729 const int order = 1 ) { 1730 using Finite_Difference_Arrays::shifts; 1731 using Finite_Difference_Arrays::weights; 1732 1733 Real one(1.0); 1734 1735 Real tol = std::sqrt(ROL_EPSILON<Real>()); 1736 1737 int numSteps = steps.size(); 1738 int numVals = 4; 1739 std::vector<Real> tmp(numVals); 1740 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp); 1741 1742 // Temporary vectors. 1743 Ptr<Vector<Real>> AJdif = hv.clone(); 1744 Ptr<Vector<Real>> AJp = hv.clone(); 1745 Ptr<Vector<Real>> AHpv = hv.clone(); 1746 Ptr<Vector<Real>> AJnew = hv.clone(); 1747 Ptr<Vector<Real>> unew = u.clone(); 1748 1749 // Save the format state of the original outStream. 1750 nullstream oldFormatState; 1751 oldFormatState.copyfmt(outStream); 1752 1753 // Apply adjoint Jacobian to p. 1754 this->update(u,z); 1755 this->applyAdjointJacobian_2(*AJp, p, u, z, tol); 1756 1757 // Apply adjoint Hessian at (u,z), in direction v, to p. 1758 this->applyAdjointHessian_12(*AHpv, p, v, u, z, tol); 1759 Real normAHpv = AHpv->norm(); 1760 1761 for (int i=0; i<numSteps; i++) { 1762 1763 Real eta = steps[i]; 1764 1765 // Apply adjoint Jacobian to p at (u+eta*v,z). 1766 unew->set(u); 1767 1768 AJdif->set(*AJp); 1769 AJdif->scale(weights[order-1][0]); 1770 1771 for(int j=0; j<order; ++j) { 1772 1773 unew->axpy(eta*shifts[order-1][j],v); 1774 1775 if( weights[order-1][j+1] != 0 ) { 1776 this->update(*unew,z); 1777 this->applyAdjointJacobian_2(*AJnew, p, *unew, z, tol); 1778 AJdif->axpy(weights[order-1][j+1],*AJnew); 1779 } 1780 } 1781 1782 AJdif->scale(one/eta); 1783 1784 // Compute norms of Jacobian-vector products, finite-difference approximations, and error. 1785 ahpvCheck[i][0] = eta; 1786 ahpvCheck[i][1] = normAHpv; 1787 ahpvCheck[i][2] = AJdif->norm(); 1788 AJdif->axpy(-one, *AHpv); 1789 ahpvCheck[i][3] = AJdif->norm(); 1790 1791 if (printToStream) { 1792 std::stringstream hist; 1793 if (i==0) { 1794 hist << std::right 1795 << std::setw(20) << "Step size" 1796 << std::setw(20) << "norm(adj(H)(u,v))" 1797 << std::setw(20) << "norm(FD approx)" 1798 << std::setw(20) << "norm(abs error)" 1799 << "\n" 1800 << std::setw(20) << "---------" 1801 << std::setw(20) << "-----------------" 1802 << std::setw(20) << "---------------" 1803 << std::setw(20) << "---------------" 1804 << "\n"; 1805 } 1806 hist << std::scientific << std::setprecision(11) << std::right 1807 << std::setw(20) << ahpvCheck[i][0] 1808 << std::setw(20) << ahpvCheck[i][1] 1809 << std::setw(20) << ahpvCheck[i][2] 1810 << std::setw(20) << ahpvCheck[i][3] 1811 << "\n"; 1812 outStream << hist.str(); 1813 } 1814 1815 } 1816 1817 // Reset format state of outStream. 1818 outStream.copyfmt(oldFormatState); 1819 1820 return ahpvCheck; 1821 } // checkApplyAdjointHessian_12 1822 checkApplyAdjointHessian_22(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const bool printToStream=true,std::ostream & outStream=std::cout,const int numSteps=ROL_NUM_CHECKDERIV_STEPS,const int order=1)1823 std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u, 1824 const Vector<Real> &z, 1825 const Vector<Real> &p, 1826 const Vector<Real> &v, 1827 const Vector<Real> &hv, 1828 const bool printToStream = true, 1829 std::ostream & outStream = std::cout, 1830 const int numSteps = ROL_NUM_CHECKDERIV_STEPS, 1831 const int order = 1 ) { 1832 std::vector<Real> steps(numSteps); 1833 for(int i=0;i<numSteps;++i) { 1834 steps[i] = pow(10,-i); 1835 } 1836 1837 return checkApplyAdjointHessian_22(u,z,p,v,hv,steps,printToStream,outStream,order); 1838 1839 } 1840 checkApplyAdjointHessian_22(const Vector<Real> & u,const Vector<Real> & z,const Vector<Real> & p,const Vector<Real> & v,const Vector<Real> & hv,const std::vector<Real> & steps,const bool printToStream=true,std::ostream & outStream=std::cout,const int order=1)1841 std::vector<std::vector<Real>> checkApplyAdjointHessian_22(const Vector<Real> &u, 1842 const Vector<Real> &z, 1843 const Vector<Real> &p, 1844 const Vector<Real> &v, 1845 const Vector<Real> &hv, 1846 const std::vector<Real> &steps, 1847 const bool printToStream = true, 1848 std::ostream & outStream = std::cout, 1849 const int order = 1 ) { 1850 using Finite_Difference_Arrays::shifts; 1851 using Finite_Difference_Arrays::weights; 1852 1853 Real one(1.0); 1854 1855 Real tol = std::sqrt(ROL_EPSILON<Real>()); 1856 1857 int numSteps = steps.size(); 1858 int numVals = 4; 1859 std::vector<Real> tmp(numVals); 1860 std::vector<std::vector<Real>> ahpvCheck(numSteps, tmp); 1861 1862 // Temporary vectors. 1863 Ptr<Vector<Real>> AJdif = hv.clone(); 1864 Ptr<Vector<Real>> AJp = hv.clone(); 1865 Ptr<Vector<Real>> AHpv = hv.clone(); 1866 Ptr<Vector<Real>> AJnew = hv.clone(); 1867 Ptr<Vector<Real>> znew = z.clone(); 1868 1869 // Save the format state of the original outStream. 1870 nullstream oldFormatState; 1871 oldFormatState.copyfmt(outStream); 1872 1873 // Apply adjoint Jacobian to p. 1874 this->update(u,z); 1875 this->applyAdjointJacobian_2(*AJp, p, u, z, tol); 1876 1877 // Apply adjoint Hessian at (u,z), in direction v, to p. 1878 this->applyAdjointHessian_22(*AHpv, p, v, u, z, tol); 1879 Real normAHpv = AHpv->norm(); 1880 1881 for (int i=0; i<numSteps; i++) { 1882 1883 Real eta = steps[i]; 1884 1885 // Apply adjoint Jacobian to p at (u,z+eta*v). 1886 znew->set(z); 1887 1888 AJdif->set(*AJp); 1889 AJdif->scale(weights[order-1][0]); 1890 1891 for(int j=0; j<order; ++j) { 1892 1893 znew->axpy(eta*shifts[order-1][j],v); 1894 1895 if( weights[order-1][j+1] != 0 ) { 1896 this->update(u,*znew); 1897 this->applyAdjointJacobian_2(*AJnew, p, u, *znew, tol); 1898 AJdif->axpy(weights[order-1][j+1],*AJnew); 1899 } 1900 } 1901 1902 AJdif->scale(one/eta); 1903 1904 // Compute norms of Jacobian-vector products, finite-difference approximations, and error. 1905 ahpvCheck[i][0] = eta; 1906 ahpvCheck[i][1] = normAHpv; 1907 ahpvCheck[i][2] = AJdif->norm(); 1908 AJdif->axpy(-one, *AHpv); 1909 ahpvCheck[i][3] = AJdif->norm(); 1910 1911 if (printToStream) { 1912 std::stringstream hist; 1913 if (i==0) { 1914 hist << std::right 1915 << std::setw(20) << "Step size" 1916 << std::setw(20) << "norm(adj(H)(u,v))" 1917 << std::setw(20) << "norm(FD approx)" 1918 << std::setw(20) << "norm(abs error)" 1919 << "\n" 1920 << std::setw(20) << "---------" 1921 << std::setw(20) << "-----------------" 1922 << std::setw(20) << "---------------" 1923 << std::setw(20) << "---------------" 1924 << "\n"; 1925 } 1926 hist << std::scientific << std::setprecision(11) << std::right 1927 << std::setw(20) << ahpvCheck[i][0] 1928 << std::setw(20) << ahpvCheck[i][1] 1929 << std::setw(20) << ahpvCheck[i][2] 1930 << std::setw(20) << ahpvCheck[i][3] 1931 << "\n"; 1932 outStream << hist.str(); 1933 } 1934 1935 } 1936 1937 // Reset format state of outStream. 1938 outStream.copyfmt(oldFormatState); 1939 1940 return ahpvCheck; 1941 } // checkApplyAdjointHessian_22 1942 1943 }; // class Constraint_SimOpt 1944 1945 } // namespace ROL 1946 1947 #endif 1948