EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraMultiPointModelEval4DOpt.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 
45 #include "Teuchos_ScalarTraits.hpp"
46 #include "Epetra_SerialComm.h"
47 #include "Epetra_CrsMatrix.h"
48 
49 #include "Epetra_MpiComm.h"
50 
51 namespace {
52 
53 inline double sqr( const double& s ) { return s*s; }
54 
55 } // namespace
56 
58  Teuchos::RCP<Epetra_Comm> epetra_comm
59  ,const double xt0
60  ,const double xt1
61  ,const double pt0
62  ,const double pt1
63  ,const double d
64  ,const double x00
65  ,const double x01
66  ,const double p00
67  ,const double p01
68  ,const double q0
69  )
70  :isInitialized_(false), epetra_comm_(epetra_comm),
71  xt0_(xt0),xt1_(xt1),pt0_(pt0),pt1_(pt1),d_(d)
72 {
73  using Teuchos::rcp;
74 
75  const int nx = 2, np = 2, ng = 1, nq = 1;
76 
77  map_x_ = rcp(new Epetra_Map(nx,0,*epetra_comm_));
78  map_p_ = rcp(new Epetra_Map(np,0,*epetra_comm_));
79  map_q_ = rcp(new Epetra_Map(nq,0,*epetra_comm_));
80  map_g_ = rcp(new Epetra_Map(ng,0,*epetra_comm_));
81 
82  //const double inf = infiniteBound();
83  const double inf = 1e+50;
84 
85  xL_ = rcp(new Epetra_Vector(*map_x_)); xL_->PutScalar(-inf);
86  xU_ = rcp(new Epetra_Vector(*map_x_)); xU_->PutScalar(+inf);
87  x0_ = rcp(new Epetra_Vector(*map_x_)); (*x0_)[0] = x00; (*x0_)[1] = x01;
88  pL_ = rcp(new Epetra_Vector(*map_p_)); pL_->PutScalar(-inf);
89  pU_ = rcp(new Epetra_Vector(*map_p_)); pU_->PutScalar(+inf);
90  p0_ = rcp(new Epetra_Vector(*map_p_)); (*p0_)[0] = p00; (*p0_)[1] = p01;
91  gL_ = rcp(new Epetra_Vector(*map_g_)); gL_->PutScalar(-inf);
92  gU_ = rcp(new Epetra_Vector(*map_g_)); gU_->PutScalar(+inf);
93  q_ = rcp(new Epetra_Vector(*map_q_)); (*q_)[0] = q0;
94  qL_ = rcp(new Epetra_Vector(*map_q_)); (*qL_)[0] = -inf;
95  qU_ = rcp(new Epetra_Vector(*map_q_)); (*qU_)[0] = 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 
111  double pL0, double pL1, double pU0, double pU1
112  )
113 {
114  (*pL_)[0] = pL0;
115  (*pL_)[1] = pL1;
116  (*pU_)[0] = pU0;
117  (*pU_)[1] = pU1;
118 }
119 
121  double xL0, double xL1, double xU0, double xU1
122  )
123 {
124  (*xL_)[0] = xL0;
125  (*xL_)[1] = xL1;
126  (*xU_)[0] = xU0;
127  (*xU_)[1] = xU1;
128 }
129 
130 // Overridden from EpetraExt::ModelEvaluator
131 
134 {
135  return map_x_;
136 }
137 
140 {
141  return map_x_;
142 }
143 
146 {
148  if (l==0) return map_p_;
149  else return map_q_;
150 }
151 
154 {
156  return map_g_;
157 }
158 
161 {
162  return x0_;
163 }
164 
167 {
169  if (l==0) return p0_;
170  else return q_;
171 }
172 
175 {
176  return xL_;
177 }
178 
181 {
182  return xU_;
183 }
184 
187 {
189  if (l==0) return pL_;
190  else return qL_;
191 }
192 
195 {
197  if (l==0) return pU_;
198  else return qU_;
199 }
200 
203 {
204  return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
205 }
206 
209 {
210  InArgsSetup inArgs;
211  inArgs.setModelEvalDescription(this->description());
212  inArgs.set_Np(2);
213  inArgs.setSupports(IN_ARG_x,true);
214  return inArgs;
215 }
216 
219 {
220  OutArgsSetup outArgs;
221  outArgs.setModelEvalDescription(this->description());
222  outArgs.set_Np_Ng(2,1);
223  outArgs.setSupports(OUT_ARG_f,true);
224  outArgs.setSupports(OUT_ARG_W,true);
225  outArgs.set_W_properties(
229  ,true // supportsAdjoint
230  )
231  );
233  outArgs.set_DfDp_properties(
237  ,true // supportsAdjoint
238  )
239  );
241  outArgs.set_DgDx_properties(
245  ,true // supportsAdjoint
246  )
247  );
249  outArgs.set_DgDp_properties(
253  ,true // supportsAdjoint
254  )
255  );
256  return outArgs;
257 }
258 
259 void EpetraMultiPointModelEval4DOpt::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
260 {
261  using Teuchos::dyn_cast;
262  using Teuchos::rcp_dynamic_cast;
263  //
264  // Get the input arguments
265  //
266  Teuchos::RCP<const Epetra_Vector> p_in = inArgs.get_p(0);
267  Teuchos::RCP<const Epetra_Vector> q_in = inArgs.get_p(1);
268  const Epetra_Vector &p = (p_in.get() ? *p_in : *p0_);
269  const Epetra_Vector &q = (q_in.get() ? *q_in : *q_);
270  const Epetra_Vector &x = *inArgs.get_x();
271  //
272  // Get the output arguments
273  //
274  Epetra_Vector *f_out = outArgs.get_f().get();
275  Epetra_Vector *g_out = outArgs.get_g(0).get();
276  Epetra_Operator *W_out = outArgs.get_W().get();
277  Epetra_MultiVector *DfDp_out = get_DfDp_mv(0,outArgs).get();
278  Epetra_MultiVector *DgDx_trans_out = get_DgDx_mv(0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
279  Epetra_MultiVector *DgDp_trans_out = get_DgDp_mv(0,0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
280  //
281  // Compute the functions
282  //
283  if(f_out) {
284  Epetra_Vector &f = *f_out;
285  // zero-based indexing for Epetra!
286  //q[0] added for multipoint problem!
287  f[0] = x[0] + x[1]*x[1] - p[0] + q[0];
288  f[1] = d_ * ( x[0]*x[0] - x[1] - p[1] );
289  }
290  if(g_out) {
291  Epetra_Vector &g = *g_out;
292  g[0] = 0.5 * ( sqr(x[0]-xt0_) + sqr(x[1]-xt1_) + sqr(p[0]-pt0_) + sqr(p[1]-pt1_) );
293  }
294  if(W_out) {
295  Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
296  DfDx.PutScalar(0.0);
297  //
298  // Fill W = DfDx
299  //
300  // W = DfDx = [ 1.0, 2*x[1] ]
301  // [ 2*d*x[0], -d ]
302  //
303  double values[2];
304  int indexes[2];
305  // Row [0]
306  values[0] = 1.0; indexes[0] = 0;
307  values[1] = 2.0*x[1]; indexes[1] = 1;
308  DfDx.SumIntoGlobalValues( 0, 2, values, indexes );
309  // Row [1]
310  values[0] = 2.0*d_*x[0]; indexes[0] = 0;
311  values[1] = -d_; indexes[1] = 1;
312  DfDx.SumIntoGlobalValues( 1, 2, values, indexes );
313  }
314  if(DfDp_out) {
315  Epetra_MultiVector &DfDp = *DfDp_out;
316  DfDp.PutScalar(0.0);
317  DfDp[0][0] = -1.0;
318  DfDp[1][1] = -d_;
319  }
320  if(DgDx_trans_out) {
321  Epetra_Vector &DgDx_trans = *(*DgDx_trans_out)(0);
322  DgDx_trans[0] = x[0]-xt0_;
323  DgDx_trans[1] = x[1]-xt1_;
324  }
325  if(DgDp_trans_out) {
326  Epetra_Vector &DgDp_trans = *(*DgDp_trans_out)(0);
327  DgDp_trans[0] = p[0]-pt0_;
328  DgDp_trans[1] = p[1]-pt1_;
329  }
330 }
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< const Epetra_Vector > get_p_init(int l) const
Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
Teuchos::RCP< const Epetra_Map > get_x_map() const
Teuchos::RCP< Epetra_MultiVector > get_DgDp_mv(const int j, const int l, const ModelEvaluator::OutArgs &outArgs, const ModelEvaluator::EDerivativeMultiVectorOrientation mvOrientation)
Teuchos::RCP< Epetra_Operator > get_W() const
Teuchos::RCP< const Epetra_Map > map_p_
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
void setSupports(EInArgsMembers arg, bool supports=true)
Teuchos::RCP< Epetra_Operator > create_W() const
T * get() const
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
.
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
Teuchos::RCP< const Epetra_Map > get_f_map() const
Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
.
int PutScalar(double ScalarConstant)
void set_DgDp_properties(int j, int l, const DerivativeProperties &properties)
Teuchos::RCP< const Epetra_Map > map_g_
void set_DfDp_properties(int l, const DerivativeProperties &properties)
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)
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()
void setModelEvalDescription(const std::string &modelEvalDescription)
void set_p_bounds(double pL0, double pL1, double pU0, double pU1)
Teuchos::RCP< const Epetra_Map > map_x_
Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
EpetraMultiPointModelEval4DOpt(Teuchos::RCP< Epetra_Comm > epetra_comm, const double xt0=1.0, const double xt1=1.0, const double pt0=2.0, const double pt1=0.0, const double d=10.0, const double x00=1.0, const double x01=1.0, const double p00=2.0, const double p01=0.0, const double q0=0.0)
Teuchos::RCP< const Epetra_Vector > get_x_lower_bounds() const
Teuchos::RCP< const Epetra_Comm > epetra_comm_
void set_DgDx_properties(int j, const DerivativeProperties &properties)
void set_x_bounds(double xL0, double xL1, double xU0, double xU1)
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Teuchos::RCP< Epetra_CrsGraph > W_graph_
Evaluation< Epetra_Vector > get_f() const
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
Teuchos::RCP< const Epetra_Map > map_q_
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)