EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraModelEval4DOpt.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // EpetraExt: Epetra Extended - Linear Algebra Services Package
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #include "EpetraModelEval4DOpt.hpp"
45 #include "Teuchos_ScalarTraits.hpp"
46 #include "Epetra_SerialComm.h"
47 #include "Epetra_CrsMatrix.h"
48 
49 #ifdef HAVE_MPI
50 # include "Epetra_MpiComm.h"
51 #else
52 # include "Epetra_SerialComm.h"
53 #endif
54 
55 namespace {
56 
57 inline double sqr( const double& s ) { return s*s; }
58 
59 } // namespace
60 
62  const double xt0
63  ,const double xt1
64  ,const double pt0
65  ,const double pt1
66  ,const double d
67  ,const double x00
68  ,const double x01
69  ,const double p00
70  ,const double p01
71  )
72  :xt0_(xt0),xt1_(xt1),pt0_(pt0),pt1_(pt1),d_(d),
73  isInitialized_(false),supportDerivs_(true)
74 {
75  using Teuchos::rcp;
76 
78 
79  const int nx = 2, np = 2, ng = 1;
80 
81  map_x_ = rcp(new Epetra_Map(nx,0,*epetra_comm_));
82  map_p_ = rcp(new Epetra_Map(np,0,*epetra_comm_));
83  map_g_ = rcp(new Epetra_Map(ng,0,*epetra_comm_));
84 
85  //const double inf = infiniteBound();
86  const double inf = 1e+50;
87 
88  xL_ = rcp(new Epetra_Vector(*map_x_)); xL_->PutScalar(-inf);
89  xU_ = rcp(new Epetra_Vector(*map_x_)); xU_->PutScalar(+inf);
90  x0_ = rcp(new Epetra_Vector(*map_x_)); (*x0_)[0] = x00; (*x0_)[1] = x01;
91  pL_ = rcp(new Epetra_Vector(*map_p_)); pL_->PutScalar(-inf);
92  pU_ = rcp(new Epetra_Vector(*map_p_)); pU_->PutScalar(+inf);
93  p0_ = rcp(new Epetra_Vector(*map_p_)); (*p0_)[0] = p00; (*p0_)[1] = p01;
94  gL_ = rcp(new Epetra_Vector(*map_g_)); gL_->PutScalar(-inf);
95  gU_ = rcp(new Epetra_Vector(*map_g_)); gU_->PutScalar(+inf);
96 
97  // Initialize the graph for W CrsMatrix object
98  W_graph_ = rcp(new Epetra_CrsGraph(::Copy,*map_x_,nx));
99  {
100  int indices[nx] = { 0, 1 };
101  for( int i = 0; i < nx; ++i )
102  W_graph_->InsertGlobalIndices(i,nx,indices);
103  }
105 
106  isInitialized_ = true;
107 
108 }
109 
110 
111 void EpetraModelEval4DOpt::setSupportDerivs( bool supportDerivs )
112 {
113  supportDerivs_ = supportDerivs;
114 }
115 
116 
118  double pL0, double pL1, double pU0, double pU1
119  )
120 {
121  (*pL_)[0] = pL0;
122  (*pL_)[1] = pL1;
123  (*pU_)[0] = pU0;
124  (*pU_)[1] = pU1;
125 }
126 
128  double xL0, double xL1, double xU0, double xU1
129  )
130 {
131  (*xL_)[0] = xL0;
132  (*xL_)[1] = xL1;
133  (*xU_)[0] = xU0;
134  (*xU_)[1] = xU1;
135 }
136 
137 // Overridden from EpetraExt::ModelEvaluator
138 
141 {
142  return map_x_;
143 }
144 
147 {
148  return map_x_;
149 }
150 
153 {
155  return map_p_;
156 }
157 
160 {
162  return map_g_;
163 }
164 
167 {
168  return x0_;
169 }
170 
173 {
175  return p0_;
176 }
177 
180 {
181  return xL_;
182 }
183 
186 {
187  return xU_;
188 }
189 
192 {
194  return pL_;
195 }
196 
199 {
201  return pU_;
202 }
203 
206 {
207  return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
208 }
209 
212 {
213  InArgsSetup inArgs;
214  inArgs.setModelEvalDescription(this->description());
215  inArgs.set_Np(1);
216  inArgs.setSupports(IN_ARG_x,true);
217  return inArgs;
218 }
219 
222 {
223  OutArgsSetup outArgs;
224  outArgs.setModelEvalDescription(this->description());
225  outArgs.set_Np_Ng(1,1);
226  outArgs.setSupports(OUT_ARG_f,true);
227  outArgs.setSupports(OUT_ARG_W,true);
228  outArgs.set_W_properties(
232  ,true // supportsAdjoint
233  )
234  );
235  if (supportDerivs_) {
237  outArgs.set_DfDp_properties(
241  ,true // supportsAdjoint
242  )
243  );
245  outArgs.set_DgDx_properties(
249  ,true // supportsAdjoint
250  )
251  );
253  outArgs.set_DgDp_properties(
257  ,true // supportsAdjoint
258  )
259  );
260  }
261  return outArgs;
262 }
263 
264 
266  const InArgs& inArgs, const OutArgs& outArgs
267  ) const
268 {
269  using Teuchos::dyn_cast;
270  using Teuchos::rcp_dynamic_cast;
271  //
272  // Get the input arguments
273  //
274  Teuchos::RCP<const Epetra_Vector> p_in = inArgs.get_p(0);
275  const Epetra_Vector &p = (p_in.get() ? *p_in : *p0_);
276  const Epetra_Vector &x = *inArgs.get_x();
277  //
278  // Get the output arguments
279  //
280  Epetra_Vector *f_out = outArgs.get_f().get();
281  Epetra_Vector *g_out = outArgs.get_g(0).get();
282  Epetra_Operator *W_out = outArgs.get_W().get();
283  Epetra_MultiVector *DfDp_out = supportDerivs_ ? get_DfDp_mv(0,outArgs).get() : 0;
284  Epetra_MultiVector *DgDx_trans_out = supportDerivs_ ? get_DgDx_mv(0,outArgs,DERIV_TRANS_MV_BY_ROW).get() : 0;
285  Epetra_MultiVector *DgDp_trans_out = supportDerivs_ ? get_DgDp_mv(0,0,outArgs,DERIV_TRANS_MV_BY_ROW).get() : 0;
286  //
287  // Compute the functions
288  //
289  if(f_out) {
290  Epetra_Vector &f = *f_out;
291  // zero-based indexing for Epetra!
292  f[0] = x[0] + x[1]*x[1] - p[0];
293  f[1] = d_ * ( x[0]*x[0] - x[1] - p[1] );
294  }
295  if(g_out) {
296  Epetra_Vector &g = *g_out;
297  g[0] = 0.5 * ( sqr(x[0]-xt0_) + sqr(x[1]-xt1_) + sqr(p[0]-pt0_) + sqr(p[1]-pt1_) );
298  }
299  if(W_out) {
300  Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
301  DfDx.PutScalar(0.0);
302  //
303  // Fill W = DfDx
304  //
305  // W = DfDx = [ 1.0, 2*x[1] ]
306  // [ 2*d*x[0], -d ]
307  //
308  double values[2];
309  int indexes[2];
310  // Row [0]
311  values[0] = 1.0; indexes[0] = 0;
312  values[1] = 2.0*x[1]; indexes[1] = 1;
313  DfDx.SumIntoGlobalValues( 0, 2, values, indexes );
314  // Row [1]
315  values[0] = 2.0*d_*x[0]; indexes[0] = 0;
316  values[1] = -d_; indexes[1] = 1;
317  DfDx.SumIntoGlobalValues( 1, 2, values, indexes );
318  }
319  if(DfDp_out) {
320  Epetra_MultiVector &DfDp = *DfDp_out;
321  DfDp.PutScalar(0.0);
322  DfDp[0][0] = -1.0;
323  DfDp[1][1] = -d_;
324  }
325  if(DgDx_trans_out) {
326  Epetra_Vector &DgDx_trans = *(*DgDx_trans_out)(0);
327  DgDx_trans[0] = x[0]-xt0_;
328  DgDx_trans[1] = x[1]-xt1_;
329  }
330  if(DgDp_trans_out) {
331  Epetra_Vector &DgDp_trans = *(*DgDp_trans_out)(0);
332  DgDp_trans[0] = p[0]-pt0_;
333  DgDp_trans[1] = p[1]-pt1_;
334  }
335 }
void set_x_bounds(double xL0, double xL1, double xU0, double xU1)
Teuchos::RCP< const Epetra_Map > get_x_map() const
Teuchos::RCP< const Epetra_Map > map_p_
void set_W_properties(const DerivativeProperties &properties)
void f()
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 &lt;= j &amp;&amp; j &lt; this-&gt;Ng().
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< Epetra_Vector > p0_
Teuchos::RCP< Epetra_Vector > pL_
Teuchos::RCP< const Epetra_Map > get_f_map() const
Teuchos::RCP< Epetra_Vector > xL_
Teuchos::RCP< Epetra_CrsGraph > W_graph_
Teuchos::RCP< Epetra_MultiVector > get_DgDp_mv(const int j, const int l, const ModelEvaluator::OutArgs &outArgs, const ModelEvaluator::EDerivativeMultiVectorOrientation mvOrientation)
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Teuchos::RCP< Epetra_Operator > get_W() const
Teuchos::RCP< const Epetra_Map > map_x_
Teuchos::RCP< const Epetra_Map > map_g_
Teuchos::RCP< Epetra_Vector > xU_
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
void setSupports(EInArgsMembers arg, bool supports=true)
void set_p_bounds(double pL0, double pL1, double pU0, double pU1)
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
.
T * get() const
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
OutArgs createOutArgs() const
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
Teuchos::RCP< Epetra_Vector > gU_
Teuchos::RCP< const Epetra_Vector > get_x_init() const
int PutScalar(double ScalarConstant)
void set_DgDp_properties(int j, int l, const DerivativeProperties &properties)
Teuchos::RCP< const Epetra_Comm > epetra_comm_
void set_DfDp_properties(int l, const DerivativeProperties &properties)
Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
void setModelEvalDescription(const std::string &modelEvalDescription)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< Epetra_MultiVector > get_DfDp_mv(const int l, const ModelEvaluator::OutArgs &outArgs)
Teuchos::RCP< Epetra_Vector > gL_
T_To & dyn_cast(T_From &from)
Teuchos::RCP< Epetra_MultiVector > get_DgDx_mv(const int j, const ModelEvaluator::OutArgs &outArgs, const ModelEvaluator::EDerivativeMultiVectorOrientation mvOrientation)
virtual std::string description() const
void g()
Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< Epetra_Vector > pU_
Teuchos::RCP< Epetra_Vector > x0_
void setSupportDerivs(bool supportDerivs)
Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
.
void set_DgDx_properties(int j, const DerivativeProperties &properties)
Teuchos::RCP< Epetra_Operator > create_W() const
Evaluation< Epetra_Vector > get_f() const
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
Teuchos::RCP< const Epetra_Vector > get_x_lower_bounds() const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)