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

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