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