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