example_01b.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#include "ROL_StdVector.hpp"
#include "Sacado.hpp"
 
using namespace ROL;
 
template<class ScalarT>
class FunctionZakharov {
 
  
 
  public:
 
    ScalarT eval(const std::vector<ScalarT> &x);
 
};
 
template<class ScalarT>
ScalarT FunctionZakharov<ScalarT>::eval(const std::vector<ScalarT> & x) {
 
    typedef typename std::vector<ScalarT>::size_type luint;
 
    ScalarT xdotx = 0;
    ScalarT kdotx = 0;
    ScalarT J = 0;
   
    // Compute dot products 
    for(luint i=0; i<x.size(); ++i) {
        xdotx += pow(x[i],2);       // (k,x)
        kdotx += ScalarT(i+1)*x[i];  // (x,x)
    }
 
    // Sum terms in objective function
    J = xdotx + pow(kdotx,2)/4.0 + pow(kdotx,4)/16.0;
    
    return J;
}
 
 
template<class Real>
class Zakharov_Sacado_Objective : public Objective<Real> {
 
  typedef std::vector<Real>  vector;
  typedef Vector<Real>       V;
  typedef StdVector<Real>    SV;
 
  typedef Sacado::Fad::DFad<Real>          GradType;
  typedef Sacado::Fad::SFad<Real,1>        DirDerivType;
  typedef Sacado::Fad::DFad<DirDerivType>  HessVecType;
 
  typedef typename vector::size_type luint;
 
  // In C++11, we could use template typedefs:
  // template <typename T>       using  GradTypeT = Sacado::Fad::DFad<T>;
  // typedef Sacado::Fad::SFad<Real,1>  DirDerivType;
  // typedef GradTypeT<Real>            GradType;
  // typedef GradTypeT<DirDerivType>    HessVecType;
 
  private:
 
    FunctionZakharov<Real>        zfunc_; 
    FunctionZakharov<GradType>    zfuncGrad_; 
    FunctionZakharov<HessVecType> zfuncHessVec_; 
 
    ROL::Ptr<const vector> getVector( const V& x ) {
      
      return dynamic_cast<const SV&>(x).getVector();
    }
 
    ROL::Ptr<vector> getVector( V& x ) {
      
      return dynamic_cast<SV&>(x).getVector();
    }
 
  public:
 
    Zakharov_Sacado_Objective() {}
 
    /* Evaluate the objective function at x */
    Real value( const Vector<Real> &x, Real &tol ) {
      
      ROL::Ptr<const vector> xp = getVector(x);      
      return zfunc_.eval(*xp);
    }
 
    /* Evaluate the gradient at x */
    void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
      
      ROL::Ptr<const vector> xp = getVector(x);
      ROL::Ptr<vector> gp = getVector(g);
 
      luint n = xp->size();
 
      std::vector<GradType> x_grad(n);
 
      for(luint i=0; i<n; ++i) {
        x_grad[i] = (*xp)[i]; // Set values x(i).
        x_grad[i].diff(i,n);  // Choose canonical directions.
      }
 
      GradType J_grad = zfuncGrad_.eval(x_grad);
 
      for(luint i=0; i<n; ++i) {
        (*gp)[i] = J_grad.dx(i);
      }
 
    }
 
    /* Compute the action of the Hessian evaluated at x on a vector v */
    void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
 
      
      ROL::Ptr<vector> hvp = getVector(hv);
      ROL::Ptr<const vector> vp = getVector(v);
      ROL::Ptr<const vector> xp = getVector(x);
 
      luint n = xp->size();
 
      std::vector<HessVecType>  x_hessvec(n);
 
      for(luint i=0; i<n; ++i) {
        DirDerivType tmp(1,(*xp)[i]);  // Set values x(i).
        tmp.fastAccessDx(0)= (*vp)[i]; // Set direction values v(i).
        x_hessvec[i] = tmp;            // Use tmp to define hessvec-able x.
        x_hessvec[i].diff(i,n);        // Choose directions.
      }
 
      // Compute Hessian-vector product (and other currently irrelevant things).
      HessVecType J_hessvec = zfuncHessVec_.eval(x_hessvec);
 
      for(luint i=0; i<n; ++i) {
        (*hvp)[i]   =  (J_hessvec.dx(i)).fastAccessDx(0);
        // hessvec  =  get gradient      get dir deriv
      }
    }
};