EpetraExt Development
Loading...
Searching...
No Matches
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
51namespace {
52
53inline 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 }
104 W_graph_->FillComplete();
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
132Teuchos::RCP<const Epetra_Map>
134{
135 return map_x_;
136}
137
138Teuchos::RCP<const Epetra_Map>
140{
141 return map_x_;
142}
143
144Teuchos::RCP<const Epetra_Map>
146{
147 TEUCHOS_TEST_FOR_EXCEPT(l>1);
148 if (l==0) return map_p_;
149 else return map_q_;
150}
151
152Teuchos::RCP<const Epetra_Map>
154{
155 TEUCHOS_TEST_FOR_EXCEPT(j!=0);
156 return map_g_;
157}
158
159Teuchos::RCP<const Epetra_Vector>
161{
162 return x0_;
163}
164
165Teuchos::RCP<const Epetra_Vector>
167{
168 TEUCHOS_TEST_FOR_EXCEPT(l>1);
169 if (l==0) return p0_;
170 else return q_;
171}
172
173Teuchos::RCP<const Epetra_Vector>
175{
176 return xL_;
177}
178
179Teuchos::RCP<const Epetra_Vector>
181{
182 return xU_;
183}
184
185Teuchos::RCP<const Epetra_Vector>
187{
188 TEUCHOS_TEST_FOR_EXCEPT(l>1);
189 if (l==0) return pL_;
190 else return qL_;
191}
192
193Teuchos::RCP<const Epetra_Vector>
195{
196 TEUCHOS_TEST_FOR_EXCEPT(l>1);
197 if (l==0) return pU_;
198 else return qU_;
199}
200
201Teuchos::RCP<Epetra_Operator>
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
259void 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 setSupports(EInArgsMembers arg, bool supports=true)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
void setModelEvalDescription(const std::string &modelEvalDescription)
void set_W_properties(const DerivativeProperties &properties)
void set_DfDp_properties(int l, const DerivativeProperties &properties)
void set_DgDp_properties(int j, int l, const DerivativeProperties &properties)
void set_DgDx_properties(int j, const DerivativeProperties &properties)
void setSupports(EOutArgsMembers arg, bool supports=true)
Teuchos::RCP< Epetra_Operator > get_W() const
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 <= j && j < this->Ng().
Evaluation< Epetra_Vector > get_f() const
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
Teuchos::RCP< const Epetra_Map > get_x_map() const
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) 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_Map > get_p_map(int l) const
\breif .
void set_x_bounds(double xL0, double xL1, double xU0, double xU1)
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
\breif .
void set_p_bounds(double pL0, double pL1, double pU0, double pU1)
Teuchos::RCP< const Epetra_Vector > get_x_lower_bounds() const
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int PutScalar(double ScalarConstant)
int PutScalar(double ScalarConstant)