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

Generated on Wed May 12 21:39:33 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7