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 /** \file 45 \brief Contains definitions for W. Hock and K. Schittkowski 39th test function. 46 \author Created by G. von Winckel 47 */ 48 49 50 #ifndef ROL_HS39_HPP 51 #define ROL_HS39_HPP 52 53 #include "ROL_StdVector.hpp" 54 #include "ROL_TestProblem.hpp" 55 #include "ROL_Bounds.hpp" 56 #include "ROL_Constraint_Partitioned.hpp" 57 #include "ROL_Types.hpp" 58 59 namespace ROL { 60 namespace ZOO { 61 62 /** \brief W. Hock and K. Schittkowski 39th test function. 63 * 64 * Exact solution x* = (1, 0, 0, 0) 65 * f(x*) = -1 66 */ 67 68 template<class Real> 69 class Objective_HS39 : public Objective<Real> { 70 71 typedef std::vector<Real> vector; 72 73 typedef StdVector<Real> SV; 74 75 76 77 public: 78 value(const Vector<Real> & x,Real & tol)79 Real value( const Vector<Real> &x, Real &tol ) { 80 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 81 return -(*xp)[0]; 82 } 83 gradient(Vector<Real> & g,const Vector<Real> & x,Real & tol)84 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) { 85 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 86 ROL::Ptr<vector> gp = dynamic_cast<SV&>(g).getVector(); 87 88 (*gp)[0] = -1.0; 89 (*gp)[1] = 0.0; 90 (*gp)[2] = 0.0; 91 (*gp)[3] = 0.0; 92 93 } 94 hessVec(Vector<Real> & hv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)95 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) { 96 hv.zero(); 97 } 98 }; 99 100 // First of two equality constraints 101 template<class Real> 102 class Constraint_HS39a : public Constraint<Real> { 103 104 typedef std::vector<Real> vector; 105 typedef StdVector<Real> SV; 106 107 public: Constraint_HS39a(void)108 Constraint_HS39a(void) {} 109 value(Vector<Real> & c,const Vector<Real> & x,Real & tol)110 void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) { 111 112 ROL::Ptr<vector> cp = dynamic_cast<SV&>(c).getVector(); 113 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 114 115 (*cp)[0] = (*xp)[1]-std::pow((*xp)[0],3)-std::pow((*xp)[2],2); 116 } 117 applyJacobian(Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)118 void applyJacobian(Vector<Real> &jv, const Vector<Real> &v, 119 const Vector<Real> &x, Real &tol) { 120 121 ROL::Ptr<vector> jvp = dynamic_cast<SV&>(jv).getVector(); 122 ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector(); 123 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 124 125 (*jvp)[0] = (*vp)[1] - 3.0*(*xp)[0]*(*xp)[0]*(*vp)[0] - 2.0*(*xp)[2]*(*vp)[2]; 126 127 } 128 applyAdjointJacobian(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)129 void applyAdjointJacobian( Vector<Real> &ajv, const Vector<Real> &v, 130 const Vector<Real> &x, Real &tol ) { 131 132 ROL::Ptr<vector> ajvp = dynamic_cast<SV&>(ajv).getVector(); 133 ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector(); 134 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 135 136 (*ajvp)[0] = -3.0*(*xp)[0]*(*xp)[0]*(*vp)[0]; 137 (*ajvp)[1] = (*vp)[0]; 138 (*ajvp)[2] = -2.0*(*xp)[2]*(*vp)[0]; 139 (*ajvp)[3] = 0.0; 140 } 141 applyAdjointHessian(Vector<Real> & ahuv,const Vector<Real> & u,const Vector<Real> & v,const Vector<Real> & x,Real & tol)142 void applyAdjointHessian(Vector<Real> &ahuv, const Vector<Real> &u, 143 const Vector<Real> &v, const Vector<Real> &x, 144 Real &tol) { 145 146 ROL::Ptr<vector> ahuvp = dynamic_cast<SV&>(ahuv).getVector(); 147 ROL::Ptr<const vector> up = dynamic_cast<const SV&>(u).getVector(); 148 ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector(); 149 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 150 151 (*ahuvp)[0] = -6.0*(*up)[0]*(*xp)[0]*(*vp)[0]; 152 (*ahuvp)[1] = 0.0; 153 (*ahuvp)[2] = -2.0*(*up)[0]*(*vp)[2]; 154 (*ahuvp)[3] = 0.0; 155 156 } 157 158 159 }; 160 161 // Second of two equality constraints 162 template<class Real> 163 class Constraint_HS39b : public Constraint<Real> { 164 165 typedef std::vector<Real> vector; 166 typedef StdVector<Real> SV; 167 168 public: Constraint_HS39b(void)169 Constraint_HS39b(void) {} 170 value(Vector<Real> & c,const Vector<Real> & x,Real & tol)171 void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) { 172 ROL::Ptr<vector> cp = dynamic_cast<SV&>(c).getVector(); 173 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 174 175 (*cp)[0] = std::pow((*xp)[0],2)-(*xp)[1]-std::pow((*xp)[3],2); 176 } 177 applyJacobian(Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)178 void applyJacobian(Vector<Real> &jv, const Vector<Real> &v, 179 const Vector<Real> &x, Real &tol) { 180 181 ROL::Ptr<vector> jvp = dynamic_cast<SV&>(jv).getVector(); 182 ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector(); 183 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 184 185 (*jvp)[0] = 2.0*(*xp)[0]*(*vp)[0] - (*vp)[1] - 2.0*(*xp)[3]*(*vp)[3]; 186 187 } 188 applyAdjointJacobian(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)189 void applyAdjointJacobian( Vector<Real> &ajv, const Vector<Real> &v, 190 const Vector<Real> &x, Real &tol ) { 191 192 ROL::Ptr<vector> ajvp = dynamic_cast<SV&>(ajv).getVector(); 193 ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector(); 194 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 195 196 (*ajvp)[0] = 2.0*(*xp)[0]*(*vp)[0]; 197 (*ajvp)[1] = -(*vp)[0]; 198 (*ajvp)[2] = 0.0; 199 (*ajvp)[3] = -2.0*(*vp)[0]*(*xp)[3]; 200 } 201 applyAdjointHessian(Vector<Real> & ahuv,const Vector<Real> & u,const Vector<Real> & v,const Vector<Real> & x,Real & tol)202 void applyAdjointHessian(Vector<Real> &ahuv, const Vector<Real> &u, 203 const Vector<Real> &v, const Vector<Real> &x, 204 Real &tol) { 205 206 ROL::Ptr<vector> ahuvp = dynamic_cast<SV&>(ahuv).getVector(); 207 ROL::Ptr<const vector> up = dynamic_cast<const SV&>(u).getVector(); 208 ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector(); 209 ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 210 211 // (*cp)[0] = std::pow((*xp)[0],2)-(*xp)[1]-std::pow((*xp)[3],2); 212 213 (*ahuvp)[0] = 2.0*(*up)[0]*(*vp)[0]; 214 (*ahuvp)[1] = 0.0; 215 (*ahuvp)[2] = 0.0; 216 (*ahuvp)[3] = -2.0*(*up)[0]*(*vp)[3]; 217 } 218 219 220 }; 221 222 template<class Real> 223 class getHS39 : public TestProblem<Real> { 224 public: getHS39(void)225 getHS39(void) {} 226 getObjective(void) const227 Ptr<Objective<Real>> getObjective(void) const { 228 // Instantiate Objective Function 229 return ROL::makePtr<Objective_HS39<Real>>(); 230 } 231 getInitialGuess(void) const232 Ptr<Vector<Real>> getInitialGuess(void) const { 233 // Problem dimension 234 int n = 4; 235 // Get Initial Guess 236 ROL::Ptr<std::vector<Real> > x0p = ROL::makePtr<std::vector<Real>>(n,2.0); 237 return ROL::makePtr<StdVector<Real>>(x0p); 238 } 239 getSolution(const int i=0) const240 Ptr<Vector<Real>> getSolution(const int i = 0) const { 241 // Problem dimension 242 int n = 4; 243 // Get Solution 244 ROL::Ptr<std::vector<Real> > xp = ROL::makePtr<std::vector<Real>>(n,0.0); 245 (*xp)[0] = static_cast<Real>(1); 246 (*xp)[1] = static_cast<Real>(1); 247 (*xp)[2] = static_cast<Real>(0); 248 (*xp)[3] = static_cast<Real>(0); 249 return ROL::makePtr<StdVector<Real>>(xp); 250 } 251 getEqualityConstraint(void) const252 Ptr<Constraint<Real>> getEqualityConstraint(void) const { 253 std::vector<Ptr<Constraint<Real>>> cvec(2); 254 cvec[0] = makePtr<Constraint_HS39a<Real>>(); 255 cvec[1] = makePtr<Constraint_HS39b<Real>>(); 256 return ROL::makePtr<Constraint_Partitioned<Real>>(cvec); 257 } 258 getEqualityMultiplier(void) const259 Ptr<Vector<Real>> getEqualityMultiplier(void) const { 260 std::vector<Ptr<Vector<Real>>> lvec(2); 261 lvec[0] = makePtr<StdVector<Real>>(makePtr<std::vector<Real>>(1,0.0)); 262 lvec[1] = makePtr<StdVector<Real>>(makePtr<std::vector<Real>>(1,0.0)); 263 return ROL::makePtr<PartitionedVector<Real>>(lvec); 264 } 265 }; 266 267 268 269 } // End ZOO Namespace 270 } // End ROL Namespace 271 272 #endif 273