EpetraExt Development
EpetraModelEval4DOpt.cpp
Go to the documentation of this file.
00001 #include "EpetraModelEval4DOpt.hpp"
00002 #include "Teuchos_ScalarTraits.hpp"
00003 #include "Epetra_SerialComm.h"
00004 #include "Epetra_CrsMatrix.h"
00005 
00006 #ifdef HAVE_MPI
00007 #  include "Epetra_MpiComm.h"
00008 #else
00009 #  include "Epetra_SerialComm.h"
00010 #endif
00011 
00012 namespace {
00013 
00014 inline double sqr( const double& s ) { return s*s; }
00015 
00016 } // namespace
00017 
00018 EpetraModelEval4DOpt::EpetraModelEval4DOpt(
00019   const double         xt0
00020   ,const double        xt1
00021   ,const double        pt0
00022   ,const double        pt1
00023   ,const double        d
00024   ,const double        x00
00025   ,const double        x01
00026   ,const double        p00
00027   ,const double        p01
00028   )
00029   :xt0_(xt0),xt1_(xt1),pt0_(pt0),pt1_(pt1),d_(d),
00030    isInitialized_(false),supportDerivs_(true)
00031 {
00032   using Teuchos::rcp;
00033 
00034   epetra_comm_ = rcp(new Epetra_SerialComm());
00035 
00036   const int nx = 2, np = 2, ng = 1;
00037 
00038   map_x_ = rcp(new Epetra_Map(nx,0,*epetra_comm_));
00039   map_p_ = rcp(new Epetra_Map(np,0,*epetra_comm_));
00040   map_g_ = rcp(new Epetra_Map(ng,0,*epetra_comm_));
00041 
00042   //const double inf = infiniteBound();
00043   const double inf = 1e+50;
00044 
00045   xL_ = rcp(new Epetra_Vector(*map_x_));  xL_->PutScalar(-inf);
00046   xU_ = rcp(new Epetra_Vector(*map_x_));  xU_->PutScalar(+inf);
00047   x0_ = rcp(new Epetra_Vector(*map_x_));  (*x0_)[0] = x00; (*x0_)[1] = x01;
00048   pL_ = rcp(new Epetra_Vector(*map_p_));  pL_->PutScalar(-inf);
00049   pU_ = rcp(new Epetra_Vector(*map_p_));  pU_->PutScalar(+inf);
00050   p0_ = rcp(new Epetra_Vector(*map_p_));  (*p0_)[0] = p00; (*p0_)[1] = p01;
00051   gL_ = rcp(new Epetra_Vector(*map_g_));  gL_->PutScalar(-inf);
00052   gU_ = rcp(new Epetra_Vector(*map_g_));  gU_->PutScalar(+inf);
00053 
00054   // Initialize the graph for W CrsMatrix object
00055   W_graph_ = rcp(new Epetra_CrsGraph(::Copy,*map_x_,nx));
00056   {
00057     int indices[nx] = { 0, 1 };
00058     for( int i = 0; i < nx; ++i )
00059       W_graph_->InsertGlobalIndices(i,nx,indices);
00060   }
00061   W_graph_->FillComplete();
00062 
00063   isInitialized_ = true;
00064 
00065 }
00066 
00067 
00068 void EpetraModelEval4DOpt::setSupportDerivs( bool supportDerivs )
00069 {
00070   supportDerivs_ = supportDerivs;
00071 }
00072 
00073 
00074 void EpetraModelEval4DOpt::set_p_bounds(
00075   double pL0, double pL1, double pU0, double pU1
00076   )
00077 {
00078   (*pL_)[0] = pL0;
00079   (*pL_)[1] = pL1;
00080   (*pU_)[0] = pU0;
00081   (*pU_)[1] = pU1;
00082 }
00083 
00084 void EpetraModelEval4DOpt::set_x_bounds(
00085   double xL0, double xL1, double xU0, double xU1
00086   )
00087 {
00088   (*xL_)[0] = xL0;
00089   (*xL_)[1] = xL1;
00090   (*xU_)[0] = xU0;
00091   (*xU_)[1] = xU1;
00092 }
00093 
00094 // Overridden from EpetraExt::ModelEvaluator
00095 
00096 Teuchos::RefCountPtr<const Epetra_Map>
00097 EpetraModelEval4DOpt::get_x_map() const
00098 {
00099   return map_x_;
00100 }
00101 
00102 Teuchos::RefCountPtr<const Epetra_Map>
00103 EpetraModelEval4DOpt::get_f_map() const
00104 {
00105   return map_x_;
00106 }
00107 
00108 Teuchos::RefCountPtr<const Epetra_Map>
00109 EpetraModelEval4DOpt::get_p_map(int l) const
00110 {
00111   TEST_FOR_EXCEPT(l!=0);
00112   return map_p_;
00113 }
00114 
00115 Teuchos::RefCountPtr<const Epetra_Map>
00116 EpetraModelEval4DOpt::get_g_map(int j) const
00117 {
00118   TEST_FOR_EXCEPT(j!=0);
00119   return map_g_;
00120 }
00121 
00122 Teuchos::RefCountPtr<const Epetra_Vector>
00123 EpetraModelEval4DOpt::get_x_init() const
00124 {
00125   return x0_;
00126 }
00127 
00128 Teuchos::RefCountPtr<const Epetra_Vector>
00129 EpetraModelEval4DOpt::get_p_init(int l) const
00130 {
00131   TEST_FOR_EXCEPT(l!=0);
00132   return p0_;
00133 }
00134 
00135 Teuchos::RefCountPtr<const Epetra_Vector>
00136 EpetraModelEval4DOpt::get_x_lower_bounds() const
00137 {
00138   return xL_;
00139 }
00140 
00141 Teuchos::RefCountPtr<const Epetra_Vector>
00142 EpetraModelEval4DOpt::get_x_upper_bounds() const
00143 {
00144   return xU_;
00145 }
00146 
00147 Teuchos::RefCountPtr<const Epetra_Vector>
00148 EpetraModelEval4DOpt::get_p_lower_bounds(int l) const
00149 {
00150   TEST_FOR_EXCEPT(l!=0);
00151   return pL_;
00152 }
00153 
00154 Teuchos::RefCountPtr<const Epetra_Vector>
00155 EpetraModelEval4DOpt::get_p_upper_bounds(int l) const
00156 {
00157   TEST_FOR_EXCEPT(l!=0);
00158   return pU_;
00159 }
00160 
00161 Teuchos::RefCountPtr<Epetra_Operator>
00162 EpetraModelEval4DOpt::create_W() const
00163 {
00164   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00165 }
00166 
00167 EpetraExt::ModelEvaluator::InArgs
00168 EpetraModelEval4DOpt::createInArgs() const
00169 {
00170   InArgsSetup inArgs;
00171   inArgs.setModelEvalDescription(this->description());
00172   inArgs.set_Np(1);
00173   inArgs.setSupports(IN_ARG_x,true);
00174   return inArgs;
00175 }
00176 
00177 EpetraExt::ModelEvaluator::OutArgs
00178 EpetraModelEval4DOpt::createOutArgs() const
00179 {
00180   OutArgsSetup outArgs;
00181   outArgs.setModelEvalDescription(this->description());
00182   outArgs.set_Np_Ng(1,1);
00183   outArgs.setSupports(OUT_ARG_f,true);
00184   outArgs.setSupports(OUT_ARG_W,true);
00185   outArgs.set_W_properties(
00186     DerivativeProperties(
00187       DERIV_LINEARITY_NONCONST
00188       ,DERIV_RANK_FULL
00189       ,true // supportsAdjoint
00190       )
00191     );
00192   if (supportDerivs_) {
00193     outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
00194     outArgs.set_DfDp_properties(
00195       0,DerivativeProperties(
00196         DERIV_LINEARITY_CONST
00197         ,DERIV_RANK_DEFICIENT
00198         ,true // supportsAdjoint
00199         )
00200       );
00201     outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
00202     outArgs.set_DgDx_properties(
00203       0,DerivativeProperties(
00204         DERIV_LINEARITY_NONCONST
00205         ,DERIV_RANK_DEFICIENT
00206         ,true // supportsAdjoint
00207         )
00208       );
00209     outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_TRANS_MV_BY_ROW);
00210     outArgs.set_DgDp_properties(
00211       0,0,DerivativeProperties(
00212         DERIV_LINEARITY_NONCONST
00213         ,DERIV_RANK_DEFICIENT
00214         ,true // supportsAdjoint
00215         )
00216       );
00217   }
00218   return outArgs;
00219 }
00220 
00221 
00222 void EpetraModelEval4DOpt::evalModel(
00223   const InArgs& inArgs, const OutArgs& outArgs
00224   ) const
00225 {
00226   using Teuchos::dyn_cast;
00227   using Teuchos::rcp_dynamic_cast;
00228   //
00229   // Get the input arguments
00230   //
00231   Teuchos::RefCountPtr<const Epetra_Vector> p_in = inArgs.get_p(0);
00232   const Epetra_Vector &p = (p_in.get() ? *p_in : *p0_);
00233   const Epetra_Vector &x = *inArgs.get_x();
00234   //
00235   // Get the output arguments
00236   //
00237   Epetra_Vector       *f_out = outArgs.get_f().get();
00238   Epetra_Vector       *g_out = outArgs.get_g(0).get();
00239   Epetra_Operator     *W_out = outArgs.get_W().get();
00240   Epetra_MultiVector  *DfDp_out = supportDerivs_ ? get_DfDp_mv(0,outArgs).get() : 0;
00241   Epetra_MultiVector  *DgDx_trans_out = supportDerivs_ ? get_DgDx_mv(0,outArgs,DERIV_TRANS_MV_BY_ROW).get() : 0;
00242   Epetra_MultiVector  *DgDp_trans_out = supportDerivs_ ? get_DgDp_mv(0,0,outArgs,DERIV_TRANS_MV_BY_ROW).get() : 0;
00243   //
00244   // Compute the functions
00245   //
00246   if(f_out) {
00247     Epetra_Vector &f = *f_out;
00248     // zero-based indexing for Epetra!
00249     f[0] =        x[0]      + x[1]*x[1] - p[0];
00250     f[1] = d_ * ( x[0]*x[0] - x[1]      - p[1] );
00251   }
00252   if(g_out) {
00253     Epetra_Vector &g = *g_out;
00254     g[0] = 0.5 * ( sqr(x[0]-xt0_) + sqr(x[1]-xt1_) + sqr(p[0]-pt0_) + sqr(p[1]-pt1_) );
00255   }
00256   if(W_out) {
00257     Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
00258     DfDx.PutScalar(0.0);
00259     //
00260     // Fill W = DfDx
00261     //
00262     // W = DfDx = [      1.0,  2*x[1] ]
00263     //            [ 2*d*x[0],     -d  ]
00264     //
00265     double values[2];
00266     int indexes[2];
00267     // Row [0]
00268     values[0] = 1.0;           indexes[0] = 0;
00269     values[1] = 2.0*x[1];      indexes[1] = 1;
00270     DfDx.SumIntoGlobalValues( 0, 2, values, indexes );
00271     // Row [1]
00272     values[0] = 2.0*d_*x[0];   indexes[0] = 0;
00273     values[1] = -d_;           indexes[1] = 1;
00274     DfDx.SumIntoGlobalValues( 1, 2, values, indexes );
00275   }
00276   if(DfDp_out) {
00277     Epetra_MultiVector &DfDp = *DfDp_out;
00278     DfDp.PutScalar(0.0);
00279     DfDp[0][0] = -1.0;
00280     DfDp[1][1] = -d_;
00281   }
00282   if(DgDx_trans_out) {
00283     Epetra_Vector &DgDx_trans = *(*DgDx_trans_out)(0);
00284     DgDx_trans[0] = x[0]-xt0_;
00285     DgDx_trans[1] = x[1]-xt1_;
00286   }
00287   if(DgDp_trans_out) {
00288     Epetra_Vector &DgDp_trans = *(*DgDp_trans_out)(0);
00289     DgDp_trans[0] = p[0]-pt0_;
00290     DgDp_trans[1] = p[1]-pt1_;
00291   }
00292 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines