Thyra Version of the Day
Thyra_DiagonalQuadraticResponseOnlyModelEvaluator_def.hpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
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 Roscoe A. Bartlett (bartlettra@ornl.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 #ifndef OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
00045 #define OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
00046 
00047 
00048 #include "Thyra_DiagonalQuadraticResponseOnlyModelEvaluator_decl.hpp"
00049 #include "Thyra_DiagonalScalarProd.hpp"
00050 #include "Thyra_VectorStdOps.hpp"
00051 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00052 #include "Thyra_DetachedSpmdVectorView.hpp"
00053 #include "Teuchos_DefaultComm.hpp"
00054 #include "Teuchos_CommHelpers.hpp"
00055 #include "Teuchos_Assert.hpp"
00056 
00057 
00058 namespace Thyra {
00059 
00060 
00061 //
00062 // DiagonalQuadraticResponseOnlyModelEvaluator
00063 //
00064 
00065 
00066 // Constructors, Initialization, Misc.
00067 
00068 
00069 template<class Scalar>
00070 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::DiagonalQuadraticResponseOnlyModelEvaluator(
00071   const int localDim,
00072   const RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm
00073   )
00074   :Np_(1), Ng_(1), comm_(comm), localDim_(localDim),
00075    nonlinearTermFactor_(0.0), g_offset_(0.0)
00076 {
00077 
00078   typedef ScalarTraits<Scalar> ST;
00079   using Thyra::createMember;
00080 
00081   TEUCHOS_ASSERT( localDim > 0 );
00082 
00083   // Get the comm
00084   if (is_null(comm_)) {
00085     comm_ = Teuchos::DefaultComm<Thyra::Ordinal>::getComm();
00086   }
00087 
00088   // Locally replicated space for g
00089   g_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm_, 1, 1);
00090 
00091   // Distributed space for p
00092   p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm_, localDim, -1);
00093 
00094   // Default solution
00095   const RCP<Thyra::VectorBase<Scalar> > ps = createMember<Scalar>(p_space_);
00096   V_S(ps.ptr(), ST::zero());
00097   ps_ = ps;
00098 
00099   // Default diagonal
00100   const RCP<Thyra::VectorBase<Scalar> > diag = createMember<Scalar>(p_space_);
00101   V_S(diag.ptr(), ST::one());
00102   diag_ = diag;
00103   diag_bar_ = diag;
00104 
00105   // Default inner product
00106   const RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_);
00107   V_S(s_bar.ptr(), ST::one());
00108   s_bar_ = s_bar;
00109 
00110   // Default response offset
00111   g_offset_ = ST::zero();
00112 
00113 }
00114 
00115 
00116 template<class Scalar>
00117 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setSolutionVector(
00118   const RCP<const Thyra::VectorBase<Scalar> > &ps)
00119 {
00120   ps_ = ps.assert_not_null();
00121 }
00122 
00123 
00124 template<class Scalar>
00125 const RCP<const Thyra::VectorBase<Scalar> >
00126 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::getSolutionVector() const
00127 {
00128   return ps_;
00129 }
00130 
00131 
00132 template<class Scalar>
00133 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setDiagonalVector(
00134   const RCP<const Thyra::VectorBase<Scalar> > &diag)
00135 {
00136   diag_ = diag;
00137   diag_bar_ = diag;
00138 }
00139 
00140 
00141 template<class Scalar>
00142 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setDiagonalBarVector(
00143   const RCP<const Thyra::VectorBase<Scalar> > &diag_bar)
00144 {
00145 
00146   typedef Teuchos::ScalarTraits<Scalar> ST;
00147   using Teuchos::rcp_dynamic_cast;
00148   using Thyra::createMember;
00149   using Thyra::ele_wise_divide;
00150   using Thyra::V_S;
00151 
00152   diag_bar_ = diag_bar.assert_not_null();
00153 
00154   // Reset the scalar product for p_space!
00155 
00156   RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_->clone());
00157   // NOTE: We have to clone the vector space in order to avoid creating a
00158   // circular reference between the space and the vector that defines the
00159   // scalar product for the vector space.
00160 
00161   // s_bar[i] = diag[i] / diag_bar[i] 
00162   V_S( s_bar.ptr(), ST::zero() );
00163   ele_wise_divide( ST::one(), *diag_, *diag_bar_, s_bar.ptr() );
00164   s_bar_ = s_bar;
00165   
00166   const RCP<Thyra::ScalarProdVectorSpaceBase<Scalar> > sp_p_space =
00167     rcp_dynamic_cast<Thyra::ScalarProdVectorSpaceBase<Scalar> >(p_space_, true);
00168   sp_p_space->setScalarProd(diagonalScalarProd<Scalar>(s_bar_));
00169 
00170 }
00171 
00172 
00173 template<class Scalar>
00174 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setNonlinearTermFactor(
00175   const Scalar &nonlinearTermFactor)
00176 {
00177   nonlinearTermFactor_ = nonlinearTermFactor;
00178 }
00179 
00180 
00181 template<class Scalar>
00182 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::setScalarOffset(
00183   const Scalar &g_offset)
00184 {
00185   g_offset_ = g_offset;
00186 }
00187 
00188 
00189 // Public functions overridden from ModelEvaulator
00190 
00191 
00192 template<class Scalar>
00193 int DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::Np() const
00194 {
00195   return Np_;
00196 }
00197 
00198 
00199 template<class Scalar>
00200 int DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::Ng() const
00201 {
00202   return Ng_;
00203 }
00204 
00205 
00206 template<class Scalar>
00207 RCP<const Thyra::VectorSpaceBase<Scalar> >
00208 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::get_p_space(int l) const
00209 {
00210 #ifdef TEUCHOS_DEBUG
00211   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
00212 #endif
00213   return p_space_;
00214 }
00215 
00216 
00217 template<class Scalar>
00218 RCP<const Thyra::VectorSpaceBase<Scalar> >
00219 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::get_g_space(int j) const
00220 {
00221 #ifdef TEUCHOS_DEBUG
00222   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
00223 #endif
00224   return g_space_;
00225 }
00226 
00227 
00228 template<class Scalar>
00229 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00230 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::createInArgs() const
00231 {
00232   typedef Thyra::ModelEvaluatorBase MEB;
00233   MEB::InArgsSetup<Scalar> inArgs;
00234   inArgs.setModelEvalDescription(this->description());
00235   inArgs.set_Np(Np_);
00236   return inArgs;
00237 }
00238 
00239 
00240 // Private functions overridden from ModelEvaulatorDefaultBase
00241 
00242 
00243 template<class Scalar>
00244 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00245 DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::createOutArgsImpl() const
00246 {
00247   typedef Thyra::ModelEvaluatorBase MEB;
00248   typedef MEB::DerivativeSupport DS;
00249   MEB::OutArgsSetup<Scalar> outArgs;
00250   outArgs.setModelEvalDescription(this->description());
00251   outArgs.set_Np_Ng(Np_,Ng_);
00252   outArgs.setSupports(MEB::OUT_ARG_DgDp, 0 ,0, MEB::DERIV_TRANS_MV_BY_ROW);
00253   return outArgs;
00254 }
00255 
00256 
00257 template<class Scalar>
00258 void DiagonalQuadraticResponseOnlyModelEvaluator<Scalar>::evalModelImpl(
00259   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
00260   const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
00261   ) const
00262 {
00263 
00264   using Teuchos::as;
00265   using Teuchos::outArg;
00266   typedef Teuchos::ScalarTraits<Scalar> ST;
00267   using Thyra::get_mv;
00268   using Thyra::ConstDetachedSpmdVectorView;
00269   using Thyra::DetachedSpmdVectorView;
00270   typedef Thyra::ModelEvaluatorBase MEB;
00271   typedef MEB::DerivativeMultiVector<Scalar> DMV;
00272 
00273   const ConstDetachedSpmdVectorView<Scalar> p(inArgs.get_p(0));
00274   const ConstDetachedSpmdVectorView<Scalar> ps(ps_);
00275   const ConstDetachedSpmdVectorView<Scalar> diag(diag_);
00276   const ConstDetachedSpmdVectorView<Scalar> s_bar(s_bar_);
00277 
00278   // g(p)
00279   if (!is_null(outArgs.get_g(0))) {
00280     Scalar g_val = ST::zero();
00281     for (Ordinal i = 0; i < p.subDim(); ++i) {
00282       const Scalar p_ps = p[i] - ps[i];
00283       g_val += diag[i] * p_ps*p_ps;
00284       if (nonlinearTermFactor_ != ST::zero()) {
00285         g_val += nonlinearTermFactor_ * p_ps * p_ps * p_ps;
00286       }
00287     }
00288     Scalar global_g_val;
00289     Teuchos::reduceAll<Ordinal, Scalar>(*comm_, Teuchos::REDUCE_SUM, 
00290       g_val, outArg(global_g_val) );
00291     DetachedSpmdVectorView<Scalar>(outArgs.get_g(0))[0] =
00292       as<Scalar>(0.5) * global_g_val + g_offset_;
00293   }
00294 
00295   // DgDp[i]
00296   if (!outArgs.get_DgDp(0,0).isEmpty()) {
00297     const RCP<Thyra::MultiVectorBase<Scalar> > DgDp_trans_mv =
00298       get_mv<Scalar>(outArgs.get_DgDp(0,0), "DgDp^T", MEB::DERIV_TRANS_MV_BY_ROW);
00299     const DetachedSpmdVectorView<Scalar> DgDp_grad(DgDp_trans_mv->col(0));
00300     for (Thyra::Ordinal i = 0; i < p.subDim(); ++i) {
00301       const Scalar p_ps = p[i] - ps[i];
00302       Scalar DgDp_grad_i = diag[i] * p_ps;
00303       if (nonlinearTermFactor_ != ST::zero()) {
00304         DgDp_grad_i += as<Scalar>(1.5) * nonlinearTermFactor_ * p_ps * p_ps;
00305       }
00306       DgDp_grad[i] = DgDp_grad_i / s_bar[i];
00307 
00308     }
00309   }
00310   
00311 }
00312 
00313 
00314 } // namespace Thyra
00315 
00316 
00317 #endif // OPTIPACK_DIAGONAL_QUADRATIC_RESPONSE_ONLY_MODEL_EVALUATOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines