EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
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
55namespace {
56
57inline 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 }
104 W_graph_->FillComplete();
105
106 isInitialized_ = true;
107
108}
109
110
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
139Teuchos::RCP<const Epetra_Map>
141{
142 return map_x_;
143}
144
145Teuchos::RCP<const Epetra_Map>
147{
148 return map_x_;
149}
150
151Teuchos::RCP<const Epetra_Map>
153{
154 TEUCHOS_TEST_FOR_EXCEPT(l!=0);
155 return map_p_;
156}
157
158Teuchos::RCP<const Epetra_Map>
160{
161 TEUCHOS_TEST_FOR_EXCEPT(j!=0);
162 return map_g_;
163}
164
165Teuchos::RCP<const Epetra_Vector>
167{
168 return x0_;
169}
170
171Teuchos::RCP<const Epetra_Vector>
173{
174 TEUCHOS_TEST_FOR_EXCEPT(l!=0);
175 return p0_;
176}
177
178Teuchos::RCP<const Epetra_Vector>
180{
181 return xL_;
182}
183
184Teuchos::RCP<const Epetra_Vector>
186{
187 return xU_;
188}
189
190Teuchos::RCP<const Epetra_Vector>
192{
193 TEUCHOS_TEST_FOR_EXCEPT(l!=0);
194 return pL_;
195}
196
197Teuchos::RCP<const Epetra_Vector>
199{
200 TEUCHOS_TEST_FOR_EXCEPT(l!=0);
201 return pU_;
202}
203
204Teuchos::RCP<Epetra_Operator>
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_) {
236 outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
237 outArgs.set_DfDp_properties(
241 ,true // supportsAdjoint
242 )
243 );
244 outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
245 outArgs.set_DgDx_properties(
249 ,true // supportsAdjoint
250 )
251 );
252 outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_TRANS_MV_BY_ROW);
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 setModelEvalDescription(const std::string &modelEvalDescription)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< Epetra_Vector > gL_
Teuchos::RCP< Epetra_Vector > pL_
Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
Teuchos::RCP< const Epetra_Comm > epetra_comm_
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Teuchos::RCP< const Epetra_Map > map_x_
Teuchos::RCP< const Epetra_Map > get_f_map() const
void set_x_bounds(double xL0, double xL1, double xU0, double xU1)
Teuchos::RCP< const Epetra_Vector > get_x_lower_bounds() const
Teuchos::RCP< const Epetra_Map > map_p_
Teuchos::RCP< Epetra_Vector > xL_
void set_p_bounds(double pL0, double pL1, double pU0, double pU1)
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const
\breif .
Teuchos::RCP< Epetra_Vector > pU_
Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< Epetra_Vector > xU_
Teuchos::RCP< Epetra_Vector > x0_
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
void setSupportDerivs(bool supportDerivs)
Teuchos::RCP< const Epetra_Map > map_g_
Teuchos::RCP< Epetra_Vector > p0_
Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
\breif .
Teuchos::RCP< Epetra_CrsGraph > W_graph_
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Teuchos::RCP< const Epetra_Map > get_x_map() const
Teuchos::RCP< Epetra_Vector > gU_
int PutScalar(double ScalarConstant)