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                           LINEAR,DX,FASTACCESSDX)     \
00063 namespace Sacado {              \
00064   namespace ELRFad {              \
00065                   \
00066     template <typename ExprT>           \
00067     class OP {};              \
00068                   \
00069     template <typename ExprT>           \
00070     class Expr< OP<ExprT> > {           \
00071     public:               \
00072                   \
00073       typedef typename ExprT::value_type value_type;      \
00074                   \
00075       typedef typename ExprT::base_expr_type base_expr_type;    \
00076                   \
00077       static const int num_args = ExprT::num_args;      \
00078                         \
00079       static const bool is_linear = LINEAR;       \
00080                   \
00081       Expr(const ExprT& expr_) : expr(expr_)  {}      \
00082                   \
00083       int size() const { return expr.size(); }        \
00084                   \
00085       template <int Arg>            \
00086       bool isActive() const { return expr.template isActive<Arg>(); } \
00087                   \
00088       bool isActive2(int j) const { return expr.template isActive2(j); } \
00089                   \
00090       bool updateValue() const { return expr.updateValue(); }   \
00091                   \
00092       value_type val() const {            \
00093   return VALUE;             \
00094       }                 \
00095                   \
00096       void computePartials(const value_type& bar,     \
00097          value_type partials[]) const {   \
00098   expr.computePartials(ADJOINT, partials);      \
00099       }                 \
00100                   \
00101       void getTangents(int i, value_type dots[]) const {    \
00102   expr.getTangents(i, dots); }          \
00103                   \
00104       template <int Arg>            \
00105       const value_type& getTangent(int i) const {     \
00106   return expr.template getTangent<Arg>(i);      \
00107       }                 \
00108                   \
00109       bool isLinear() const {                                           \
00110         return LINEAR;                                                  \
00111       }                                                                 \
00112                                                                         \
00113       bool hasFastAccess() const {                                      \
00114         return expr.hasFastAccess();                                    \
00115       }                                                                 \
00116                                                                         \
00117       const value_type dx(int i) const {                                \
00118         return DX;                                                      \
00119       }                                                                 \
00120                                                                         \
00121       const value_type fastAccessDx(int i) const {                      \
00122         return FASTACCESSDX;                                            \
00123       }                                                                 \
00124                   \
00125       const value_type* getDx(int j) const {        \
00126   return expr.getDx(j);           \
00127       }                 \
00128                   \
00129       const base_expr_type& getArg(int j) const {     \
00130   return expr.getArg(j);            \
00131       }                 \
00132                   \
00133       int numActiveArgs() const {         \
00134   return expr.numActiveArgs();          \
00135       }                 \
00136                   \
00137       void computeActivePartials(const value_type& bar,     \
00138          value_type *partials) const {    \
00139         expr.computePartials(ADJOINT, partials);      \
00140       }                 \
00141                                                                         \
00142     protected:                \
00143                   \
00144       const ExprT& expr;            \
00145     };                  \
00146                   \
00147     template <typename T>           \
00148     inline Expr< OP< Expr<T> > >          \
00149     OPNAME (const Expr<T>& expr)          \
00150     {                 \
00151       typedef OP< Expr<T> > expr_t;         \
00152                         \
00153       return Expr<expr_t>(expr);          \
00154     }                 \
00155   }                 \
00156 }
00157 
00158 FAD_UNARYOP_MACRO(operator+,
00159       UnaryPlusOp, 
00160       expr.val(),
00161       bar,
00162                   true,
00163                   expr.dx(i),
00164                   expr.fastAccessDx(i))
00165 FAD_UNARYOP_MACRO(operator-,
00166       UnaryMinusOp, 
00167       -expr.val(),
00168       -bar,
00169                   true,
00170                   -expr.dx(i),
00171                   -expr.fastAccessDx(i))
00172 FAD_UNARYOP_MACRO(exp,
00173       ExpOp, 
00174       std::exp(expr.val()),
00175       bar*std::exp(expr.val()),
00176                   false,
00177                   std::exp(expr.val())*expr.dx(i),
00178                   std::exp(expr.val())*expr.fastAccessDx(i))
00179 FAD_UNARYOP_MACRO(log,
00180       LogOp, 
00181       std::log(expr.val()),
00182       bar/expr.val(),
00183                   false,
00184                   expr.dx(i)/expr.val(),
00185                   expr.fastAccessDx(i)/expr.val())
00186 FAD_UNARYOP_MACRO(log10,
00187       Log10Op, 
00188       std::log10(expr.val()),
00189       bar/( std::log(value_type(10.))*expr.val() ),
00190                   false,
00191                   expr.dx(i)/( std::log(value_type(10))*expr.val()),
00192                   expr.fastAccessDx(i) / ( std::log(value_type(10))*expr.val()))
00193 FAD_UNARYOP_MACRO(sqrt,
00194       SqrtOp, 
00195       std::sqrt(expr.val()),
00196       value_type(0.5)*bar/std::sqrt(expr.val()),
00197                   false,
00198                   expr.dx(i)/(value_type(2)* std::sqrt(expr.val())),
00199                   expr.fastAccessDx(i)/(value_type(2)* std::sqrt(expr.val())))
00200 FAD_UNARYOP_MACRO(cos,
00201       CosOp, 
00202       std::cos(expr.val()),
00203       -bar*std::sin(expr.val()),
00204                   false,
00205                   -expr.dx(i)* std::sin(expr.val()),
00206                   -expr.fastAccessDx(i)* std::sin(expr.val()))
00207 FAD_UNARYOP_MACRO(sin,
00208       SinOp, 
00209       std::sin(expr.val()),
00210       bar*std::cos(expr.val()),
00211                   false,
00212                   expr.dx(i)* std::cos(expr.val()),
00213                   expr.fastAccessDx(i)* std::cos(expr.val()))
00214 FAD_UNARYOP_MACRO(tan,
00215       TanOp, 
00216       std::tan(expr.val()),
00217       bar*(value_type(1.)+ std::tan(expr.val())*std::tan(expr.val())),
00218                   false,
00219                   expr.dx(i)*
00220                     (value_type(1)+ std::tan(expr.val())* std::tan(expr.val())),
00221                   expr.fastAccessDx(i)*
00222                     (value_type(1)+ std::tan(expr.val())* std::tan(expr.val())))
00223 FAD_UNARYOP_MACRO(acos,
00224       ACosOp, 
00225       std::acos(expr.val()),
00226       -bar/std::sqrt(value_type(1.)-expr.val()*expr.val()),
00227                   false,
00228                   -expr.dx(i)/ std::sqrt(value_type(1)-expr.val()*expr.val()),
00229                   -expr.fastAccessDx(i) /
00230                     std::sqrt(value_type(1)-expr.val()*expr.val()))
00231 FAD_UNARYOP_MACRO(asin,
00232       ASinOp, 
00233       std::asin(expr.val()),
00234       bar/std::sqrt(value_type(1.)-expr.val()*expr.val()),
00235                   false,
00236                   expr.dx(i)/ std::sqrt(value_type(1)-expr.val()*expr.val()),
00237                   expr.fastAccessDx(i) /
00238                     std::sqrt(value_type(1)-expr.val()*expr.val()))
00239 FAD_UNARYOP_MACRO(atan,
00240       ATanOp, 
00241       std::atan(expr.val()),
00242       bar/(value_type(1.)+expr.val()*expr.val()),
00243                   false,
00244                   expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
00245                   expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
00246 FAD_UNARYOP_MACRO(cosh,
00247       CoshOp, 
00248       std::cosh(expr.val()),
00249       bar*std::sinh(expr.val()),
00250                   false,
00251                   expr.dx(i)* std::sinh(expr.val()),
00252                   expr.fastAccessDx(i)* std::sinh(expr.val()))
00253 FAD_UNARYOP_MACRO(sinh,
00254       SinhOp, 
00255       std::sinh(expr.val()),
00256       bar*std::cosh(expr.val()),
00257                   false,
00258                   expr.dx(i)* std::cosh(expr.val()),
00259                   expr.fastAccessDx(i)* std::cosh(expr.val()))
00260 FAD_UNARYOP_MACRO(tanh,
00261       TanhOp, 
00262       std::tanh(expr.val()),
00263       bar/(std::cosh(expr.val())*std::cosh(expr.val())),
00264                   false,
00265                   expr.dx(i)/( std::cosh(expr.val())* std::cosh(expr.val())),
00266                   expr.fastAccessDx(i) /
00267                     ( std::cosh(expr.val())* std::cosh(expr.val())))
00268 FAD_UNARYOP_MACRO(acosh,
00269       ACoshOp, 
00270       std::acosh(expr.val()),
00271       bar/std::sqrt((expr.val()-value_type(1.)) * 
00272         (expr.val()+value_type(1.))),
00273                   false,
00274                   expr.dx(i)/ std::sqrt((expr.val()-value_type(1)) *
00275                                        (expr.val()+value_type(1))),
00276                   expr.fastAccessDx(i)/ std::sqrt((expr.val()-value_type(1)) *
00277                                                  (expr.val()+value_type(1))))
00278 FAD_UNARYOP_MACRO(asinh,
00279       ASinhOp, 
00280       std::asinh(expr.val()),
00281       bar/std::sqrt(value_type(1.)+expr.val()*expr.val()),
00282                   false,
00283                   expr.dx(i)/ std::sqrt(value_type(1)+expr.val()*expr.val()),
00284                   expr.fastAccessDx(i)/ std::sqrt(value_type(1)+
00285                                                  expr.val()*expr.val()))
00286 FAD_UNARYOP_MACRO(atanh,
00287       ATanhOp, 
00288       std::atanh(expr.val()),
00289       bar/(value_type(1.)-expr.val()*expr.val()),
00290                   false,
00291                   expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
00292                   expr.fastAccessDx(i)/(value_type(1)-
00293                                                  expr.val()*expr.val()))
00294 FAD_UNARYOP_MACRO(abs,
00295       AbsOp, 
00296       std::abs(expr.val()),
00297       (expr.val() >= value_type(0.)) ? bar : value_type(-bar),
00298                   false,
00299                   expr.val() >= 0 ? value_type(+expr.dx(i)) :
00300                     value_type(-expr.dx(i)),
00301                   expr.val() >= 0 ? value_type(+expr.fastAccessDx(i)) :
00302                     value_type(-expr.fastAccessDx(i)))
00303 FAD_UNARYOP_MACRO(fabs,
00304       FAbsOp, 
00305       std::fabs(expr.val()),
00306       (expr.val() >= value_type(0.)) ? bar : value_type(-bar),
00307                   false,
00308                   expr.val() >= 0 ? value_type(+expr.dx(i)) :
00309                     value_type(-expr.dx(i)),
00310                   expr.val() >= 0 ? value_type(+expr.fastAccessDx(i)) :
00311                     value_type(-expr.fastAccessDx(i)))
00312 
00313 #undef FAD_UNARYOP_MACRO
00314 
00315 #define FAD_BINARYOP_MACRO(           \
00316   OPNAME,OP,VALUE,LADJOINT,RADJOINT,          \
00317   LINEAR,CONST_LINEAR_1, CONST_LINEAR_2,        \
00318   LINEAR_2,CONST_LINEAR_1_2, CONST_LINEAR_2_2,        \
00319   DX,FASTACCESSDX,CONST_DX_1,CONST_DX_2,        \
00320   CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2)                    \
00321 namespace Sacado {              \
00322   namespace ELRFad {              \
00323                   \
00324     template <typename ExprT1, typename ExprT2>       \
00325     class OP {};              \
00326                   \
00327     template <typename ExprT1, typename ExprT2>       \
00328     class Expr< OP<ExprT1,ExprT2> > {         \
00329                   \
00330     public:               \
00331                   \
00332       typedef typename ExprT1::value_type value_type_1;     \
00333       typedef typename ExprT2::value_type value_type_2;     \
00334       typedef typename Sacado::Promote<value_type_1,      \
00335                value_type_2>::type value_type;  \
00336                   \
00337       typedef typename ExprT1::base_expr_type base_expr_type_1;   \
00338       typedef typename ExprT2::base_expr_type base_expr_type_2;   \
00339       typedef typename ExprPromote<base_expr_type_1,      \
00340            base_expr_type_2>::type base_expr_type; \
00341                   \
00342       static const int num_args1 = ExprT1::num_args;      \
00343       static const int num_args2 = ExprT2::num_args;      \
00344       static const int num_args = num_args1 + num_args2;    \
00345                   \
00346       static const bool is_linear = LINEAR_2;       \
00347                   \
00348       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :    \
00349   expr1(expr1_), expr2(expr2_) {}         \
00350                   \
00351       int size() const {            \
00352   int sz1 = expr1.size(), sz2 = expr2.size();     \
00353   return sz1 > sz2 ? sz1 : sz2;         \
00354       }                 \
00355                   \
00356       template <int Arg> bool isActive() const {      \
00357   if (Arg < num_args1)            \
00358     return expr1.template isActive<Arg>();      \
00359   else                \
00360     return expr2.template isActive<Arg-num_args1>();    \
00361       }                 \
00362                   \
00363       bool isActive2(int j) const {         \
00364   if (j < num_args1)            \
00365     return expr1.template isActive2(j);       \
00366   else                \
00367     return expr2.template isActive2(j);       \
00368       }                 \
00369                   \
00370       bool updateValue() const {          \
00371   return expr1.updateValue() && expr2.updateValue();    \
00372       }                 \
00373                   \
00374       value_type val() const {            \
00375   return VALUE;             \
00376       }                 \
00377                   \
00378       void computePartials(const value_type& bar,     \
00379          value_type partials[]) const {   \
00380   if (num_args1 > 0)            \
00381     expr1.computePartials(LADJOINT, partials);      \
00382   if (num_args2 > 0)            \
00383     expr2.computePartials(RADJOINT, partials+num_args1);    \
00384       }                 \
00385                   \
00386       void getTangents(int i, value_type dots[]) const {    \
00387   expr1.getTangents(i, dots);         \
00388   expr2.getTangents(i, dots+num_args1);       \
00389       }                 \
00390                   \
00391       template <int Arg> const value_type& getTangent(int i) const {  \
00392   if (Arg < num_args1)            \
00393     return expr1.template getTangent<Arg>(i);     \
00394   else                \
00395     return expr2.template getTangent<Arg-num_args1>(i);   \
00396       }                 \
00397                   \
00398       bool isLinear() const {                                           \
00399         return LINEAR;                                                  \
00400       }                                                                 \
00401                                                                         \
00402       bool hasFastAccess() const {                                      \
00403         return expr1.hasFastAccess() && expr2.hasFastAccess();          \
00404       }                                                                 \
00405                                                                         \
00406       const value_type dx(int i) const {                                \
00407         return DX;                                                      \
00408       }                                                                 \
00409                                                                         \
00410       const value_type fastAccessDx(int i) const {                      \
00411         return FASTACCESSDX;                                            \
00412       }                                                                 \
00413                                                                         \
00414       const value_type* getDx(int j) const {        \
00415   if (j < num_args1)            \
00416     return expr1.getDx(j);          \
00417   else                \
00418     return expr2.getDx(j-num_args1);        \
00419       }                 \
00420                   \
00421       const base_expr_type& getArg(int j) const {     \
00422   if (j < num_args1)            \
00423     return expr1.getArg(j);         \
00424   else                \
00425     return expr2.getArg(j-num_args1);       \
00426       }                 \
00427                   \
00428       int numActiveArgs() const {         \
00429   return expr1.numActiveArgs() + expr2.numActiveArgs();   \
00430       }                 \
00431                   \
00432       void computeActivePartials(const value_type& bar,     \
00433          value_type *partials) const {    \
00434   if (expr1.numActiveArgs() > 0)          \
00435     expr1.computePartials(LADJOINT, partials);      \
00436   if (expr2.numActiveArgs() > 0)          \
00437     expr2.computePartials(RADJOINT, partials+expr2.numActiveArgs()); \
00438       }                 \
00439     protected:                \
00440                   \
00441       typename ExprConstRef<ExprT1>::type expr1;      \
00442       typename ExprConstRef<ExprT2>::type expr2;      \
00443                   \
00444     };                  \
00445                                                                         \
00446     template <typename ExprT1, typename T2>       \
00447     class Expr< OP<ExprT1, ConstExpr<T2> > > {        \
00448                   \
00449     public:               \
00450                   \
00451       typedef ConstExpr<T2> ExprT2;                                     \
00452       typedef typename ExprT1::value_type value_type_1;     \
00453       typedef typename ExprT2::value_type value_type_2;     \
00454       typedef typename Sacado::Promote<value_type_1,      \
00455                value_type_2>::type value_type;  \
00456                   \
00457       typedef typename ExprT1::base_expr_type base_expr_type_1;   \
00458       typedef typename ExprT2::base_expr_type base_expr_type_2;   \
00459       typedef typename ExprPromote<base_expr_type_1,      \
00460            base_expr_type_2>::type base_expr_type; \
00461                   \
00462       static const int num_args = ExprT1::num_args;     \
00463                   \
00464       static const bool is_linear = CONST_LINEAR_2_2;     \
00465                   \
00466       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :    \
00467   expr1(expr1_), expr2(expr2_) {}         \
00468                   \
00469       int size() const {            \
00470   return expr1.size();            \
00471       }                 \
00472                   \
00473       template <int Arg> bool isActive() const {      \
00474   return expr1.template isActive<Arg>();        \
00475       }                 \
00476                   \
00477       bool isActive2(int j) const { return expr1.template isActive2(j); } \
00478                   \
00479       bool updateValue() const {          \
00480   return expr1.updateValue();         \
00481       }                 \
00482                   \
00483       value_type val() const {            \
00484   return VALUE;             \
00485       }                 \
00486                   \
00487       void computePartials(const value_type& bar,     \
00488          value_type partials[]) const {   \
00489   expr1.computePartials(LADJOINT, partials);      \
00490       }                 \
00491                   \
00492       void getTangents(int i, value_type dots[]) const {    \
00493   expr1.getTangents(i, dots);         \
00494       }                 \
00495                   \
00496       template <int Arg> const value_type& getTangent(int i) const {  \
00497   return expr1.template getTangent<Arg>(i);     \
00498       }                 \
00499                   \
00500       bool isLinear() const {                                           \
00501         return CONST_LINEAR_2;                                          \
00502       }                                                                 \
00503                                                                         \
00504       bool hasFastAccess() const {                                      \
00505         return expr1.hasFastAccess();                                   \
00506       }                                                                 \
00507                                                                         \
00508       const value_type dx(int i) const {                                \
00509         return CONST_DX_2;                                              \
00510       }                                                                 \
00511                                                                         \
00512       const value_type fastAccessDx(int i) const {                      \
00513         return CONST_FASTACCESSDX_2;                                    \
00514       }                                                                 \
00515                   \
00516       const value_type* getDx(int j) const {        \
00517   return expr1.getDx(j);            \
00518       }                 \
00519                   \
00520       const base_expr_type& getArg(int j) const {     \
00521   return expr1.getArg(j);           \
00522       }                 \
00523                   \
00524       int numActiveArgs() const {         \
00525   return expr1.numActiveArgs();         \
00526       }                 \
00527                   \
00528       void computeActivePartials(const value_type& bar,     \
00529          value_type *partials) const {    \
00530         expr1.computePartials(LADJOINT, partials);      \
00531       }                 \
00532                                                                         \
00533     protected:                \
00534                   \
00535       typename ExprConstRef<ExprT1>::type expr1;      \
00536       typename ExprConstRef<ExprT2>::type expr2;      \
00537                   \
00538     };                  \
00539                                                                         \
00540     template <typename T1, typename ExprT2>       \
00541     class Expr< OP<ConstExpr<T1>,ExprT2> > {        \
00542                   \
00543     public:               \
00544                   \
00545       typedef ConstExpr<T1> ExprT1;                                     \
00546       typedef typename ExprT1::value_type value_type_1;     \
00547       typedef typename ExprT2::value_type value_type_2;     \
00548       typedef typename Sacado::Promote<value_type_1,      \
00549                value_type_2>::type value_type;  \
00550                   \
00551       typedef typename ExprT1::base_expr_type base_expr_type_1;   \
00552       typedef typename ExprT2::base_expr_type base_expr_type_2;   \
00553       typedef typename ExprPromote<base_expr_type_1,      \
00554            base_expr_type_2>::type base_expr_type; \
00555                   \
00556       static const int num_args = ExprT2::num_args;     \
00557                   \
00558       static const bool is_linear = CONST_LINEAR_1_2;     \
00559                   \
00560       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :    \
00561   expr1(expr1_), expr2(expr2_) {}         \
00562                   \
00563       int size() const {            \
00564   return expr2.size();            \
00565       }                 \
00566                   \
00567       template <int Arg> bool isActive() const {      \
00568   return expr2.template isActive<Arg>();        \
00569       }                 \
00570                   \
00571       bool isActive2(int j) const { return expr2.template isActive2(j); } \
00572                   \
00573       bool updateValue() const {          \
00574   return expr2.updateValue();         \
00575       }                 \
00576                   \
00577       value_type val() const {            \
00578   return VALUE;             \
00579       }                 \
00580                   \
00581       void computePartials(const value_type& bar,     \
00582          value_type partials[]) const {   \
00583   expr2.computePartials(RADJOINT, partials);      \
00584       }                 \
00585                   \
00586       void getTangents(int i, value_type dots[]) const {    \
00587   expr2.getTangents(i, dots);         \
00588       }                 \
00589                   \
00590       template <int Arg> const value_type& getTangent(int i) const {  \
00591   return expr2.template getTangent<Arg>(i);     \
00592       }                 \
00593                   \
00594       bool isLinear() const {                                           \
00595         return CONST_LINEAR_1;                                          \
00596       }                                                                 \
00597                                                                         \
00598       bool hasFastAccess() const {                                      \
00599         return expr2.hasFastAccess();                                   \
00600       }                                                                 \
00601                                                                         \
00602       const value_type dx(int i) const {                                \
00603         return CONST_DX_1;                                              \
00604       }                                                                 \
00605                                                                         \
00606       const value_type fastAccessDx(int i) const {                      \
00607         return CONST_FASTACCESSDX_1;                                    \
00608       }                                                                 \
00609                   \
00610       const value_type* getDx(int j) const {        \
00611   return expr2.getDx(j);            \
00612       }                 \
00613                                                                         \
00614                   \
00615       const base_expr_type& getArg(int j) const {     \
00616   return expr2.getArg(j);           \
00617       }                 \
00618                   \
00619       int numActiveArgs() const {         \
00620   return expr2.numActiveArgs();         \
00621       }                 \
00622                   \
00623       void computeActivePartials(const value_type& bar,     \
00624          value_type *partials) const {    \
00625         expr2.computePartials(RADJOINT, partials);      \
00626       }                 \
00627     protected:                \
00628                   \
00629       typename ExprConstRef<ExprT1>::type expr1;      \
00630       typename ExprConstRef<ExprT2>::type expr2;      \
00631                   \
00632     };                  \
00633                   \
00634     template <typename T1, typename T2>         \
00635     inline Expr< OP< Expr<T1>, Expr<T2> > >       \
00636     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
00637     {                 \
00638       typedef OP< Expr<T1>, Expr<T2> > expr_t;        \
00639                       \
00640       return Expr<expr_t>(expr1, expr2);        \
00641     }                 \
00642                   \
00643     template <typename T>           \
00644     inline Expr< OP< Expr<T>, Expr<T> > >       \
00645     OPNAME (const Expr<T>& expr1, const Expr<T>& expr2)     \
00646     {                 \
00647       typedef OP< Expr<T>, Expr<T> > expr_t;        \
00648                       \
00649       return Expr<expr_t>(expr1, expr2);        \
00650     }                 \
00651                   \
00652     template <typename T>           \
00653     inline Expr< OP< ConstExpr<typename Expr<T>::value_type>,   \
00654          Expr<T> > >          \
00655     OPNAME (const typename Expr<T>::value_type& c,      \
00656       const Expr<T>& expr)          \
00657     {                 \
00658       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00659       typedef OP< ConstT, Expr<T> > expr_t;       \
00660                   \
00661       return Expr<expr_t>(ConstT(c), expr);       \
00662     }                 \
00663                   \
00664     template <typename T>           \
00665     inline Expr< OP< Expr<T>,           \
00666          ConstExpr<typename Expr<T>::value_type> > >  \
00667     OPNAME (const Expr<T>& expr,          \
00668       const typename Expr<T>::value_type& c)      \
00669     {                 \
00670       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00671       typedef OP< Expr<T>, ConstT > expr_t;       \
00672                   \
00673       return Expr<expr_t>(expr, ConstT(c));       \
00674     }                 \
00675   }                 \
00676 }
00677 
00678 
00679 FAD_BINARYOP_MACRO(operator+,
00680        AdditionOp, 
00681        expr1.val() + expr2.val(),
00682        bar,
00683        bar,
00684                    expr1.isLinear() && expr2.isLinear(),
00685                    expr2.isLinear(),
00686                    expr1.isLinear(),
00687        ExprT1::is_linear && ExprT2::is_linear,
00688        ExprT2::is_linear,
00689        ExprT1::is_linear,
00690                    expr1.dx(i) + expr2.dx(i),
00691                    expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
00692                    expr2.dx(i),
00693                    expr1.dx(i),
00694                    expr2.fastAccessDx(i),
00695                    expr1.fastAccessDx(i))
00696 FAD_BINARYOP_MACRO(operator-,
00697        SubtractionOp, 
00698        expr1.val() - expr2.val(),
00699        bar,
00700        -bar,
00701                    expr1.isLinear() && expr2.isLinear(),
00702                    expr2.isLinear(),
00703                    expr1.isLinear(),
00704        ExprT1::is_linear && ExprT2::is_linear,
00705        ExprT2::is_linear,
00706        ExprT1::is_linear,
00707                    expr1.dx(i) - expr2.dx(i),
00708                    expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
00709                    -expr2.dx(i),
00710                    expr1.dx(i),
00711                    -expr2.fastAccessDx(i),
00712                    expr1.fastAccessDx(i))
00713 FAD_BINARYOP_MACRO(operator*,
00714        MultiplicationOp, 
00715        expr1.val() * expr2.val(),
00716        bar*expr2.val(),
00717        bar*expr1.val(),
00718                    false,
00719                    expr2.isLinear(),
00720                    expr1.isLinear(),
00721        false,
00722        ExprT2::is_linear,
00723        ExprT1::is_linear,
00724                    expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
00725                    expr1.val()*expr2.fastAccessDx(i) + 
00726                      expr1.fastAccessDx(i)*expr2.val(),
00727                    expr1.val()*expr2.dx(i),
00728                    expr1.dx(i)*expr2.val(),
00729                    expr1.val()*expr2.fastAccessDx(i),
00730                    expr1.fastAccessDx(i)*expr2.val())
00731 FAD_BINARYOP_MACRO(operator/,
00732        DivisionOp, 
00733        expr1.val() / expr2.val(),
00734        bar/expr2.val(),
00735        -bar*expr1.val()/(expr2.val()*expr2.val()),
00736                    false,
00737                    false,
00738                    expr1.isLinear(),
00739        false,
00740        false,
00741        ExprT1::is_linear,
00742                    (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
00743                      (expr2.val()*expr2.val()),
00744                    (expr1.fastAccessDx(i)*expr2.val() -
00745                       expr2.fastAccessDx(i)*expr1.val()) /
00746                       (expr2.val()*expr2.val()),
00747                    -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
00748                    expr1.dx(i)/expr2.val(),
00749                    -expr2.fastAccessDx(i)*expr1.val() / (expr2.val()*expr2.val()),
00750                    expr1.fastAccessDx(i)/expr2.val())
00751 FAD_BINARYOP_MACRO(atan2,
00752        Atan2Op,
00753        std::atan2(expr1.val(), expr2.val()),
00754        bar*expr2.val()/
00755        (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00756        -bar*expr1.val()/
00757        (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00758                    false,
00759                    false,
00760                    false,
00761        false,
00762                    false,
00763                    false,
00764                    (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/                              (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00765                    (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
00766                      (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00767                    (-expr1.val()*expr2.dx(i)) / (expr1.val()*expr1.val() + expr2.val()*expr2.val()),                   
00768                    (expr2.val()*expr1.dx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
00769                    (-expr1.val()*expr2.fastAccessDx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),                   
00770                    (expr2.val()*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()))
00771 FAD_BINARYOP_MACRO(pow,
00772        PowerOp,
00773        std::pow(expr1.val(), expr2.val()),
00774        expr1.val() == value_type(0) ? value_type(0) : value_type(bar*std::pow(expr1.val(),expr2.val())*expr2.val()/expr1.val()),
00775                    expr1.val() == value_type(0) ? value_type(0) : value_type(bar*std::pow(expr1.val(),expr2.val())*std::log(expr1.val())),
00776                    false,
00777                    false,
00778                    false,
00779        false,
00780                    false,
00781                    false,
00782                    expr1.val() == value_type(0) ? value_type(0) : value_type((expr2.dx(i)*std::log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*std::pow(expr1.val(),expr2.val())),
00783                    expr1.val() == value_type(0) ? value_type(0.0) : value_type((expr2.fastAccessDx(i)*std::log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*std::pow(expr1.val(),expr2.val())),
00784                    expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.dx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
00785                    expr1.val() == value_type(0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)/expr1.val()*std::pow(expr1.val(),expr2.val())),
00786                    expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.fastAccessDx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
00787                    expr1.val() == value_type(0) ? value_type(0.0) : value_type(expr2.val()*expr1.fastAccessDx(i)/expr1.val()*std::pow(expr1.val(),expr2.val())))
00788 FAD_BINARYOP_MACRO(max,
00789                    MaxOp,
00790                    std::max(expr1.val(), expr2.val()),
00791                    expr1.val() >= expr2.val() ? bar : value_type(0.),
00792                    expr2.val() > expr1.val() ? bar : value_type(0.),
00793                    expr1.isLinear() && expr2.isLinear(),
00794                    expr2.isLinear(),
00795                    expr1.isLinear(),
00796        ExprT1::is_linear && ExprT2::is_linear,
00797        ExprT2::is_linear,
00798        ExprT1::is_linear,
00799                    expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00800                    expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) :
00801                                                 expr2.fastAccessDx(i),
00802                    expr1.val() >= expr2.val() ? value_type(0) : expr2.dx(i),
00803                    expr1.val() >= expr2.val() ? expr1.dx(i) : value_type(0),
00804                    expr1.val() >= expr2.val() ? value_type(0) : 
00805                                                 expr2.fastAccessDx(i),
00806                    expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) : 
00807                                                 value_type(0))
00808 FAD_BINARYOP_MACRO(min,
00809                    MinOp,
00810                    std::min(expr1.val(), expr2.val()),
00811                    expr1.val() <= expr2.val() ? bar : value_type(0.),
00812                    expr2.val() < expr1.val() ? bar : value_type(0.),
00813                    expr1.isLinear() && expr2.isLinear(),
00814                    expr2.isLinear(),
00815                    expr1.isLinear(),
00816        ExprT1::is_linear && ExprT2::is_linear,
00817        ExprT2::is_linear,
00818        ExprT1::is_linear,
00819                    expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00820                    expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) :
00821                                                 expr2.fastAccessDx(i),
00822                    expr1.val() <= expr2.val() ? value_type(0) : expr2.dx(i),
00823                    expr1.val() <= expr2.val() ? expr1.dx(i) : value_type(0),
00824                    expr1.val() <= expr2.val() ? value_type(0) : 
00825                                                 expr2.fastAccessDx(i),
00826                    expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) : 
00827                                                 value_type(0))
00828 
00829 #undef FAD_BINARYOP_MACRO
00830 
00831 //-------------------------- Relational Operators -----------------------
00832 
00833 #define FAD_RELOP_MACRO(OP)           \
00834 namespace Sacado {              \
00835   namespace ELRFad {              \
00836     template <typename ExprT1, typename ExprT2>       \
00837     inline bool               \
00838     operator OP (const Expr<ExprT1>& expr1,       \
00839      const Expr<ExprT2>& expr2)       \
00840     {                 \
00841       return expr1.val() OP expr2.val();        \
00842     }                 \
00843                   \
00844     template <typename ExprT2>            \
00845     inline bool               \
00846     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00847      const Expr<ExprT2>& expr2)       \
00848     {                 \
00849       return a OP expr2.val();            \
00850     }                 \
00851                   \
00852     template <typename ExprT1>            \
00853     inline bool               \
00854     operator OP (const Expr<ExprT1>& expr1,       \
00855      const typename Expr<ExprT1>::value_type& b)    \
00856     {                 \
00857       return expr1.val() OP b;            \
00858     }                 \
00859   }                 \
00860 }
00861 
00862 FAD_RELOP_MACRO(==)
00863 FAD_RELOP_MACRO(!=)
00864 FAD_RELOP_MACRO(<)
00865 FAD_RELOP_MACRO(>)
00866 FAD_RELOP_MACRO(<=)
00867 FAD_RELOP_MACRO(>=)
00868 FAD_RELOP_MACRO(<<=)
00869 FAD_RELOP_MACRO(>>=)
00870 FAD_RELOP_MACRO(&)
00871 FAD_RELOP_MACRO(|)
00872 
00873 #undef FAD_RELOP_MACRO
00874 
00875 namespace Sacado {
00876 
00877   namespace ELRFad {
00878 
00879     template <typename ExprT>
00880     inline bool operator ! (const Expr<ExprT>& expr) 
00881     {
00882       return ! expr.val();
00883     }
00884 
00885   } // namespace ELRFad
00886 
00887 } // namespace Sacado
00888 
00889 //-------------------------- Boolean Operators -----------------------
00890 namespace Sacado {
00891 
00892   namespace ELRFad {
00893 
00894     template <typename ExprT>
00895     bool toBool(const Expr<ExprT>& x) {
00896       bool is_zero = (x.val() == 0.0);
00897       for (int i=0; i<x.size(); i++)
00898   is_zero = is_zero && (x.dx(i) == 0.0);
00899       return !is_zero;
00900     }
00901 
00902   } // namespace Fad
00903 
00904 } // namespace Sacado
00905 
00906 #define FAD_BOOL_MACRO(OP)            \
00907 namespace Sacado {              \
00908   namespace ELRFad {              \
00909     template <typename ExprT1, typename ExprT2>       \
00910     inline bool               \
00911     operator OP (const Expr<ExprT1>& expr1,       \
00912      const Expr<ExprT2>& expr2)       \
00913     {                 \
00914       return toBool(expr1) OP toBool(expr2);        \
00915     }                 \
00916                   \
00917     template <typename ExprT2>            \
00918     inline bool               \
00919     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00920      const Expr<ExprT2>& expr2)       \
00921     {                 \
00922       return a OP toBool(expr2);          \
00923     }                 \
00924                   \
00925     template <typename ExprT1>            \
00926     inline bool               \
00927     operator OP (const Expr<ExprT1>& expr1,       \
00928      const typename Expr<ExprT1>::value_type& b)    \
00929     {                 \
00930       return toBool(expr1) OP b;          \
00931     }                 \
00932   }                 \
00933 }
00934 
00935 FAD_BOOL_MACRO(&&)
00936 FAD_BOOL_MACRO(||)
00937 
00938 #undef FAD_BOOL_MACRO
00939 
00940 //-------------------------- I/O Operators -----------------------
00941 
00942 namespace Sacado {
00943 
00944   namespace ELRFad {
00945 
00946     template <typename ExprT>
00947     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
00948       typedef typename Expr<ExprT>::base_expr_type base_expr_type;
00949       return os << base_expr_type(x);
00950     }
00951 
00952   } // namespace Fad
00953 
00954 } // namespace Sacado
00955 
00956 
00957 #endif // SACADO_FAD_OPS_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines