Sacado Package Browser (Single Doxygen Collection): Sacado_ELRFad

00001 // $Id: Sacado_ELRFad_Ops.hpp,v 1.1.2.3 2009/03/05 19:41:29 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     public:               \
00069                   \
00070       typedef typename ExprT::value_type value_type;      \
00071                   \
00072       OP() {}                           \
00073                   \
00074       static value_type computeValue(const ExprT& expr) {   \
00075   return VALUE;             \
00076       }                 \
00077                   \
00078       static value_type computeAdjoint(const value_type& bar,   \
00079                const ExprT& expr) {   \
00080   return ADJOINT;             \
00081       }                 \
00082     };                  \
00083                   \
00084     template <typename T>           \
00085     inline Expr< UnaryExpr< Expr<T>, OP > >       \
00086     OPNAME (const Expr<T>& expr)          \
00087     {                 \
00088       typedef UnaryExpr< Expr<T>, OP > expr_t;        \
00089                         \
00090       return Expr<expr_t>(expr_t(expr));        \
00091     }                 \
00092   }                 \
00093 }
00094 
00095 FAD_UNARYOP_MACRO(operator+,
00096       UnaryPlusOp, 
00097       expr.val(),
00098       bar)
00099 FAD_UNARYOP_MACRO(operator-,
00100       UnaryMinusOp, 
00101       -expr.val(),
00102       -bar)
00103 FAD_UNARYOP_MACRO(exp,
00104       ExpOp, 
00105       std::exp(expr.val()),
00106       bar*std::exp(expr.val()))
00107 FAD_UNARYOP_MACRO(log,
00108       LogOp, 
00109       std::log(expr.val()),
00110       bar/expr.val())
00111 FAD_UNARYOP_MACRO(log10,
00112       Log10Op, 
00113       std::log10(expr.val()),
00114       bar/( std::log(value_type(10.))*expr.val() ))
00115 FAD_UNARYOP_MACRO(sqrt,
00116       SqrtOp, 
00117       std::sqrt(expr.val()),
00118       value_type(0.5)*bar/std::sqrt(expr.val()))
00119 FAD_UNARYOP_MACRO(cos,
00120       CosOp, 
00121       std::cos(expr.val()),
00122       -bar*std::sin(expr.val()))
00123 FAD_UNARYOP_MACRO(sin,
00124       SinOp, 
00125       std::sin(expr.val()),
00126       bar*std::cos(expr.val()))
00127 FAD_UNARYOP_MACRO(tan,
00128       TanOp, 
00129       std::tan(expr.val()),
00130       bar*(value_type(1.)+ std::tan(expr.val())*std::tan(expr.val())))
00131 FAD_UNARYOP_MACRO(acos,
00132       ACosOp, 
00133       std::acos(expr.val()),
00134       -bar/std::sqrt(value_type(1.)-expr.val()*expr.val()))
00135 FAD_UNARYOP_MACRO(asin,
00136       ASinOp, 
00137       std::asin(expr.val()),
00138       bar/std::sqrt(value_type(1.)-expr.val()*expr.val()))
00139 FAD_UNARYOP_MACRO(atan,
00140       ATanOp, 
00141       std::atan(expr.val()),
00142       bar/(value_type(1.)+expr.val()*expr.val()))
00143 FAD_UNARYOP_MACRO(cosh,
00144       CoshOp, 
00145       std::cosh(expr.val()),
00146       bar*std::sinh(expr.val()))
00147 FAD_UNARYOP_MACRO(sinh,
00148       SinhOp, 
00149       std::sinh(expr.val()),
00150       bar*std::cosh(expr.val()))
00151 FAD_UNARYOP_MACRO(tanh,
00152       TanhOp, 
00153       std::tanh(expr.val()),
00154       bar/(std::cosh(expr.val())*std::cosh(expr.val())))
00155 FAD_UNARYOP_MACRO(acosh,
00156       ACoshOp, 
00157       acosh(expr.val()),
00158       bar/std::sqrt((expr.val()-value_type(1.)) * 
00159         (expr.val()+value_type(1.))))
00160 FAD_UNARYOP_MACRO(asinh,
00161       ASinhOp, 
00162       asinh(expr.val()),
00163       bar/std::sqrt(value_type(1.)+expr.val()*expr.val()))
00164 FAD_UNARYOP_MACRO(atanh,
00165       ATanhOp, 
00166       atanh(expr.val()),
00167       bar/(value_type(1.)-expr.val()*expr.val()))
00168 FAD_UNARYOP_MACRO(abs,
00169       AbsOp, 
00170       std::abs(expr.val()),
00171       expr.val() >= 0 ? bar : -bar)
00172 FAD_UNARYOP_MACRO(fabs,
00173       FAbsOp, 
00174       std::fabs(expr.val()),
00175       expr.val() >= 0 ? bar : -bar)
00176 
00177 #undef FAD_UNARYOP_MACRO
00178 
00179 #define FAD_BINARYOP_MACRO(OPNAME,OP,VALUE,LADJOINT,RADJOINT)   \
00180 namespace Sacado {              \
00181   namespace ELRFad {              \
00182                   \
00183     template <typename ExprT1, typename ExprT2>       \
00184     class OP {                \
00185                   \
00186     public:               \
00187                   \
00188       typedef typename ExprT1::value_type value_type_1;     \
00189       typedef typename ExprT2::value_type value_type_2;     \
00190       typedef typename Sacado::Promote<value_type_1,      \
00191                value_type_2>::type value_type;  \
00192                   \
00193       OP() {}                                             \
00194                   \
00195       static value_type             \
00196       computeValue(const ExprT1& expr1, const ExprT2& expr2) {          \
00197   return VALUE;             \
00198       }                 \
00199                   \
00200       static value_type             \
00201       computeLeftAdjoint(const value_type& bar,       \
00202          const ExprT1& expr1,       \
00203          const ExprT2& expr2) {       \
00204   return LADJOINT;            \
00205       }                 \
00206                   \
00207       static value_type             \
00208       computeRightAdjoint(const value_type& bar,      \
00209         const ExprT1& expr1,        \
00210         const ExprT2& expr2) {      \
00211   return RADJOINT;            \
00212       }                 \
00213     };                  \
00214                   \
00215     template <typename T1, typename T2>         \
00216     inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, OP > >     \
00217     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
00218     {                 \
00219       typedef BinaryExpr< Expr<T1>, Expr<T2>, OP > expr_t;    \
00220                       \
00221       return Expr<expr_t>(expr_t(expr1, expr2));      \
00222     }                 \
00223                   \
00224     template <typename T>           \
00225     inline Expr< BinaryExpr< Expr<T>, Expr<T>, OP > >     \
00226     OPNAME (const Expr<T>& expr1, const Expr<T>& expr2)     \
00227     {                 \
00228       typedef BinaryExpr< Expr<T>, Expr<T>, OP > expr_t;    \
00229                       \
00230       return Expr<expr_t>(expr_t(expr1, expr2));      \
00231     }                 \
00232                   \
00233     template <typename T>           \
00234     inline Expr< BinaryExpr< ConstExpr<typename Expr<T>::value_type>, \
00235            Expr<T>, OP > >        \
00236     OPNAME (const typename Expr<T>::value_type& c,      \
00237       const Expr<T>& expr)          \
00238     {                 \
00239       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00240       typedef BinaryExpr< ConstT, Expr<T>, OP > expr_t;     \
00241                   \
00242       return Expr<expr_t>(expr_t(ConstT(c), expr));     \
00243     }                 \
00244                   \
00245     template <typename T>           \
00246     inline Expr< BinaryExpr< Expr<T>,         \
00247            ConstExpr<typename Expr<T>::value_type>, \
00248            OP > >         \
00249     OPNAME (const Expr<T>& expr,          \
00250       const typename Expr<T>::value_type& c)      \
00251     {                 \
00252       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00253       typedef BinaryExpr< Expr<T>, ConstT, OP > expr_t;     \
00254                   \
00255       return Expr<expr_t>(expr_t(expr, ConstT(c)));     \
00256     }                 \
00257   }                 \
00258 }
00259 
00260 
00261 FAD_BINARYOP_MACRO(operator+,
00262        AdditionOp, 
00263        expr1.val() + expr2.val(),
00264        bar,
00265        bar)
00266 FAD_BINARYOP_MACRO(operator-,
00267        SubtractionOp, 
00268        expr1.val() - expr2.val(),
00269        bar,
00270        -bar)
00271 FAD_BINARYOP_MACRO(operator*,
00272        MultiplicationOp, 
00273        expr1.val() * expr2.val(),
00274        bar*expr2.val(),
00275        bar*expr1.val())
00276 FAD_BINARYOP_MACRO(operator/,
00277        DivisionOp, 
00278        expr1.val() / expr2.val(),
00279        bar/expr2.val(),
00280        -bar*expr1.val()/(expr2.val()*expr2.val()))
00281 FAD_BINARYOP_MACRO(atan2,
00282        Atan2Op,
00283        std::atan2(expr1.val(), expr2.val()),
00284        bar*expr2.val()/
00285        (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00286        bar*expr1.val()/
00287        (expr1.val()*expr1.val() + expr2.val()*expr2.val()))
00288 FAD_BINARYOP_MACRO(pow,
00289        PowerOp,
00290        std::pow(expr1.val(), expr2.val()),
00291        bar*std::pow(expr1.val(),expr2.val())*expr2.val()/
00292        expr1.val(),
00293        bar*std::pow(expr1.val(),expr2.val())*std::log(expr1.val()))
00294 FAD_BINARYOP_MACRO(max,
00295                    MaxOp,
00296                    std::max(expr1.val(), expr2.val()),
00297                    expr1.val() >= expr2.val() ? bar : value_type(0.),
00298                    expr2.val() > expr1.val() ? bar : value_type(0.))
00299 FAD_BINARYOP_MACRO(min,
00300                    MinOp,
00301                    std::min(expr1.val(), expr2.val()),
00302                    expr1.val() <= expr2.val() ? bar : value_type(0.),
00303                    expr2.val() < expr1.val() ? bar : value_type(0.))
00304 
00305 #undef FAD_BINARYOP_MACRO
00306 
00307 //-------------------------- Relational Operators -----------------------
00308 
00309 #define FAD_RELOP_MACRO(OP)           \
00310 namespace Sacado {              \
00311   namespace ELRFad {              \
00312     template <typename ExprT1, typename ExprT2>       \
00313     inline bool               \
00314     operator OP (const Expr<ExprT1>& expr1,       \
00315      const Expr<ExprT2>& expr2)       \
00316     {                 \
00317       return expr1.val() OP expr2.val();        \
00318     }                 \
00319                   \
00320     template <typename ExprT2>            \
00321     inline bool               \
00322     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00323      const Expr<ExprT2>& expr2)       \
00324     {                 \
00325       return a OP expr2.val();            \
00326     }                 \
00327                   \
00328     template <typename ExprT1>            \
00329     inline bool               \
00330     operator OP (const Expr<ExprT1>& expr1,       \
00331      const typename Expr<ExprT1>::value_type& b)    \
00332     {                 \
00333       return expr1.val() OP b;            \
00334     }                 \
00335   }                 \
00336 }
00337 
00338 FAD_RELOP_MACRO(==)
00339 FAD_RELOP_MACRO(!=)
00340 FAD_RELOP_MACRO(<)
00341 FAD_RELOP_MACRO(>)
00342 FAD_RELOP_MACRO(<=)
00343 FAD_RELOP_MACRO(>=)
00344 FAD_RELOP_MACRO(<<=)
00345 FAD_RELOP_MACRO(>>=)
00346 FAD_RELOP_MACRO(&)
00347 FAD_RELOP_MACRO(|)
00348 
00349 #undef FAD_RELOP_MACRO
00350 
00351 namespace Sacado {
00352 
00353   namespace ELRFad {
00354 
00355     template <typename ExprT>
00356     inline bool operator ! (const Expr<ExprT>& expr) 
00357     {
00358       return ! expr.val();
00359     }
00360 
00361   } // namespace ELRFad
00362 
00363 } // namespace Sacado
00364 
00365 //-------------------------- I/O Operators -----------------------
00366 
00367 namespace Sacado {
00368 
00369   namespace ELRFad {
00370 
00371     template <typename ExprT>
00372     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
00373       typedef typename Expr<ExprT>::base_expr_type base_expr_type;
00374       return os << base_expr_type(x);
00375     }
00376 
00377   } // namespace Fad
00378 
00379 } // namespace Sacado
00380 
00381 
00382 #endif // SACADO_FAD_OPS_HPP