Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DiagonalQuadraticResponseOnlyModelEvaluator_def.hpp
Go to the documentation of this file.
00001 #ifndef OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
00002 #define OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
00003 
00004 
00005 #include "Thyra_DiagonalQuadraticResponseOnlyModelEvaluator_decl.hpp"
00006 #include "Thyra_DiagonalScalarProd.hpp"
00007 #include "Thyra_VectorStdOps.hpp"
00008 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00009 #include "Thyra_DetachedSpmdVectorView.hpp"
00010 #include "Teuchos_DefaultComm.hpp"
00011 #include "Teuchos_CommHelpers.hpp"
00012 #include "Teuchos_Assert.hpp"
00013 
00014 
00015 namespace Thyra {
00016 
00017 
00018 //
00019 // DiagonalQuadraticResponseOnlyModelEvaluator
00020 //
00021 
00022 
00023 // Constructors, Initialization, Misc.
00024 
00025 
00026 template<class Scalar>
00027 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::DiagonalQuadraticResponseOnlyModelEvaluator(
00028   const int localDim,
00029   const RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm
00030   )
00031   :Np_(1), Ng_(1), comm_(comm), localDim_(localDim),
00032    nonlinearTermFactor_(0.0), g_offset_(0.0)
00033 {
00034 
00035   typedef ScalarTraits<Scalar> ST;
00036   using Thyra::createMember;
00037 
00038   TEUCHOS_ASSERT( localDim > 0 );
00039 
00040   // Get the comm
00041   if (is_null(comm_)) {
00042     comm_ = Teuchos::DefaultComm<Thyra::Ordinal>::getComm();
00043   }
00044 
00045   // Locally replicated space for g
00046   g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm_, 1, 1);
00047 
00048   // Distributed space for p
00049   p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm_, localDim, -1);
00050 
00051   // Default solution
00052   const RCP<Thyra::VectorBase<Scalar> > ps = createMember<Scalar>(p_space_);
00053   V_S(ps.ptr(), ST::zero());
00054   ps_ = ps;
00055 
00056   // Default diagonal
00057   const RCP<Thyra::VectorBase<Scalar> > diag = createMember<Scalar>(p_space_);
00058   V_S(diag.ptr(), ST::one());
00059   diag_ = diag;
00060   diag_bar_ = diag;
00061 
00062   // Default inner product
00063   const RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_);
00064   V_S(s_bar.ptr(), ST::one());
00065   s_bar_ = s_bar;
00066 
00067   // Default response offset
00068   g_offset_ = ST::zero();
00069 
00070 }
00071 
00072 
00073 template<class Scalar>
00074 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setSolutionVector(
00075   const RCP<const Thyra::VectorBase<Scalar> > &ps)
00076 {
00077   ps_ = ps.assert_not_null();
00078 }
00079 
00080 
00081 template<class Scalar>
00082 const RCP<const Thyra::VectorBase<Scalar> >
00083 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::getSolutionVector() const
00084 {
00085   return ps_;
00086 }
00087 
00088 
00089 template<class Scalar>
00090 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setDiagonalVector(
00091   const RCP<const Thyra::VectorBase<Scalar> > &diag)
00092 {
00093   diag_ = diag;
00094   diag_bar_ = diag;
00095 }
00096 
00097 
00098 template<class Scalar>
00099 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setDiagonalBarVector(
00100   const RCP<const Thyra::VectorBase<Scalar> > &diag_bar)
00101 {
00102 
00103   typedef Teuchos::ScalarTraits<Scalar> ST;
00104   using Teuchos::rcp_dynamic_cast;
00105   using Thyra::createMember;
00106   using Thyra::ele_wise_divide;
00107   using Thyra::V_S;
00108 
00109   diag_bar_ = diag_bar.assert_not_null();
00110 
00111   // Reset the scalar product for p_space!
00112 
00113   RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_->clone());
00114   // NOTE: We have to clone the vector space in order to avoid creating a
00115   // circular reference between the space and the vector that defines the
00116   // scalar product for the vector space.
00117 
00118   // s_bar[i] = diag[i] / diag_bar[i] 
00119   V_S( s_bar.ptr(), ST::zero() );
00120   ele_wise_divide( ST::one(), *diag_, *diag_bar_, s_bar.ptr() );
00121   s_bar_ = s_bar;
00122   
00123   const RCP<Thyra::ScalarProdVectorSpaceBase<Scalar> > sp_p_space =
00124     rcp_dynamic_cast<Thyra::ScalarProdVectorSpaceBase<Scalar> >(p_space_, true);
00125   sp_p_space->setScalarProd(diagonalScalarProd<Scalar>(s_bar_));
00126 
00127 }
00128 
00129 
00130 template<class Scalar>
00131 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setNonlinearTermFactor(
00132   const Scalar &nonlinearTermFactor)
00133 {
00134   nonlinearTermFactor_ = nonlinearTermFactor;
00135 }
00136 
00137 
00138 template<class Scalar>
00139 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setScalarOffset(
00140   const Scalar &g_offset)
00141 {
00142   g_offset_ = g_offset;
00143 }
00144 
00145 
00146 // Public functions overridden from ModelEvaulator
00147 
00148 
00149 template<class Scalar>
00150 int DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::Np() const
00151 {
00152   return Np_;
00153 }
00154 
00155 
00156 template<class Scalar>
00157 int DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::Ng() const
00158 {
00159   return Ng_;
00160 }
00161 
00162 
00163 template<class Scalar>
00164 RCP<const Thyra::VectorSpaceBase<Scalar> >
00165 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::get_p_space(int l) const
00166 {
00167 #ifdef TEUCHOS_DEBUG
00168   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
00169 #endif
00170   return p_space_;
00171 }
00172 
00173 
00174 template<class Scalar>
00175 RCP<const Thyra::VectorSpaceBase<Scalar> >
00176 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::get_g_space(int j) const
00177 {
00178 #ifdef TEUCHOS_DEBUG
00179   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
00180 #endif
00181   return g_space_;
00182 }
00183 
00184 
00185 template<class Scalar>
00186 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00187 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::createInArgs() const
00188 {
00189   typedef Thyra::ModelEvaluatorBase MEB;
00190   MEB::InArgsSetup<Scalar> inArgs;
00191   inArgs.setModelEvalDescription(this->description());
00192   inArgs.set_Np(Np_);
00193   return inArgs;
00194 }
00195 
00196 
00197 // Private functions overridden from ModelEvaulatorDefaultBase
00198 
00199 
00200 template<class Scalar>
00201 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00202 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::createOutArgsImpl() const
00203 {
00204   typedef Thyra::ModelEvaluatorBase MEB;
00205   typedef MEB::DerivativeSupport DS;
00206   MEB::OutArgsSetup<Scalar> outArgs;
00207   outArgs.setModelEvalDescription(this->description());
00208   outArgs.set_Np_Ng(Np_,Ng_);
00209   outArgs.setSupports(MEB::OUT_ARG_DgDp, 0 ,0, MEB::DERIV_TRANS_MV_BY_ROW);
00210   return outArgs;
00211 }
00212 
00213 
00214 template<class Scalar>
00215 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::evalModelImpl(
00216   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
00217   const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
00218   ) const
00219 {
00220 
00221   using Teuchos::as;
00222   using Teuchos::outArg;
00223   typedef Teuchos::ScalarTraits<Scalar> ST;
00224   using Thyra::get_mv;
00225   using Thyra::ConstDetachedSpmdVectorView;
00226   using Thyra::DetachedSpmdVectorView;
00227   typedef Thyra::Ordinal Ordinal;
00228   typedef Thyra::ModelEvaluatorBase MEB;
00229   typedef MEB::DerivativeMultiVector<Scalar> DMV;
00230 
00231   const ConstDetachedSpmdVectorView<Scalar> p(inArgs.get_p(0));
00232   const ConstDetachedSpmdVectorView<Scalar> ps(ps_);
00233   const ConstDetachedSpmdVectorView<Scalar> diag(diag_);
00234   const ConstDetachedSpmdVectorView<Scalar> s_bar(s_bar_);
00235 
00236   // g(p)
00237   if (!is_null(outArgs.get_g(0))) {
00238     Scalar g_val = ST::zero();
00239     for (Ordinal i = 0; i < p.subDim(); ++i) {
00240       const Scalar p_ps = p[i] - ps[i];
00241       g_val += diag[i] * p_ps*p_ps;
00242       if (nonlinearTermFactor_ != ST::zero()) {
00243         g_val += nonlinearTermFactor_ * p_ps * p_ps * p_ps;
00244       }
00245     }
00246     Scalar global_g_val;
00247     Teuchos::reduceAll<Ordinal, Scalar>(*comm_, Teuchos::REDUCE_SUM, 
00248       g_val, outArg(global_g_val) );
00249     DetachedSpmdVectorView<Scalar>(outArgs.get_g(0))[0] =
00250       as<Scalar>(0.5) * global_g_val + g_offset_;
00251   }
00252 
00253   // DgDp[i]
00254   if (!outArgs.get_DgDp(0,0).isEmpty()) {
00255     const RCP<Thyra::MultiVectorBase<Scalar> > DgDp_trans_mv =
00256       get_mv<Scalar>(outArgs.get_DgDp(0,0), "DgDp^T", MEB::DERIV_TRANS_MV_BY_ROW);
00257     const DetachedSpmdVectorView<Scalar> DgDp_grad(DgDp_trans_mv->col(0));
00258     for (Thyra::Ordinal i = 0; i < p.subDim(); ++i) {
00259       const Scalar p_ps = p[i] - ps[i];
00260       Scalar DgDp_grad_i = diag[i] * p_ps;
00261       if (nonlinearTermFactor_ != ST::zero()) {
00262         DgDp_grad_i += as<Scalar>(1.5) * nonlinearTermFactor_ * p_ps * p_ps;
00263       }
00264       DgDp_grad[i] = DgDp_grad_i / s_bar[i];
00265 
00266     }
00267   }
00268   
00269 }
00270 
00271 
00272 } // namespace Thyra
00273 
00274 
00275 #endif // OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines