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 45 #pragma once 46 #ifndef ROL_SEMISMOOTHNEWTONDUALMODEL_HPP 47 #define ROL_SEMISMOOTHNEWTONDUALMODEL_HPP 48 49 #include "ROL_TrustRegionModel.hpp" 50 #include "ROL_InactiveSetVector.hpp" 51 #include "ROL_BoundConstraint.hpp" 52 #include "ROL_VectorWorkspace.hpp" 53 54 /** @ingroup func_group 55 \class ROL::SemismoothNewtonDualModel 56 \brief Implements the dual variable model function for a 57 semismooth Newton step. 58 59 Reference: 60 Konstantin Pieper dissertation "Finite element discretization and efficient 61 numerical solution of elliptic and parabolic sparse control problems." 62 http://nbn-resolving.de/urn/resolver.pl?urn:nbn:de:bvb:91-diss-20150420-1241413-1-4 63 ----- 64 */ 65 66 67 namespace ROL { 68 69 70 template<class Real> 71 class SemismoothNewtonDualModel : public TrustRegionModel<Real> { 72 73 using V = Vector<Real>; 74 using VPrim = InactiveSet_PrimalVector<Real>; 75 using VDual = InactiveSet_DualVector<Real>; 76 77 using Obj = Objective<Real>; 78 using Sec = Secant<Real>; 79 using Bnd = BoundConstraint<Real>; 80 81 private: 82 83 class ProjectedObjective : public Objective<Real> { 84 private: 85 Obj& objPrimal_; 86 Bnd& bnd_; 87 Ptr<V> primalVec_; 88 89 public: ProjectedObjective(Obj & objPrimal,Bnd & bnd,const Ptr<V> & primalVec)90 ProjectedObjective( Obj& objPrimal, Bnd& bnd, const Ptr<V>& primalVec ) : 91 objPrimal_(objPrimal), bnd_(bnd), primalVec_( primalVec ) {} 92 value(const V & p,Real & tol)93 Real value( const V& p, Real& tol ) override { 94 primalVec_->set(p); 95 bnd_.project(*primalVec_); 96 return objPrimal_->value(*primalVec_, tol); 97 } 98 gradient(V & g,const V & p,Real & tol)99 void gradient( V& g, const V& p, Real& tol ) override { 100 primalVec_->set(p); 101 bnd_.project(*primalVec_); 102 objPrimal_->gradient(g,*primalVec_, tol); 103 } 104 hessVec(V & hv,const V & v,const V & p,Real & tol)105 void hessVec( V& hv, const V& v, const V& p, Real& tol ) override { 106 primalVec_->set(p); 107 bnd_.project(*primalVec_); 108 objPrimal_->hessVec(hv,v,*primalVec_, tol); 109 } 110 111 }; // ProjectedObjective 112 113 ProjectedObjective projObj_; 114 Bnd bnd_; 115 Sec secant_; 116 Ptr<V> p_, g_, x_; 117 Ptr<V> ones_; 118 Ptr<VPrim> s_; 119 Real alpha_; 120 121 VectorWorkspace<Real> workspace_; 122 123 124 public: 125 SemismoothNewtonDualModel(Obj & obj,Bnd & bnd,const V & p,const V & g,const Real alpha)126 SemismoothNewtonDualModel( Obj& obj, Bnd& bnd, const V& p, const V& g, const Real alpha ) : 127 TrustRegionModel( obj, p, g, false ), bnd_( bnd ), 128 p_( p.clone() ), g_( p.dual().clone() ), x_( p.clone() ), ones_( p.clone() ), 129 s_( p.clone(), ones_, p_, bnd_ ), projObj_( obj, bnd, p_ ), alpha_(alpha) { 130 131 ones_->setScalar( Real(1.0) ); 132 } 133 134 value(const V & s,Real & tol)135 Real value( const V& s, Real& tol ) { 136 137 auto hs = workspace_.clone(*g_); 138 139 gradient(*g_,s,tol); 140 hessVec(*hs,s,s,tol); 141 hs->scale( 0.5 ); 142 hs->plus(*g_); 143 s_->set(s); 144 return s_->dot(*hs); 145 } 146 gradient(V & g,const V & s,Real & tol)147 void gradient( V& g, const V& s, Real& tol ) { 148 projObj_->gradient(g,*p_,tol); 149 g.axpy(alpha_,*p_); 150 } 151 hessVec(V & hv,const V & v,const V & s,Real & tol)152 void hessVec( V& hv, const V& v, const V& s, Real& tol ) { 153 auto vprune_ = workspace_.copy(v); 154 bnd_->pruneActive( *vprune_, *p_ ); 155 projObj_->hessVec( hv, *vprune_, *p_, tol ); 156 hv.axpy(alpha_,v); 157 } 158 update(const V & p,bool flag=true,int iter=-1)159 void update( const V& p, bool flag = true, int iter = -1 ) { 160 p_->set(p); 161 auto x = this->getIterate(); 162 } 163 164 } // namespace ROL 165 166