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

Generated on Thu Sep 18 12:31:44 2008 for EpetraExt by doxygen 1.3.9.1