EpetraMultiPointModelEval4DOpt.cpp

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

Generated on Tue Oct 20 12:45:30 2009 for EpetraExt by doxygen 1.4.7