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

Generated on Wed May 12 21:24:46 2010 for EpetraExt by  doxygen 1.4.7