EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraModelEval4DOpt.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 "EpetraModelEval4DOpt.hpp"
00045 #include "Teuchos_ScalarTraits.hpp"
00046 #include "Epetra_SerialComm.h"
00047 #include "Epetra_CrsMatrix.h"
00048 
00049 #ifdef HAVE_MPI
00050 #  include "Epetra_MpiComm.h"
00051 #else
00052 #  include "Epetra_SerialComm.h"
00053 #endif
00054 
00055 namespace {
00056 
00057 inline double sqr( const double& s ) { return s*s; }
00058 
00059 } // namespace
00060 
00061 EpetraModelEval4DOpt::EpetraModelEval4DOpt(
00062   const double         xt0
00063   ,const double        xt1
00064   ,const double        pt0
00065   ,const double        pt1
00066   ,const double        d
00067   ,const double        x00
00068   ,const double        x01
00069   ,const double        p00
00070   ,const double        p01
00071   )
00072   :xt0_(xt0),xt1_(xt1),pt0_(pt0),pt1_(pt1),d_(d),
00073    isInitialized_(false),supportDerivs_(true)
00074 {
00075   using Teuchos::rcp;
00076 
00077   epetra_comm_ = rcp(new Epetra_SerialComm());
00078 
00079   const int nx = 2, np = 2, ng = 1;
00080 
00081   map_x_ = rcp(new Epetra_Map(nx,0,*epetra_comm_));
00082   map_p_ = rcp(new Epetra_Map(np,0,*epetra_comm_));
00083   map_g_ = rcp(new Epetra_Map(ng,0,*epetra_comm_));
00084 
00085   //const double inf = infiniteBound();
00086   const double inf = 1e+50;
00087 
00088   xL_ = rcp(new Epetra_Vector(*map_x_));  xL_->PutScalar(-inf);
00089   xU_ = rcp(new Epetra_Vector(*map_x_));  xU_->PutScalar(+inf);
00090   x0_ = rcp(new Epetra_Vector(*map_x_));  (*x0_)[0] = x00; (*x0_)[1] = x01;
00091   pL_ = rcp(new Epetra_Vector(*map_p_));  pL_->PutScalar(-inf);
00092   pU_ = rcp(new Epetra_Vector(*map_p_));  pU_->PutScalar(+inf);
00093   p0_ = rcp(new Epetra_Vector(*map_p_));  (*p0_)[0] = p00; (*p0_)[1] = p01;
00094   gL_ = rcp(new Epetra_Vector(*map_g_));  gL_->PutScalar(-inf);
00095   gU_ = rcp(new Epetra_Vector(*map_g_));  gU_->PutScalar(+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 
00111 void EpetraModelEval4DOpt::setSupportDerivs( bool supportDerivs )
00112 {
00113   supportDerivs_ = supportDerivs;
00114 }
00115 
00116 
00117 void EpetraModelEval4DOpt::set_p_bounds(
00118   double pL0, double pL1, double pU0, double pU1
00119   )
00120 {
00121   (*pL_)[0] = pL0;
00122   (*pL_)[1] = pL1;
00123   (*pU_)[0] = pU0;
00124   (*pU_)[1] = pU1;
00125 }
00126 
00127 void EpetraModelEval4DOpt::set_x_bounds(
00128   double xL0, double xL1, double xU0, double xU1
00129   )
00130 {
00131   (*xL_)[0] = xL0;
00132   (*xL_)[1] = xL1;
00133   (*xU_)[0] = xU0;
00134   (*xU_)[1] = xU1;
00135 }
00136 
00137 // Overridden from EpetraExt::ModelEvaluator
00138 
00139 Teuchos::RefCountPtr<const Epetra_Map>
00140 EpetraModelEval4DOpt::get_x_map() const
00141 {
00142   return map_x_;
00143 }
00144 
00145 Teuchos::RefCountPtr<const Epetra_Map>
00146 EpetraModelEval4DOpt::get_f_map() const
00147 {
00148   return map_x_;
00149 }
00150 
00151 Teuchos::RefCountPtr<const Epetra_Map>
00152 EpetraModelEval4DOpt::get_p_map(int l) const
00153 {
00154   TEUCHOS_TEST_FOR_EXCEPT(l!=0);
00155   return map_p_;
00156 }
00157 
00158 Teuchos::RefCountPtr<const Epetra_Map>
00159 EpetraModelEval4DOpt::get_g_map(int j) const
00160 {
00161   TEUCHOS_TEST_FOR_EXCEPT(j!=0);
00162   return map_g_;
00163 }
00164 
00165 Teuchos::RefCountPtr<const Epetra_Vector>
00166 EpetraModelEval4DOpt::get_x_init() const
00167 {
00168   return x0_;
00169 }
00170 
00171 Teuchos::RefCountPtr<const Epetra_Vector>
00172 EpetraModelEval4DOpt::get_p_init(int l) const
00173 {
00174   TEUCHOS_TEST_FOR_EXCEPT(l!=0);
00175   return p0_;
00176 }
00177 
00178 Teuchos::RefCountPtr<const Epetra_Vector>
00179 EpetraModelEval4DOpt::get_x_lower_bounds() const
00180 {
00181   return xL_;
00182 }
00183 
00184 Teuchos::RefCountPtr<const Epetra_Vector>
00185 EpetraModelEval4DOpt::get_x_upper_bounds() const
00186 {
00187   return xU_;
00188 }
00189 
00190 Teuchos::RefCountPtr<const Epetra_Vector>
00191 EpetraModelEval4DOpt::get_p_lower_bounds(int l) const
00192 {
00193   TEUCHOS_TEST_FOR_EXCEPT(l!=0);
00194   return pL_;
00195 }
00196 
00197 Teuchos::RefCountPtr<const Epetra_Vector>
00198 EpetraModelEval4DOpt::get_p_upper_bounds(int l) const
00199 {
00200   TEUCHOS_TEST_FOR_EXCEPT(l!=0);
00201   return pU_;
00202 }
00203 
00204 Teuchos::RefCountPtr<Epetra_Operator>
00205 EpetraModelEval4DOpt::create_W() const
00206 {
00207   return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00208 }
00209 
00210 EpetraExt::ModelEvaluator::InArgs
00211 EpetraModelEval4DOpt::createInArgs() const
00212 {
00213   InArgsSetup inArgs;
00214   inArgs.setModelEvalDescription(this->description());
00215   inArgs.set_Np(1);
00216   inArgs.setSupports(IN_ARG_x,true);
00217   return inArgs;
00218 }
00219 
00220 EpetraExt::ModelEvaluator::OutArgs
00221 EpetraModelEval4DOpt::createOutArgs() const
00222 {
00223   OutArgsSetup outArgs;
00224   outArgs.setModelEvalDescription(this->description());
00225   outArgs.set_Np_Ng(1,1);
00226   outArgs.setSupports(OUT_ARG_f,true);
00227   outArgs.setSupports(OUT_ARG_W,true);
00228   outArgs.set_W_properties(
00229     DerivativeProperties(
00230       DERIV_LINEARITY_NONCONST
00231       ,DERIV_RANK_FULL
00232       ,true // supportsAdjoint
00233       )
00234     );
00235   if (supportDerivs_) {
00236     outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
00237     outArgs.set_DfDp_properties(
00238       0,DerivativeProperties(
00239         DERIV_LINEARITY_CONST
00240         ,DERIV_RANK_DEFICIENT
00241         ,true // supportsAdjoint
00242         )
00243       );
00244     outArgs.setSupports(OUT_ARG_DgDx,0,DERIV_TRANS_MV_BY_ROW);
00245     outArgs.set_DgDx_properties(
00246       0,DerivativeProperties(
00247         DERIV_LINEARITY_NONCONST
00248         ,DERIV_RANK_DEFICIENT
00249         ,true // supportsAdjoint
00250         )
00251       );
00252     outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_TRANS_MV_BY_ROW);
00253     outArgs.set_DgDp_properties(
00254       0,0,DerivativeProperties(
00255         DERIV_LINEARITY_NONCONST
00256         ,DERIV_RANK_DEFICIENT
00257         ,true // supportsAdjoint
00258         )
00259       );
00260   }
00261   return outArgs;
00262 }
00263 
00264 
00265 void EpetraModelEval4DOpt::evalModel(
00266   const InArgs& inArgs, const OutArgs& outArgs
00267   ) const
00268 {
00269   using Teuchos::dyn_cast;
00270   using Teuchos::rcp_dynamic_cast;
00271   //
00272   // Get the input arguments
00273   //
00274   Teuchos::RefCountPtr<const Epetra_Vector> p_in = inArgs.get_p(0);
00275   const Epetra_Vector &p = (p_in.get() ? *p_in : *p0_);
00276   const Epetra_Vector &x = *inArgs.get_x();
00277   //
00278   // Get the output arguments
00279   //
00280   Epetra_Vector       *f_out = outArgs.get_f().get();
00281   Epetra_Vector       *g_out = outArgs.get_g(0).get();
00282   Epetra_Operator     *W_out = outArgs.get_W().get();
00283   Epetra_MultiVector  *DfDp_out = supportDerivs_ ? get_DfDp_mv(0,outArgs).get() : 0;
00284   Epetra_MultiVector  *DgDx_trans_out = supportDerivs_ ? get_DgDx_mv(0,outArgs,DERIV_TRANS_MV_BY_ROW).get() : 0;
00285   Epetra_MultiVector  *DgDp_trans_out = supportDerivs_ ? get_DgDp_mv(0,0,outArgs,DERIV_TRANS_MV_BY_ROW).get() : 0;
00286   //
00287   // Compute the functions
00288   //
00289   if(f_out) {
00290     Epetra_Vector &f = *f_out;
00291     // zero-based indexing for Epetra!
00292     f[0] =        x[0]      + x[1]*x[1] - p[0];
00293     f[1] = d_ * ( x[0]*x[0] - x[1]      - p[1] );
00294   }
00295   if(g_out) {
00296     Epetra_Vector &g = *g_out;
00297     g[0] = 0.5 * ( sqr(x[0]-xt0_) + sqr(x[1]-xt1_) + sqr(p[0]-pt0_) + sqr(p[1]-pt1_) );
00298   }
00299   if(W_out) {
00300     Epetra_CrsMatrix &DfDx = dyn_cast<Epetra_CrsMatrix>(*W_out);
00301     DfDx.PutScalar(0.0);
00302     //
00303     // Fill W = DfDx
00304     //
00305     // W = DfDx = [      1.0,  2*x[1] ]
00306     //            [ 2*d*x[0],     -d  ]
00307     //
00308     double values[2];
00309     int indexes[2];
00310     // Row [0]
00311     values[0] = 1.0;           indexes[0] = 0;
00312     values[1] = 2.0*x[1];      indexes[1] = 1;
00313     DfDx.SumIntoGlobalValues( 0, 2, values, indexes );
00314     // Row [1]
00315     values[0] = 2.0*d_*x[0];   indexes[0] = 0;
00316     values[1] = -d_;           indexes[1] = 1;
00317     DfDx.SumIntoGlobalValues( 1, 2, values, indexes );
00318   }
00319   if(DfDp_out) {
00320     Epetra_MultiVector &DfDp = *DfDp_out;
00321     DfDp.PutScalar(0.0);
00322     DfDp[0][0] = -1.0;
00323     DfDp[1][1] = -d_;
00324   }
00325   if(DgDx_trans_out) {
00326     Epetra_Vector &DgDx_trans = *(*DgDx_trans_out)(0);
00327     DgDx_trans[0] = x[0]-xt0_;
00328     DgDx_trans[1] = x[1]-xt1_;
00329   }
00330   if(DgDp_trans_out) {
00331     Epetra_Vector &DgDp_trans = *(*DgDp_trans_out)(0);
00332     DgDp_trans[0] = p[0]-pt0_;
00333     DgDp_trans[1] = p[1]-pt1_;
00334   }
00335 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines