Sacado Package Browser (Single Doxygen Collection): Sacado_ELRFad

00001 // $Id: Sacado_ELRFad_Ops.hpp,v 1.4.2.1 2009/03/05 20:35:25 etphipp Exp $ 
00002 // $Source: /space/CVS/Trilinos/packages/sacado/src/Sacado_ELRFad_Ops.hpp,v $ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 //
00031 // The forward-mode AD classes in Sacado are a derivative work of the
00032 // expression template classes in the Fad package by Nicolas Di Cesare.  
00033 // The following banner is included in the original Fad source code:
00034 //
00035 // ************ DO NOT REMOVE THIS BANNER ****************
00036 //
00037 //  Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
00038 //  http://www.ann.jussieu.fr/~dicesare
00039 //
00040 //            CEMRACS 98 : C++ courses, 
00041 //         templates : new C++ techniques 
00042 //            for scientific computing 
00043 // 
00044 //********************************************************
00045 //
00046 //  A short implementation ( not all operators and 
00047 //  functions are overloaded ) of 1st order Automatic
00048 //  Differentiation in forward mode (FAD) using
00049 //  EXPRESSION TEMPLATES.
00050 //
00051 //********************************************************
00052 // @HEADER
00053 
00054 #ifndef SACADO_ELRFAD_OPS_HPP
00055 #define SACADO_ELRFAD_OPS_HPP
00056 
00057 #include "Sacado_ELRFad_Expression.hpp"
00058 #include <cmath>
00059 #include <algorithm>  // for std::min and std::max
00060 #include <ostream>  // for std::ostream
00061 
00062 #define FAD_UNARYOP_MACRO(OPNAME,OP,VALUE,ADJOINT)      \
00063 namespace Sacado {              \
00064   namespace ELRFad {              \
00065                   \
00066     template <typename ExprT>           \
00067     class OP {};              \
00068                   \
00069     template <typename ExprT>           \
00070     class Expr< OP<ExprT> > {           \
00071     public:               \
00072                   \
00073       typedef typename ExprT::value_type value_type;      \
00074                   \
00075       typedef typename ExprT::base_expr_type base_expr_type;    \
00076                   \
00077       static const int num_args = ExprT::num_args;      \
00078                   \
00079       Expr(const ExprT& expr_) : expr(expr_)  {}      \
00080                   \
00081       int size() const { return expr.size(); }        \
00082                   \
00083       template <int Arg>            \
00084       bool isActive() const { return expr.template isActive<Arg>(); } \
00085                   \
00086       value_type val() const {            \
00087   return VALUE;             \
00088       }                 \
00089                   \
00090       void computePartials(const value_type& bar,     \
00091          value_type partials[]) const {   \
00092   expr.computePartials(ADJOINT, partials);      \
00093       }                 \
00094                   \
00095       void getTangents(int i, value_type dots[]) const {    \
00096   expr.getTangents(i, dots); }          \
00097                   \
00098       template <int Arg>            \
00099       value_type getTangent(int i) const {        \
00100   return expr.template getTangent<Arg>(i);      \
00101       }                 \
00102                   \
00103     protected:                \
00104                   \
00105       const ExprT& expr;            \
00106     };                  \
00107                   \
00108     template <typename T>           \
00109     inline Expr< OP< Expr<T> > >          \
00110     OPNAME (const Expr<T>& expr)          \
00111     {                 \
00112       typedef OP< Expr<T> > expr_t;         \
00113                         \
00114       return Expr<expr_t>(expr);          \
00115     }                 \
00116   }                 \
00117 }
00118 
00119 FAD_UNARYOP_MACRO(operator+,
00120       UnaryPlusOp, 
00121       expr.val(),
00122       bar)
00123 FAD_UNARYOP_MACRO(operator-,
00124       UnaryMinusOp, 
00125       -expr.val(),
00126       -bar)
00127 FAD_UNARYOP_MACRO(exp,
00128       ExpOp, 
00129       std::exp(expr.val()),
00130       bar*std::exp(expr.val()))
00131 FAD_UNARYOP_MACRO(log,
00132       LogOp, 
00133       std::log(expr.val()),
00134       bar/expr.val())
00135 FAD_UNARYOP_MACRO(log10,
00136       Log10Op, 
00137       std::log10(expr.val()),
00138       bar/( std::log(value_type(10.))*expr.val() ))
00139 FAD_UNARYOP_MACRO(sqrt,
00140       SqrtOp, 
00141       std::sqrt(expr.val()),
00142       value_type(0.5)*bar/std::sqrt(expr.val()))
00143 FAD_UNARYOP_MACRO(cos,
00144       CosOp, 
00145       std::cos(expr.val()),
00146       -bar*std::sin(expr.val()))
00147 FAD_UNARYOP_MACRO(sin,
00148       SinOp, 
00149       std::sin(expr.val()),
00150       bar*std::cos(expr.val()))
00151 FAD_UNARYOP_MACRO(tan,
00152       TanOp, 
00153       std::tan(expr.val()),
00154       bar*(value_type(1.)+ std::tan(expr.val())*std::tan(expr.val())))
00155 FAD_UNARYOP_MACRO(acos,
00156       ACosOp, 
00157       std::acos(expr.val()),
00158       -bar/std::sqrt(value_type(1.)-expr.val()*expr.val()))
00159 FAD_UNARYOP_MACRO(asin,
00160       ASinOp, 
00161       std::asin(expr.val()),
00162       bar/std::sqrt(value_type(1.)-expr.val()*expr.val()))
00163 FAD_UNARYOP_MACRO(atan,
00164       ATanOp, 
00165       std::atan(expr.val()),
00166       bar/(value_type(1.)+expr.val()*expr.val()))
00167 FAD_UNARYOP_MACRO(cosh,
00168       CoshOp, 
00169       std::cosh(expr.val()),
00170       bar*std::sinh(expr.val()))
00171 FAD_UNARYOP_MACRO(sinh,
00172       SinhOp, 
00173       std::sinh(expr.val()),
00174       bar*std::cosh(expr.val()))
00175 FAD_UNARYOP_MACRO(tanh,
00176       TanhOp, 
00177       std::tanh(expr.val()),
00178       bar/(std::cosh(expr.val())*std::cosh(expr.val())))
00179 FAD_UNARYOP_MACRO(acosh,
00180       ACoshOp, 
00181       acosh(expr.val()),
00182       bar/std::sqrt((expr.val()-value_type(1.)) * 
00183         (expr.val()+value_type(1.))))
00184 FAD_UNARYOP_MACRO(asinh,
00185       ASinhOp, 
00186       asinh(expr.val()),
00187       bar/std::sqrt(value_type(1.)+expr.val()*expr.val()))
00188 FAD_UNARYOP_MACRO(atanh,
00189       ATanhOp, 
00190       atanh(expr.val()),
00191       bar/(value_type(1.)-expr.val()*expr.val()))
00192 FAD_UNARYOP_MACRO(abs,
00193       AbsOp, 
00194       std::abs(expr.val()),
00195       (expr.val() >= value_type(0.)) ? bar : value_type(-bar))
00196 FAD_UNARYOP_MACRO(fabs,
00197       FAbsOp, 
00198       std::fabs(expr.val()),
00199       (expr.val() >= value_type(0.)) ? bar : value_type(-bar))
00200 
00201 #undef FAD_UNARYOP_MACRO
00202 
00203 #define FAD_BINARYOP_MACRO(OPNAME,OP,VALUE,LADJOINT,RADJOINT)   \
00204 namespace Sacado {              \
00205   namespace ELRFad {              \
00206                   \
00207     template <typename ExprT1, typename ExprT2>       \
00208     class OP {};              \
00209                   \
00210     template <typename ExprT1, typename ExprT2>       \
00211     class Expr< OP<ExprT1,ExprT2> > {         \
00212                   \
00213     public:               \
00214                   \
00215       typedef typename ExprT1::value_type value_type_1;     \
00216       typedef typename ExprT2::value_type value_type_2;     \
00217       typedef typename Sacado::Promote<value_type_1,      \
00218                value_type_2>::type value_type;  \
00219                   \
00220       typedef typename ExprT1::base_expr_type base_expr_type_1;   \
00221       typedef typename ExprT2::base_expr_type base_expr_type_2;   \
00222       typedef typename ExprPromote<base_expr_type_1,      \
00223            base_expr_type_2>::type base_expr_type; \
00224                   \
00225       static const int num_args1 = ExprT1::num_args;      \
00226       static const int num_args2 = ExprT2::num_args;      \
00227       static const int num_args = num_args1 + num_args2;    \
00228                   \
00229       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :    \
00230   expr1(expr1_), expr2(expr2_) {}         \
00231                   \
00232       int size() const {            \
00233   int sz1 = expr1.size(), sz2 = expr2.size();     \
00234   return sz1 > sz2 ? sz1 : sz2;         \
00235       }                 \
00236                   \
00237       template <int Arg> bool isActive() const {      \
00238   if (Arg < num_args1)            \
00239     return expr1.template isActive<Arg>();      \
00240   else                \
00241     return expr2.template isActive<Arg-num_args1>();    \
00242       }                 \
00243                   \
00244       value_type val() const {            \
00245   return VALUE;             \
00246       }                 \
00247                   \
00248       void computePartials(const value_type& bar,     \
00249          value_type partials[]) const {   \
00250   if (num_args1 > 0)            \
00251     expr1.computePartials(LADJOINT, partials);      \
00252   if (num_args2 > 0)            \
00253     expr2.computePartials(RADJOINT, partials+num_args1);    \
00254       }                 \
00255                   \
00256       void getTangents(int i, value_type dots[]) const {    \
00257   expr1.getTangents(i, dots);         \
00258   expr2.getTangents(i, dots+num_args1);       \
00259       }                 \
00260                   \
00261       template <int Arg> value_type getTangent(int i) const {   \
00262   if (Arg < num_args1)            \
00263     return expr1.template getTangent<Arg>(i);     \
00264   else                \
00265     return expr2.template getTangent<Arg-num_args1>(i);   \
00266       }                 \
00267                   \
00268     protected:                \
00269                   \
00270       typename ExprConstRef<ExprT1>::type expr1;      \
00271       typename ExprConstRef<ExprT2>::type expr2;      \
00272                   \
00273     };                  \
00274                   \
00275     template <typename T1, typename T2>         \
00276     inline Expr< OP< Expr<T1>, Expr<T2> > >       \
00277     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
00278     {                 \
00279       typedef OP< Expr<T1>, Expr<T2> > expr_t;        \
00280                       \
00281       return Expr<expr_t>(expr1, expr2);        \
00282     }                 \
00283                   \
00284     template <typename T>           \
00285     inline Expr< OP< Expr<T>, Expr<T> > >       \
00286     OPNAME (const Expr<T>& expr1, const Expr<T>& expr2)     \
00287     {                 \
00288       typedef OP< Expr<T>, Expr<T> > expr_t;        \
00289                       \
00290       return Expr<expr_t>(expr1, expr2);        \
00291     }                 \
00292                   \
00293     template <typename T>           \
00294     inline Expr< OP< ConstExpr<typename Expr<T>::value_type>,   \
00295          Expr<T> > >          \
00296     OPNAME (const typename Expr<T>::value_type& c,      \
00297       const Expr<T>& expr)          \
00298     {                 \
00299       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00300       typedef OP< ConstT, Expr<T> > expr_t;       \
00301                   \
00302       return Expr<expr_t>(ConstT(c), expr);       \
00303     }                 \
00304                   \
00305     template <typename T>           \
00306     inline Expr< OP< Expr<T>,           \
00307          ConstExpr<typename Expr<T>::value_type> > >  \
00308     OPNAME (const Expr<T>& expr,          \
00309       const typename Expr<T>::value_type& c)      \
00310     {                 \
00311       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00312       typedef OP< Expr<T>, ConstT > expr_t;       \
00313                   \
00314       return Expr<expr_t>(expr, ConstT(c));       \
00315     }                 \
00316   }                 \
00317 }
00318 
00319 
00320 FAD_BINARYOP_MACRO(operator+,
00321        AdditionOp, 
00322        expr1.val() + expr2.val(),
00323        bar,
00324        bar)
00325 FAD_BINARYOP_MACRO(operator-,
00326        SubtractionOp, 
00327        expr1.val() - expr2.val(),
00328        bar,
00329        -bar)
00330 FAD_BINARYOP_MACRO(operator*,
00331        MultiplicationOp, 
00332        expr1.val() * expr2.val(),
00333        bar*expr2.val(),
00334        bar*expr1.val())
00335 FAD_BINARYOP_MACRO(operator/,
00336        DivisionOp, 
00337        expr1.val() / expr2.val(),
00338        bar/expr2.val(),
00339        -bar*expr1.val()/(expr2.val()*expr2.val()))
00340 FAD_BINARYOP_MACRO(atan2,
00341        Atan2Op,
00342        std::atan2(expr1.val(), expr2.val()),
00343        bar*expr2.val()/
00344        (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00345        bar*expr1.val()/
00346        (expr1.val()*expr1.val() + expr2.val()*expr2.val()))
00347 FAD_BINARYOP_MACRO(pow,
00348        PowerOp,
00349        std::pow(expr1.val(), expr2.val()),
00350        bar*std::pow(expr1.val(),expr2.val())*expr2.val()/
00351        expr1.val(),
00352        bar*std::pow(expr1.val(),expr2.val())*std::log(expr1.val()))
00353 FAD_BINARYOP_MACRO(max,
00354                    MaxOp,
00355                    std::max(expr1.val(), expr2.val()),
00356                    expr1.val() >= expr2.val() ? bar : value_type(0.),
00357                    expr2.val() > expr1.val() ? bar : value_type(0.))
00358 FAD_BINARYOP_MACRO(min,
00359                    MinOp,
00360                    std::min(expr1.val(), expr2.val()),
00361                    expr1.val() <= expr2.val() ? bar : value_type(0.),
00362                    expr2.val() < expr1.val() ? bar : value_type(0.))
00363 
00364 #undef FAD_BINARYOP_MACRO
00365 
00366 //-------------------------- Relational Operators -----------------------
00367 
00368 #define FAD_RELOP_MACRO(OP)           \
00369 namespace Sacado {              \
00370   namespace ELRFad {              \
00371     template <typename ExprT1, typename ExprT2>       \
00372     inline bool               \
00373     operator OP (const Expr<ExprT1>& expr1,       \
00374      const Expr<ExprT2>& expr2)       \
00375     {                 \
00376       return expr1.val() OP expr2.val();        \
00377     }                 \
00378                   \
00379     template <typename ExprT2>            \
00380     inline bool               \
00381     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00382      const Expr<ExprT2>& expr2)       \
00383     {                 \
00384       return a OP expr2.val();            \
00385     }                 \
00386                   \
00387     template <typename ExprT1>            \
00388     inline bool               \
00389     operator OP (const Expr<ExprT1>& expr1,       \
00390      const typename Expr<ExprT1>::value_type& b)    \
00391     {                 \
00392       return expr1.val() OP b;            \
00393     }                 \
00394   }                 \
00395 }
00396 
00397 FAD_RELOP_MACRO(==)
00398 FAD_RELOP_MACRO(!=)
00399 FAD_RELOP_MACRO(<)
00400 FAD_RELOP_MACRO(>)
00401 FAD_RELOP_MACRO(<=)
00402 FAD_RELOP_MACRO(>=)
00403 FAD_RELOP_MACRO(<<=)
00404 FAD_RELOP_MACRO(>>=)
00405 FAD_RELOP_MACRO(&)
00406 FAD_RELOP_MACRO(|)
00407 
00408 #undef FAD_RELOP_MACRO
00409 
00410 namespace Sacado {
00411 
00412   namespace ELRFad {
00413 
00414     template <typename ExprT>
00415     inline bool operator ! (const Expr<ExprT>& expr) 
00416     {
00417       return ! expr.val();
00418     }
00419 
00420   } // namespace ELRFad
00421 
00422 } // namespace Sacado
00423 
00424 //-------------------------- I/O Operators -----------------------
00425 
00426 namespace Sacado {
00427 
00428   namespace ELRFad {
00429 
00430     template <typename ExprT>
00431     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
00432       typedef typename Expr<ExprT>::base_expr_type base_expr_type;
00433       return os << base_expr_type(x);
00434     }
00435 
00436   } // namespace Fad
00437 
00438 } // namespace Sacado
00439 
00440 
00441 #endif // SACADO_FAD_OPS_HPP