Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_ELRFad_Ops.hpp
Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
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 "Sacado_cmath.hpp"
00059 #include <ostream>  // for std::ostream
00060 
00061 #define FAD_UNARYOP_MACRO(OPNAME,OP,VALUE,ADJOINT)      \
00062 namespace Sacado {              \
00063   namespace ELRFad {              \
00064                   \
00065     template <typename ExprT>           \
00066     class OP {};              \
00067                   \
00068     template <typename ExprT>           \
00069     class Expr< OP<ExprT> > {           \
00070     public:               \
00071                   \
00072       typedef typename ExprT::value_type value_type;      \
00073                   \
00074       typedef typename ExprT::base_expr_type base_expr_type;    \
00075                   \
00076       static const int num_args = ExprT::num_args;      \
00077                   \
00078       Expr(const ExprT& expr_) : expr(expr_)  {}      \
00079                   \
00080       int size() const { return expr.size(); }        \
00081                   \
00082       template <int Arg>            \
00083       bool isActive() const { return expr.template isActive<Arg>(); } \
00084                   \
00085       bool updateValue() const { return expr.updateValue(); }   \
00086                   \
00087       value_type val() const {            \
00088   return VALUE;             \
00089       }                 \
00090                   \
00091       void computePartials(const value_type& bar,     \
00092          value_type partials[]) const {   \
00093   expr.computePartials(ADJOINT, partials);      \
00094       }                 \
00095                   \
00096       void getTangents(int i, value_type dots[]) const {    \
00097   expr.getTangents(i, dots); }          \
00098                   \
00099       template <int Arg>            \
00100       value_type getTangent(int i) const {        \
00101   return expr.template getTangent<Arg>(i);      \
00102       }                 \
00103                   \
00104     protected:                \
00105                   \
00106       const ExprT& expr;            \
00107     };                  \
00108                   \
00109     template <typename T>           \
00110     inline Expr< OP< Expr<T> > >          \
00111     OPNAME (const Expr<T>& expr)          \
00112     {                 \
00113       typedef OP< Expr<T> > expr_t;         \
00114                         \
00115       return Expr<expr_t>(expr);          \
00116     }                 \
00117   }                 \
00118 }
00119 
00120 FAD_UNARYOP_MACRO(operator+,
00121       UnaryPlusOp, 
00122       expr.val(),
00123       bar)
00124 FAD_UNARYOP_MACRO(operator-,
00125       UnaryMinusOp, 
00126       -expr.val(),
00127       -bar)
00128 FAD_UNARYOP_MACRO(exp,
00129       ExpOp, 
00130       std::exp(expr.val()),
00131       bar*std::exp(expr.val()))
00132 FAD_UNARYOP_MACRO(log,
00133       LogOp, 
00134       std::log(expr.val()),
00135       bar/expr.val())
00136 FAD_UNARYOP_MACRO(log10,
00137       Log10Op, 
00138       std::log10(expr.val()),
00139       bar/( std::log(value_type(10.))*expr.val() ))
00140 FAD_UNARYOP_MACRO(sqrt,
00141       SqrtOp, 
00142       std::sqrt(expr.val()),
00143       value_type(0.5)*bar/std::sqrt(expr.val()))
00144 FAD_UNARYOP_MACRO(cos,
00145       CosOp, 
00146       std::cos(expr.val()),
00147       -bar*std::sin(expr.val()))
00148 FAD_UNARYOP_MACRO(sin,
00149       SinOp, 
00150       std::sin(expr.val()),
00151       bar*std::cos(expr.val()))
00152 FAD_UNARYOP_MACRO(tan,
00153       TanOp, 
00154       std::tan(expr.val()),
00155       bar*(value_type(1.)+ std::tan(expr.val())*std::tan(expr.val())))
00156 FAD_UNARYOP_MACRO(acos,
00157       ACosOp, 
00158       std::acos(expr.val()),
00159       -bar/std::sqrt(value_type(1.)-expr.val()*expr.val()))
00160 FAD_UNARYOP_MACRO(asin,
00161       ASinOp, 
00162       std::asin(expr.val()),
00163       bar/std::sqrt(value_type(1.)-expr.val()*expr.val()))
00164 FAD_UNARYOP_MACRO(atan,
00165       ATanOp, 
00166       std::atan(expr.val()),
00167       bar/(value_type(1.)+expr.val()*expr.val()))
00168 FAD_UNARYOP_MACRO(cosh,
00169       CoshOp, 
00170       std::cosh(expr.val()),
00171       bar*std::sinh(expr.val()))
00172 FAD_UNARYOP_MACRO(sinh,
00173       SinhOp, 
00174       std::sinh(expr.val()),
00175       bar*std::cosh(expr.val()))
00176 FAD_UNARYOP_MACRO(tanh,
00177       TanhOp, 
00178       std::tanh(expr.val()),
00179       bar/(std::cosh(expr.val())*std::cosh(expr.val())))
00180 FAD_UNARYOP_MACRO(acosh,
00181       ACoshOp, 
00182       std::acosh(expr.val()),
00183       bar/std::sqrt((expr.val()-value_type(1.)) * 
00184         (expr.val()+value_type(1.))))
00185 FAD_UNARYOP_MACRO(asinh,
00186       ASinhOp, 
00187       std::asinh(expr.val()),
00188       bar/std::sqrt(value_type(1.)+expr.val()*expr.val()))
00189 FAD_UNARYOP_MACRO(atanh,
00190       ATanhOp, 
00191       std::atanh(expr.val()),
00192       bar/(value_type(1.)-expr.val()*expr.val()))
00193 FAD_UNARYOP_MACRO(abs,
00194       AbsOp, 
00195       std::abs(expr.val()),
00196       (expr.val() >= value_type(0.)) ? bar : value_type(-bar))
00197 FAD_UNARYOP_MACRO(fabs,
00198       FAbsOp, 
00199       std::fabs(expr.val()),
00200       (expr.val() >= value_type(0.)) ? bar : value_type(-bar))
00201 
00202 #undef FAD_UNARYOP_MACRO
00203 
00204 #define FAD_BINARYOP_MACRO(OPNAME,OP,VALUE,LADJOINT,RADJOINT)   \
00205 namespace Sacado {              \
00206   namespace ELRFad {              \
00207                   \
00208     template <typename ExprT1, typename ExprT2>       \
00209     class OP {};              \
00210                   \
00211     template <typename ExprT1, typename ExprT2>       \
00212     class Expr< OP<ExprT1,ExprT2> > {         \
00213                   \
00214     public:               \
00215                   \
00216       typedef typename ExprT1::value_type value_type_1;     \
00217       typedef typename ExprT2::value_type value_type_2;     \
00218       typedef typename Sacado::Promote<value_type_1,      \
00219                value_type_2>::type value_type;  \
00220                   \
00221       typedef typename ExprT1::base_expr_type base_expr_type_1;   \
00222       typedef typename ExprT2::base_expr_type base_expr_type_2;   \
00223       typedef typename ExprPromote<base_expr_type_1,      \
00224            base_expr_type_2>::type base_expr_type; \
00225                   \
00226       static const int num_args1 = ExprT1::num_args;      \
00227       static const int num_args2 = ExprT2::num_args;      \
00228       static const int num_args = num_args1 + num_args2;    \
00229                   \
00230       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :    \
00231   expr1(expr1_), expr2(expr2_) {}         \
00232                   \
00233       int size() const {            \
00234   int sz1 = expr1.size(), sz2 = expr2.size();     \
00235   return sz1 > sz2 ? sz1 : sz2;         \
00236       }                 \
00237                   \
00238       template <int Arg> bool isActive() const {      \
00239   if (Arg < num_args1)            \
00240     return expr1.template isActive<Arg>();      \
00241   else                \
00242     return expr2.template isActive<Arg-num_args1>();    \
00243       }                 \
00244                   \
00245       bool updateValue() const {          \
00246   return expr1.updateValue() && expr2.updateValue();    \
00247       }                 \
00248                   \
00249       value_type val() const {            \
00250   return VALUE;             \
00251       }                 \
00252                   \
00253       void computePartials(const value_type& bar,     \
00254          value_type partials[]) const {   \
00255   if (num_args1 > 0)            \
00256     expr1.computePartials(LADJOINT, partials);      \
00257   if (num_args2 > 0)            \
00258     expr2.computePartials(RADJOINT, partials+num_args1);    \
00259       }                 \
00260                   \
00261       void getTangents(int i, value_type dots[]) const {    \
00262   expr1.getTangents(i, dots);         \
00263   expr2.getTangents(i, dots+num_args1);       \
00264       }                 \
00265                   \
00266       template <int Arg> value_type getTangent(int i) const {   \
00267   if (Arg < num_args1)            \
00268     return expr1.template getTangent<Arg>(i);     \
00269   else                \
00270     return expr2.template getTangent<Arg-num_args1>(i);   \
00271       }                 \
00272                   \
00273     protected:                \
00274                   \
00275       typename ExprConstRef<ExprT1>::type expr1;      \
00276       typename ExprConstRef<ExprT2>::type expr2;      \
00277                   \
00278     };                  \
00279                   \
00280     template <typename T1, typename T2>         \
00281     inline Expr< OP< Expr<T1>, Expr<T2> > >       \
00282     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
00283     {                 \
00284       typedef OP< Expr<T1>, Expr<T2> > expr_t;        \
00285                       \
00286       return Expr<expr_t>(expr1, expr2);        \
00287     }                 \
00288                   \
00289     template <typename T>           \
00290     inline Expr< OP< Expr<T>, Expr<T> > >       \
00291     OPNAME (const Expr<T>& expr1, const Expr<T>& expr2)     \
00292     {                 \
00293       typedef OP< Expr<T>, Expr<T> > expr_t;        \
00294                       \
00295       return Expr<expr_t>(expr1, expr2);        \
00296     }                 \
00297                   \
00298     template <typename T>           \
00299     inline Expr< OP< ConstExpr<typename Expr<T>::value_type>,   \
00300          Expr<T> > >          \
00301     OPNAME (const typename Expr<T>::value_type& c,      \
00302       const Expr<T>& expr)          \
00303     {                 \
00304       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00305       typedef OP< ConstT, Expr<T> > expr_t;       \
00306                   \
00307       return Expr<expr_t>(ConstT(c), expr);       \
00308     }                 \
00309                   \
00310     template <typename T>           \
00311     inline Expr< OP< Expr<T>,           \
00312          ConstExpr<typename Expr<T>::value_type> > >  \
00313     OPNAME (const Expr<T>& expr,          \
00314       const typename Expr<T>::value_type& c)      \
00315     {                 \
00316       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00317       typedef OP< Expr<T>, ConstT > expr_t;       \
00318                   \
00319       return Expr<expr_t>(expr, ConstT(c));       \
00320     }                 \
00321   }                 \
00322 }
00323 
00324 
00325 FAD_BINARYOP_MACRO(operator+,
00326        AdditionOp, 
00327        expr1.val() + expr2.val(),
00328        bar,
00329        bar)
00330 FAD_BINARYOP_MACRO(operator-,
00331        SubtractionOp, 
00332        expr1.val() - expr2.val(),
00333        bar,
00334        -bar)
00335 FAD_BINARYOP_MACRO(operator*,
00336        MultiplicationOp, 
00337        expr1.val() * expr2.val(),
00338        bar*expr2.val(),
00339        bar*expr1.val())
00340 FAD_BINARYOP_MACRO(operator/,
00341        DivisionOp, 
00342        expr1.val() / expr2.val(),
00343        bar/expr2.val(),
00344        -bar*expr1.val()/(expr2.val()*expr2.val()))
00345 FAD_BINARYOP_MACRO(atan2,
00346        Atan2Op,
00347        std::atan2(expr1.val(), expr2.val()),
00348        bar*expr2.val()/
00349        (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00350        -bar*expr1.val()/
00351        (expr1.val()*expr1.val() + expr2.val()*expr2.val()))
00352 FAD_BINARYOP_MACRO(pow,
00353        PowerOp,
00354        std::pow(expr1.val(), expr2.val()),
00355        expr1.val() == 0 ? value_type(0) : value_type(bar*std::pow(expr1.val(),expr2.val())*expr2.val()/expr1.val()),
00356        expr1.val() == 0 ? value_type(0) : value_type(bar*std::pow(expr1.val(),expr2.val())*std::log(expr1.val())))
00357 FAD_BINARYOP_MACRO(max,
00358                    MaxOp,
00359                    std::max(expr1.val(), expr2.val()),
00360                    expr1.val() >= expr2.val() ? bar : value_type(0.),
00361                    expr2.val() > expr1.val() ? bar : value_type(0.))
00362 FAD_BINARYOP_MACRO(min,
00363                    MinOp,
00364                    std::min(expr1.val(), expr2.val()),
00365                    expr1.val() <= expr2.val() ? bar : value_type(0.),
00366                    expr2.val() < expr1.val() ? bar : value_type(0.))
00367 
00368 #undef FAD_BINARYOP_MACRO
00369 
00370 //-------------------------- Relational Operators -----------------------
00371 
00372 #define FAD_RELOP_MACRO(OP)           \
00373 namespace Sacado {              \
00374   namespace ELRFad {              \
00375     template <typename ExprT1, typename ExprT2>       \
00376     inline bool               \
00377     operator OP (const Expr<ExprT1>& expr1,       \
00378      const Expr<ExprT2>& expr2)       \
00379     {                 \
00380       return expr1.val() OP expr2.val();        \
00381     }                 \
00382                   \
00383     template <typename ExprT2>            \
00384     inline bool               \
00385     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00386      const Expr<ExprT2>& expr2)       \
00387     {                 \
00388       return a OP expr2.val();            \
00389     }                 \
00390                   \
00391     template <typename ExprT1>            \
00392     inline bool               \
00393     operator OP (const Expr<ExprT1>& expr1,       \
00394      const typename Expr<ExprT1>::value_type& b)    \
00395     {                 \
00396       return expr1.val() OP b;            \
00397     }                 \
00398   }                 \
00399 }
00400 
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 FAD_RELOP_MACRO(<<=)
00408 FAD_RELOP_MACRO(>>=)
00409 FAD_RELOP_MACRO(&)
00410 FAD_RELOP_MACRO(|)
00411 
00412 #undef FAD_RELOP_MACRO
00413 
00414 namespace Sacado {
00415 
00416   namespace ELRFad {
00417 
00418     template <typename ExprT>
00419     inline bool operator ! (const Expr<ExprT>& expr) 
00420     {
00421       return ! expr.val();
00422     }
00423 
00424   } // namespace ELRFad
00425 
00426 } // namespace Sacado
00427 
00428 //-------------------------- I/O Operators -----------------------
00429 
00430 namespace Sacado {
00431 
00432   namespace ELRFad {
00433 
00434     template <typename ExprT>
00435     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
00436       typedef typename Expr<ExprT>::base_expr_type base_expr_type;
00437       return os << base_expr_type(x);
00438     }
00439 
00440   } // namespace Fad
00441 
00442 } // namespace Sacado
00443 
00444 
00445 #endif // SACADO_FAD_OPS_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines