EpetraExt Development
EpetraMultiPointModelEval4DOpt.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00006 //                 Copyright (2011) Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "EpetraMultiPointModelEval4DOpt.hpp"
00045 #include "Teuchos_ScalarTraits.hpp"
00046 #include "Epetra_SerialComm.h"
00047 #include "Epetra_CrsMatrix.h"
00048 
00049 #include "Epetra_MpiComm.h"
00050 
00051 namespace {
00052 
00053 inline double sqr( const double& s ) { return s*s; }
00054 
00055 } // namespace
00056 
00057 EpetraMultiPointModelEval4DOpt::EpetraMultiPointModelEval4DOpt(
00058   Teuchos::RefCountPtr<Epetra_Comm> epetra_comm
00059   ,const double         xt0
00060   ,const double        xt1
00061   ,const double        pt0
00062   ,const double        pt1
00063   ,const double        d
00064   ,const double        x00
00065   ,const double        x01
00066   ,const double        p00
00067   ,const double        p01
00068   ,const double        q0
00069   )
00070   :isInitialized_(false), epetra_comm_(epetra_comm),
00071    xt0_(xt0),xt1_(xt1),pt0_(pt0),pt1_(pt1),d_(d)
00072 {
00073   using Teuchos::rcp;
00074 
00075   const int nx = 2, np = 2, ng = 1, nq = 1;
00076 
00077   map_x_ = rcp(new Epetra_Map(nx,0,*epetra_comm_));
00078   map_p_ = rcp(new Epetra_Map(np,0,*epetra_comm_));
00079   map_q_ = rcp(new Epetra_Map(nq,0,*epetra_comm_));
00080   map_g_ = rcp(new Epetra_Map(ng,0,*epetra_comm_));
00081 
00082   //const double inf = infiniteBound();
00083   const double inf = 1e+50;
00084 
00085   xL_ = rcp(new Epetra_Vector(*map_x_));  xL_->PutScalar(-inf);
00086   xU_ = rcp(new Epetra_Vector(*map_x_));  xU_->PutScalar(+inf);
00087   x0_ = rcp(new Epetra_Vector(*map_x_));  (*x0_)[0] = x00; (*x0_)[1] = x01;
00088   pL_ = rcp(new Epetra_Vector(*map_p_));  pL_->PutScalar(-inf);
00089   pU_ = rcp(new Epetra_Vector(*map_p_));  pU_->PutScalar(+inf);
00090   p0_ = rcp(new Epetra_Vector(*map_p_));  (*p0_)[0] = p00; (*p0_)[1] = p01;
00091   gL_ = rcp(new Epetra_Vector(*map_g_));  gL_->PutScalar(-inf);
00092   gU_ = rcp(new Epetra_Vector(*map_g_));  gU_->PutScalar(+inf);
00093   q_  = rcp(new Epetra_Vector(*map_q_));  (*q_)[0] = q0;
00094   qL_ = rcp(new Epetra_Vector(*map_q_));  (*qL_)[0] = -inf;
00095   qU_ = rcp(new Epetra_Vector(*map_q_));  (*qU_)[0] =  inf;
00096 
00097   // Initialize the graph for W CrsMatrix object
00098   W_graph_ = rcp(new Epetra_CrsGraph(::Copy,*map_x_,nx));
00099   {
00100     int indices[nx] = { 0, 1 };
00101     for( int i = 0; i < nx; ++i )
00102       W_graph_->InsertGlobalIndices(i,nx,indices);
00103   }
00104   W_graph_->FillComplete();
00105 
00106   isInitialized_ = true;
00107 
00108 }
00109 
00110 void EpetraMultiPointModelEval4DOpt::set_p_bounds(
00111   double pL0, double pL1, double pU0, double pU1
00112   )
00113 {
00114   (*pL_)[0] = pL0;
00115   (*pL_)[1] = pL1;
00116   (*pU_)[0] = pU0;
00117   (*pU_)[1] = pU1;
00118 }
00119 
00120 void EpetraMultiPointModelEval4DOpt::set_x_bounds(
00121   double xL0, double xL1, double xU0, double xU1
00122   )
00123 {
00124   (*xL_)[0] = xL0;
00125   (*xL_)[1] = xL1;
00126   (*xU_)[0] = xU0;
00127   (*xU_)[1] = xU1;
00128 }
00129 
00130 // Overridden from EpetraExt::ModelEvaluator
00131 
00132 Teuchos::RefCountPtr<const Epetra_Map>
00133 EpetraMultiPointModelEval4DOpt::get_x_map() const
00134 {
00135   return map_x_;
00136 }
00137 
00138 Teuchos::RefCountPtr<const Epetra_Map>
00139 EpetraMultiPointModelEval4DOpt::get_f_map() const
00140 {
00141   return map_x_;
00142 }
00143 
00144 Teuchos::RefCountPtr<const Epetra_Map>
00145 EpetraMultiPointModelEval4DOpt::get_p_map(int l) const
00146 {
00147   TEST_FOR_EXCEPT(l>1);
00148   if (l==0) return map_p_;
00149   else return map_q_;
00150 }
00151 
00152 Teuchos::RefCountPtr<const Epetra_Map>
00153 EpetraMultiPointModelEval4DOpt::get_g_map(int j) const
00154 {
00155   TEST_FOR_EXCEPT(j!=0);
00156   return map_g_;
00157 }
00158 
00159 Teuchos::RefCountPtr<const Epetra_Vector>
00160 EpetraMultiPointModelEval4DOpt::get_x_init() const
00161 {
00162   return x0_;
00163 }
00164 
00165 Teuchos::RefCountPtr<const Epetra_Vector>
00166 EpetraMultiPointModelEval4DOpt::get_p_init(int l) const
00167 {
00168   TEST_FOR_EXCEPT(l>1);
00169   if (l==0) return p0_;
00170   else return q_;
00171 }
00172 
00173 Teuchos::RefCountPtr<const Epetra_Vector>
00174 EpetraMultiPointModelEval4DOpt::get_x_lower_bounds() const
00175 {
00176   return xL_;
00177 }
00178 
00179 Teuchos::RefCountPtr<const Epetra_Vector>
00180 EpetraMultiPointModelEval4DOpt::get_x_upper_bounds() const
00181 {
00182   return xU_;
00183 }
00184 
00185 Teuchos::RefCountPtr<const Epetra_Vector>
00186 EpetraMultiPointModelEval4DOpt::get_p_lower_bounds(int l) const
00187 {
00188   TEST_FOR_EXCEPT(l>1);
00189   if (l==0) return pL_;
00190   else      return qL_;
00191 }
00192 
00193 Teuchos::RefCountPtr<const Epetra_Vector>
00194 EpetraMultiPointModelEval4DOpt::get_p_upper_bounds(int l) const
00195 {
00196   TEST_FOR_EXCEPT(l>1);
00197   if (l==0) return pU_;
00198   else      return qU_;
00199 }
00200 
00201 Teuchos::RefCountPtr<Epetra_Operator>
00202 EpetraMultiPointModelEval4DOpt::create_W() const
00203 {
00204   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00205 }
00206 
00207 EpetraExt::ModelEvaluator::InArgs
00208 EpetraMultiPointModelEval4DOpt::createInArgs() const
00209 {
00210   InArgsSetup inArgs;
00211   inArgs.setModelEvalDescription(this->description());
00212   inArgs.set_Np(2);
00213   inArgs.setSupports(IN_ARG_x,true);
00214   return inArgs;
00215 }
00216 
00217 EpetraExt::ModelEvaluator::OutArgs
00218 EpetraMultiPointModelEval4DOpt::createOutArgs() const
00219 {
00220   OutArgsSetup outArgs;
00221   outArgs.setModelEvalDescription(this->description());
00222   outArgs.set_Np_Ng(2,1);
00223   outArgs.setSupports(OUT_ARG_f,true);
00224   outArgs.setSupports(OUT_ARG_W,true);
00225   outArgs.set_W_properties(
00226     DerivativeProperties(
00227       DERIV_LINEARITY_NONCONST
00228       ,DERIV_RANK_FULL
00229       ,true // supportsAdjoint
00230       )
00231     );
00232   outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
00233   outArgs.set_DfDp_properties(
00234     0,DerivativeProperties(
00235       DERIV_LINEARITY_CONST
00236       ,DERIV_RANK_DEFICIENT
00237       ,true // supportsAdjoint
00238       )
00239     );
00240   outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
00241   outArgs.set_DgDx_properties(
00242     0,DerivativeProperties(
00243       DERIV_LINEARITY_NONCONST
00244       ,DERIV_RANK_DEFICIENT
00245       ,true // supportsAdjoint
00246       )
00247     );
00248   outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_TRANS_MV_BY_ROW);
00249   outArgs.set_DgDp_properties(
00250     0,0,DerivativeProperties(
00251       DERIV_LINEARITY_NONCONST
00252       ,DERIV_RANK_DEFICIENT
00253       ,true // supportsAdjoint
00254       )
00255     );
00256   return outArgs;
00257 }
00258 
00259 void EpetraMultiPointModelEval4DOpt::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
00260 {
00261   using Teuchos::dyn_cast;
00262   using Teuchos::rcp_dynamic_cast;
00263   //
00264   // Get the input arguments
00265   //
00266   Teuchos::RefCountPtr<const Epetra_Vector> p_in = inArgs.get_p(0);
00267   Teuchos::RefCountPtr<const Epetra_Vector> q_in = inArgs.get_p(1);
00268   const Epetra_Vector &p = (p_in.get() ? *p_in : *p0_);
00269   const Epetra_Vector &q = (q_in.get() ? *q_in : *q_);
00270   const Epetra_Vector &x = *inArgs.get_x();
00271   //
00272   // Get the output arguments
00273   //
00274   Epetra_Vector       *f_out = outArgs.get_f().get();
00275   Epetra_Vector       *g_out = outArgs.get_g(0).get();
00276   Epetra_Operator     *W_out = outArgs.get_W().get();
00277   Epetra_MultiVector  *DfDp_out = get_DfDp_mv(0,outArgs).get();
00278   Epetra_MultiVector  *DgDx_trans_out = get_DgDx_mv(0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
00279   Epetra_MultiVector  *DgDp_trans_out = get_DgDp_mv(0,0,outArgs,DERIV_TRANS_MV_BY_ROW).get();
00280   //
00281   // Compute the functions
00282   //
00283   if(f_out) {
00284     Epetra_Vector &f = *f_out;
00285     // zero-based indexing for Epetra!
00286     //q[0] added for multipoint problem!
00287     f[0] =        x[0]      + x[1]*x[1] - p[0] + q[0];
00288     f[1] = d_ * ( x[0]*x[0] - x[1]      - p[1] );
00289   }
00290   if(g_out) {
00291     Epetra_Vector &g = *g_out;
00292     g[0] = 0.5 * ( sqr(x[0]-xt0_) + sqr(x[1]-xt1_) + sqr(p[0]-pt0_) + sqr(p[1]-pt1_) );
00293   }
00294   if(W_out) {
00295     Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
00296     DfDx.PutScalar(0.0);
00297     //
00298     // Fill W = DfDx
00299     //
00300     // W = DfDx = [      1.0,  2*x[1] ]
00301     //            [ 2*d*x[0],     -d  ]
00302     //
00303     double values[2];
00304     int indexes[2];
00305     // Row [0]
00306     values[0] = 1.0;           indexes[0] = 0;
00307     values[1] = 2.0*x[1];      indexes[1] = 1;
00308     DfDx.SumIntoGlobalValues( 0, 2, values, indexes );
00309     // Row [1]
00310     values[0] = 2.0*d_*x[0];   indexes[0] = 0;
00311     values[1] = -d_;           indexes[1] = 1;
00312     DfDx.SumIntoGlobalValues( 1, 2, values, indexes );
00313   }
00314   if(DfDp_out) {
00315     Epetra_MultiVector &DfDp = *DfDp_out;
00316     DfDp.PutScalar(0.0);
00317     DfDp[0][0] = -1.0;
00318     DfDp[1][1] = -d_;
00319   }
00320   if(DgDx_trans_out) {
00321     Epetra_Vector &DgDx_trans = *(*DgDx_trans_out)(0);
00322     DgDx_trans[0] = x[0]-xt0_;
00323     DgDx_trans[1] = x[1]-xt1_;
00324   }
00325   if(DgDp_trans_out) {
00326     Epetra_Vector &DgDp_trans = *(*DgDp_trans_out)(0);
00327     DgDp_trans[0] = p[0]-pt0_;
00328     DgDp_trans[1] = p[1]-pt1_;
00329   }
00330 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines