Sacado Package Browser (Single Doxygen Collection): Sacado_CacheFad

00001 // $Id: Sacado_CacheFad_Ops.hpp,v 1.14.2.3 2009/03/05 19:41:29 etphipp Exp $ 
00002 // $Source: /space/CVS/Trilinos/packages/sacado/src/Sacado_CacheFad_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_CACHEFAD_OPS_HPP
00055 #define SACADO_CACHEFAD_OPS_HPP
00056 
00057 #include "Sacado_CacheFad_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,PARTIAL,VALUE,DX,FASTACCESSDX)  \
00063 namespace Sacado {              \
00064   namespace CacheFad {              \
00065                   \
00066     template <typename ExprT>           \
00067     class OP {                \
00068     public:               \
00069                   \
00070       typedef typename ExprT::value_type value_type;      \
00071                   \
00072       OP(const ExprT& expr) {}            \
00073                   \
00074       value_type computeValue(const ExprT& expr) const {    \
00075   v = expr.val();             \
00076   PARTIAL;              \
00077   return VALUE;             \
00078       }                 \
00079                   \
00080       value_type computeDx(int i, const ExprT& expr) const {    \
00081   return DX;              \
00082       }                 \
00083                   \
00084       value_type computeFastAccessDx(int i, const ExprT& expr) const {  \
00085   return FASTACCESSDX;            \
00086       }                 \
00087                   \
00088     protected:                \
00089                   \
00090       mutable value_type v;           \
00091       mutable value_type a;           \
00092     };                  \
00093                   \
00094     template <typename T>           \
00095     inline Expr< UnaryExpr< Expr<T>, OP > >       \
00096     OPNAME (const Expr<T>& expr)          \
00097     {                 \
00098       typedef UnaryExpr< Expr<T>, OP > expr_t;        \
00099                         \
00100       return Expr<expr_t>(expr_t(expr));        \
00101     }                 \
00102   }                 \
00103 }
00104 
00105 FAD_UNARYOP_MACRO(operator+,
00106       UnaryPlusOp, 
00107       ;,
00108       v,
00109       expr.dx(i),
00110       expr.fastAccessDx(i))
00111 FAD_UNARYOP_MACRO(operator-,
00112       UnaryMinusOp,
00113       ;, 
00114       -v,
00115       -expr.dx(i),
00116       -expr.fastAccessDx(i))
00117 FAD_UNARYOP_MACRO(exp,
00118       ExpOp, 
00119       a = std::exp(v),
00120       a,
00121       expr.dx(i)*a,
00122       expr.fastAccessDx(i)*a)
00123 FAD_UNARYOP_MACRO(log,
00124       LogOp, 
00125       ;,
00126       std::log(v),
00127       expr.dx(i)/v,
00128       expr.fastAccessDx(i)/v)
00129 FAD_UNARYOP_MACRO(log10,
00130       Log10Op, 
00131       a = std::log(value_type(10))*v,
00132       std::log10(v),
00133       expr.dx(i)/a,
00134       expr.fastAccessDx(i)/a)
00135 FAD_UNARYOP_MACRO(sqrt,
00136       SqrtOp,
00137       a = value_type(2)*std::sqrt(v),
00138       std::sqrt(v),
00139       expr.dx(i)/a,
00140       expr.fastAccessDx(i)/a)
00141 FAD_UNARYOP_MACRO(cos,
00142       CosOp, 
00143       a = std::sin(v),
00144       std::cos(v),
00145       -expr.dx(i)*a,
00146       -expr.fastAccessDx(i)*a)
00147 FAD_UNARYOP_MACRO(sin,
00148       SinOp, 
00149       a = std::cos(v),
00150       std::sin(v),
00151       expr.dx(i)*a,
00152       expr.fastAccessDx(i)*a)
00153 FAD_UNARYOP_MACRO(tan,
00154       TanOp, 
00155       value_type t = std::tan(v); a = value_type(1)+t*t,
00156       t,
00157       expr.dx(i)*a,
00158       expr.fastAccessDx(i)*a)
00159 FAD_UNARYOP_MACRO(acos,
00160       ACosOp, 
00161       a = - std::sqrt(value_type(1)-v*v),
00162       std::acos(v),
00163       expr.dx(i)/a,
00164       expr.fastAccessDx(i)/a)
00165 FAD_UNARYOP_MACRO(asin,
00166       ASinOp, 
00167       a = std::sqrt(value_type(1)-v*v),
00168       std::asin(v),
00169       expr.dx(i)/a,
00170       expr.fastAccessDx(i)/a)
00171 FAD_UNARYOP_MACRO(atan,
00172       ATanOp, 
00173       a = (value_type(1)+v*v),
00174       std::atan(v),
00175       expr.dx(i)/a,
00176       expr.fastAccessDx(i)/a)
00177 FAD_UNARYOP_MACRO(cosh,
00178       CoshOp, 
00179       a = std::sinh(v),
00180       std::cosh(v),
00181       expr.dx(i)*a,
00182       expr.fastAccessDx(i)*a)
00183 FAD_UNARYOP_MACRO(sinh,
00184       SinhOp, 
00185       a = std::cosh(v),
00186       std::sinh(v),
00187       expr.dx(i)*a,
00188       expr.fastAccessDx(i)*a)
00189 FAD_UNARYOP_MACRO(tanh,
00190       TanhOp, 
00191       a = std::cosh(v); a = a*a,
00192       std::tanh(v),
00193       expr.dx(i)/a,
00194       expr.fastAccessDx(i)/a)
00195 FAD_UNARYOP_MACRO(acosh,
00196       ACoshOp, 
00197       a = std::sqrt((v-value_type(1))*(v+value_type(1))),
00198       acosh(v),
00199       expr.dx(i)/a,
00200       expr.fastAccessDx(i)/a)
00201 FAD_UNARYOP_MACRO(asinh,
00202       ASinhOp, 
00203       a = std::sqrt(value_type(1)+v*v),
00204       asinh(v),
00205       expr.dx(i)/a,
00206       expr.fastAccessDx(i)/a)
00207 FAD_UNARYOP_MACRO(atanh,
00208       ATanhOp, 
00209       a = value_type(1)-v*v,
00210       atanh(v),
00211       expr.dx(i)/a,
00212       expr.fastAccessDx(i)/a)
00213 FAD_UNARYOP_MACRO(abs,
00214       AbsOp, 
00215       ;,
00216       std::abs(v),
00217       v >= 0 ? value_type(+expr.dx(i)) : value_type(-expr.dx(i)),
00218       v >= 0 ? value_type(+expr.fastAccessDx(i)) : 
00219         value_type(-expr.fastAccessDx(i)))
00220 FAD_UNARYOP_MACRO(fabs,
00221       FAbsOp, 
00222       ;,
00223       std::fabs(v),
00224       v >= 0 ? value_type(+expr.dx(i)) : value_type(-expr.dx(i)),
00225       v >= 0 ? value_type(+expr.fastAccessDx(i)) : 
00226         value_type(-expr.fastAccessDx(i)))
00227 
00228 #undef FAD_UNARYOP_MACRO
00229 
00230 #define FAD_BINARYOP_MACRO(OPNAME,OP,PARTIAL,VALUE,DX,FASTACCESSDX,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
00231 namespace Sacado {              \
00232   namespace CacheFad {              \
00233                   \
00234     template <typename ExprT1, typename ExprT2>       \
00235     class OP {                \
00236                   \
00237     public:               \
00238                   \
00239       typedef typename ExprT1::value_type value_type_1;     \
00240       typedef typename ExprT2::value_type value_type_2;     \
00241       typedef typename Sacado::Promote<value_type_1,      \
00242                value_type_2>::type value_type;  \
00243                   \
00244       OP(const ExprT1& expr1, const ExprT2& expr2) {}     \
00245                   \
00246       value_type              \
00247       computeValue(const ExprT1& expr1, const ExprT2& expr2) const {  \
00248   v1 = expr1.val();           \
00249   v2 = expr2.val();           \
00250   PARTIAL;              \
00251   return VALUE;             \
00252       }                 \
00253                   \
00254       value_type              \
00255       computeDx(int i, const ExprT1& expr1,       \
00256     const ExprT2& expr2) const {        \
00257   return DX;              \
00258       }                 \
00259                   \
00260       value_type              \
00261       computeFastAccessDx(int i, const ExprT1& expr1,     \
00262         const ExprT2& expr2) const {      \
00263   return FASTACCESSDX;            \
00264       }                 \
00265                   \
00266     protected:                \
00267                   \
00268       mutable value_type_1 v1;            \
00269       mutable value_type_2 v2;            \
00270       mutable value_type a;           \
00271       mutable value_type b;           \
00272     };                  \
00273                   \
00274     template <typename ExprT1>            \
00275     class OP<ExprT1, ConstExpr<typename ExprT1::value_type> > {   \
00276                   \
00277     public:               \
00278                   \
00279       typedef typename ExprT1::value_type value_type;     \
00280       typedef ConstExpr<typename ExprT1::value_type> ExprT2;    \
00281                   \
00282       OP(const ExprT1& expr1, const ExprT2& expr2) {}     \
00283                   \
00284       value_type              \
00285       computeValue(const ExprT1& expr1, const ExprT2& expr2) const {  \
00286   v1 = expr1.val();           \
00287   v2 = expr2.val();           \
00288   PARTIAL;              \
00289   return VALUE;             \
00290       }                 \
00291                   \
00292       value_type              \
00293       computeDx(int i, const ExprT1& expr1,       \
00294     const ExprT2& expr2) const {        \
00295   return CONST_DX_2;            \
00296       }                 \
00297                   \
00298       value_type              \
00299       computeFastAccessDx(int i, const ExprT1& expr1,     \
00300         const ExprT2& expr2) const {      \
00301   return CONST_FASTACCESSDX_2;          \
00302       }                 \
00303                   \
00304     protected:                \
00305                   \
00306       mutable value_type v1;            \
00307       mutable value_type v2;            \
00308       mutable value_type a;           \
00309       mutable value_type b;           \
00310     };                  \
00311                   \
00312     template <typename ExprT2>            \
00313     class OP<ConstExpr<typename ExprT2::value_type>, ExprT2 > {   \
00314                   \
00315     public:               \
00316                   \
00317       typedef typename ExprT2::value_type value_type;     \
00318       typedef ConstExpr<typename ExprT2::value_type> ExprT1;    \
00319                   \
00320       OP(const ExprT1& expr1, const ExprT2& expr2) {}     \
00321                   \
00322       value_type              \
00323       computeValue(const ExprT1& expr1, const ExprT2& expr2) const {  \
00324   v1 = expr1.val();           \
00325   v2 = expr2.val();           \
00326   PARTIAL;              \
00327   return VALUE;             \
00328       }                 \
00329                   \
00330       value_type              \
00331       computeDx(int i, const ExprT1& expr1,       \
00332     const ExprT2& expr2) const {        \
00333   return CONST_DX_1;            \
00334       }                 \
00335                   \
00336       value_type              \
00337       computeFastAccessDx(int i, const ExprT1& expr1,     \
00338         const ExprT2& expr2) const {      \
00339   return CONST_FASTACCESSDX_1;          \
00340       }                 \
00341                   \
00342     protected:                \
00343                   \
00344       mutable value_type v1;            \
00345       mutable value_type v2;            \
00346       mutable value_type a;           \
00347       mutable value_type b;           \
00348     };                  \
00349                   \
00350     template <typename T1, typename T2>         \
00351     inline Expr< BinaryExpr< Expr<T1>, Expr<T2>, OP > >     \
00352     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
00353     {                 \
00354       typedef BinaryExpr< Expr<T1>, Expr<T2>, OP > expr_t;    \
00355                       \
00356       return Expr<expr_t>(expr_t(expr1, expr2));      \
00357     }                 \
00358                   \
00359     template <typename T>           \
00360     inline Expr< BinaryExpr< Expr<T>, Expr<T>, OP > >     \
00361     OPNAME (const Expr<T>& expr1, const Expr<T>& expr2)     \
00362     {                 \
00363       typedef BinaryExpr< Expr<T>, Expr<T>, OP > expr_t;    \
00364                       \
00365       return Expr<expr_t>(expr_t(expr1, expr2));      \
00366     }                 \
00367                   \
00368     template <typename T>           \
00369     inline Expr< BinaryExpr< ConstExpr<typename Expr<T>::value_type>, \
00370            Expr<T>, OP > >        \
00371     OPNAME (const typename Expr<T>::value_type& c,      \
00372       const Expr<T>& expr)          \
00373     {                 \
00374       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00375       typedef BinaryExpr< ConstT, Expr<T>, OP > expr_t;     \
00376                   \
00377       return Expr<expr_t>(expr_t(ConstT(c), expr));     \
00378     }                 \
00379                   \
00380     template <typename T>           \
00381     inline Expr< BinaryExpr< Expr<T>,         \
00382            ConstExpr<typename Expr<T>::value_type>, \
00383            OP > >         \
00384     OPNAME (const Expr<T>& expr,          \
00385       const typename Expr<T>::value_type& c)      \
00386     {                 \
00387       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00388       typedef BinaryExpr< Expr<T>, ConstT, OP > expr_t;     \
00389                   \
00390       return Expr<expr_t>(expr_t(expr, ConstT(c)));     \
00391     }                 \
00392   }                 \
00393 }
00394 
00395 FAD_BINARYOP_MACRO(operator+,
00396        AdditionOp, 
00397        ;,
00398        v1 + v2,
00399        expr1.dx(i) + expr2.dx(i),
00400        expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
00401        expr2.dx(i),
00402        expr1.dx(i),
00403        expr2.fastAccessDx(i),
00404        expr1.fastAccessDx(i))
00405 FAD_BINARYOP_MACRO(operator-,
00406        SubtractionOp, 
00407        ;,
00408        v1 - v2,
00409        expr1.dx(i) - expr2.dx(i),
00410        expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
00411        -expr2.dx(i),
00412        expr1.dx(i),
00413        -expr2.fastAccessDx(i),
00414        expr1.fastAccessDx(i))
00415 FAD_BINARYOP_MACRO(operator*,
00416        MultiplicationOp, 
00417        ;,
00418        v1*v2,
00419        v1*expr2.dx(i) + expr1.dx(i)*v2,
00420        v1*expr2.fastAccessDx(i) + expr1.fastAccessDx(i)*v2,
00421        v1*expr2.dx(i),
00422        expr1.dx(i)*v2,
00423        v1*expr2.fastAccessDx(i),
00424        expr1.fastAccessDx(i)*v2)
00425 FAD_BINARYOP_MACRO(operator/,
00426        DivisionOp, 
00427        value_type c = v1/v2; a = value_type(1)/v2; b = -c/v2,
00428        c,
00429        expr1.dx(i)*a + expr2.dx(i)*b,
00430        expr1.fastAccessDx(i)*a + expr2.fastAccessDx(i)*b,
00431        expr2.dx(i)*b,
00432        expr1.dx(i)*a,
00433        expr2.fastAccessDx(i)*b,
00434        expr1.fastAccessDx(i)*a)
00435 FAD_BINARYOP_MACRO(atan2,
00436        Atan2Op,
00437        a=v1*v1 + v2*v2,
00438        std::atan2(v1,v2),
00439        (expr1.dx(i)*v2 - expr2.dx(i)*v1)/a,
00440        (expr1.fastAccessDx(i)*v2 - expr2.fastAccessDx(i)*v1)/a,
00441        (-expr2.dx(i)*v1)/a,
00442        (expr1.dx(i)*v2)/a,
00443        (-expr2.fastAccessDx(i)*v1)/a,
00444        (expr1.fastAccessDx(i)*v2)/a)
00445 FAD_BINARYOP_MACRO(pow,
00446        PowerOp,
00447        value_type c= std::pow(v1,v2); a=c*v2/v1; b=c*std::log(v1),
00448        c,
00449        expr1.dx(i)*a + expr2.dx(i)*b,
00450        expr1.fastAccessDx(i)*a + expr2.fastAccessDx(i)*b,
00451        expr2.dx(i)*b,
00452        expr1.dx(i)*a,
00453        expr2.fastAccessDx(i)*b,
00454        expr1.fastAccessDx(i)*a)
00455 FAD_BINARYOP_MACRO(max,
00456        MaxOp,
00457        ;,
00458        std::max(expr1.val(), expr2.val()),
00459        expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00460        expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) : 
00461                                     expr2.fastAccessDx(i),
00462        expr1.val() >= expr2.val() ? value_type(0) : expr2.dx(i),
00463        expr1.val() >= expr2.val() ? expr1.dx(i) : value_type(0),
00464        expr1.val() >= expr2.val() ? value_type(0) : 
00465                                     expr2.fastAccessDx(i),
00466        expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) : 
00467                                     value_type(0))
00468 FAD_BINARYOP_MACRO(min,
00469        MinOp,
00470        ;,
00471        std::min(expr1.val(), expr2.val()),
00472        expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00473        expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) : 
00474                                     expr2.fastAccessDx(i),
00475        expr1.val() <= expr2.val() ? value_type(0) : expr2.dx(i),
00476        expr1.val() <= expr2.val() ? expr1.dx(i) : value_type(0),
00477        expr1.val() <= expr2.val() ? value_type(0) : 
00478                                     expr2.fastAccessDx(i),
00479        expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) : 
00480                                     value_type(0))
00481 
00482 #undef FAD_BINARYOP_MACRO
00483   
00484   // The general definition of max/min works for Fad variables too, except
00485   // we need to add a case when the argument types are different.  This 
00486   // can't conflict with the general definition, so we need to use
00487   // Substitution Failure Is Not An Error
00488 #include "Sacado_mpl_disable_if.hpp"
00489 #include "Sacado_mpl_is_same.hpp"
00490 
00491 #define FAD_SFINAE_BINARYOP_MACRO(OPNAME,OP,PARTIAL,VALUE,DX,FASTACCESSDX,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
00492 namespace Sacado {              \
00493   namespace CacheFad {              \
00494                   \
00495     template <typename ExprT1, typename ExprT2>       \
00496     class OP {                \
00497                   \
00498     public:               \
00499                   \
00500       typedef typename ExprT1::value_type value_type_1;     \
00501       typedef typename ExprT2::value_type value_type_2;     \
00502       typedef typename Sacado::Promote<value_type_1,      \
00503                value_type_2>::type value_type;  \
00504                   \
00505       OP(const ExprT1& expr1, const ExprT2& expr2) {}     \
00506                   \
00507       value_type              \
00508       computeValue(const ExprT1& expr1, const ExprT2& expr2) const {  \
00509   v1 = expr1.val();           \
00510   v2 = expr2.val();           \
00511   PARTIAL;              \
00512   return VALUE;             \
00513       }                 \
00514                   \
00515       value_type              \
00516       computeDx(int i, const ExprT1& expr1,       \
00517     const ExprT2& expr2) const {        \
00518   return DX;              \
00519       }                 \
00520                   \
00521       value_type              \
00522       computeFastAccessDx(int i, const ExprT1& expr1,     \
00523         const ExprT2& expr2) const {      \
00524   return FASTACCESSDX;            \
00525       }                 \
00526                   \
00527     protected:                \
00528                   \
00529       mutable value_type_1 v1;            \
00530       mutable value_type_2 v2;            \
00531       mutable value_type a;           \
00532       mutable value_type b;           \
00533     };                  \
00534                   \
00535     template <typename ExprT1>            \
00536     class OP<ExprT1, ConstExpr<typename ExprT1::value_type> > {   \
00537                   \
00538     public:               \
00539                   \
00540       typedef typename ExprT1::value_type value_type;     \
00541       typedef ConstExpr<typename ExprT1::value_type> ExprT2;    \
00542                   \
00543       OP(const ExprT1& expr1, const ExprT2& expr2) {}     \
00544                   \
00545       value_type              \
00546       computeValue(const ExprT1& expr1, const ExprT2& expr2) const {  \
00547   v1 = expr1.val();           \
00548   v2 = expr2.val();           \
00549   PARTIAL;              \
00550   return VALUE;             \
00551       }                 \
00552                   \
00553       value_type              \
00554       computeDx(int i, const ExprT1& expr1,       \
00555     const ExprT2& expr2) const {        \
00556   return CONST_DX_2;            \
00557       }                 \
00558                   \
00559       value_type              \
00560       computeFastAccessDx(int i, const ExprT1& expr1,     \
00561         const ExprT2& expr2) const {      \
00562   return CONST_FASTACCESSDX_2;          \
00563       }                 \
00564                   \
00565     protected:                \
00566                   \
00567       mutable value_type v1;            \
00568       mutable value_type v2;            \
00569       mutable value_type a;           \
00570       mutable value_type b;           \
00571     };                  \
00572                   \
00573     template <typename ExprT2>            \
00574     class OP<ConstExpr<typename ExprT2::value_type>, ExprT2 > {   \
00575                   \
00576     public:               \
00577                   \
00578       typedef typename ExprT2::value_type value_type;     \
00579       typedef ConstExpr<typename ExprT2::value_type> ExprT1;    \
00580                   \
00581       OP(const ExprT1& expr1, const ExprT2& expr2) {}     \
00582                   \
00583       value_type              \
00584       computeValue(const ExprT1& expr1, const ExprT2& expr2) const {  \
00585   v1 = expr1.val();           \
00586   v2 = expr2.val();           \
00587   PARTIAL;              \
00588   return VALUE;             \
00589       }                 \
00590                   \
00591       value_type              \
00592       computeDx(int i, const ExprT1& expr1,       \
00593     const ExprT2& expr2) const {        \
00594   return CONST_DX_1;            \
00595       }                 \
00596                   \
00597       value_type              \
00598       computeFastAccessDx(int i, const ExprT1& expr1,     \
00599         const ExprT2& expr2) const {      \
00600   return CONST_FASTACCESSDX_1;          \
00601       }                 \
00602                   \
00603     protected:                \
00604                   \
00605       mutable value_type v1;            \
00606       mutable value_type v2;            \
00607       mutable value_type a;           \
00608       mutable value_type b;           \
00609     };                  \
00610                   \
00611     template <typename T1, typename T2>         \
00612     inline                                                              \
00613     typename                                                            \
00614     mpl::disable_if< mpl::is_same<T1,T2>,                               \
00615                      Expr<BinaryExpr<Expr<T1>, Expr<T2>, OP> > >::type  \
00616     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
00617     {                 \
00618       typedef BinaryExpr< Expr<T1>, Expr<T2>, OP > expr_t;    \
00619                       \
00620       return Expr<expr_t>(expr_t(expr1, expr2));      \
00621     }                 \
00622                   \
00623     template <typename T>           \
00624     inline Expr< BinaryExpr< ConstExpr<typename Expr<T>::value_type>, \
00625            Expr<T>, OP > >        \
00626     OPNAME (const typename Expr<T>::value_type& c,      \
00627       const Expr<T>& expr)          \
00628     {                 \
00629       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00630       typedef BinaryExpr< ConstT, Expr<T>, OP > expr_t;     \
00631                   \
00632       return Expr<expr_t>(expr_t(ConstT(c), expr));     \
00633     }                 \
00634                   \
00635     template <typename T>           \
00636     inline Expr< BinaryExpr< Expr<T>,         \
00637            ConstExpr<typename Expr<T>::value_type>, \
00638            OP > >         \
00639     OPNAME (const Expr<T>& expr,          \
00640       const typename Expr<T>::value_type& c)      \
00641     {                 \
00642       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
00643       typedef BinaryExpr< Expr<T>, ConstT, OP > expr_t;     \
00644                   \
00645       return Expr<expr_t>(expr_t(expr, ConstT(c)));     \
00646     }                 \
00647   }                 \
00648 }
00649 
00650 // FAD_SFINAE_BINARYOP_MACRO(max,
00651 //       MaxOp,
00652 //       ;,
00653 //       std::max(expr1.val(), expr2.val()),
00654 //       expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00655 //       expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) : 
00656 //                                    expr2.fastAccessDx(i),
00657 //       expr1.val() >= expr2.val() ? value_type(0) : expr2.dx(i),
00658 //       expr1.val() >= expr2.val() ? expr1.dx(i) : value_type(0),
00659 //       expr1.val() >= expr2.val() ? value_type(0) : 
00660 //                                    expr2.fastAccessDx(i),
00661 //       expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) : 
00662 //                                    value_type(0))
00663 // FAD_SFINAE_BINARYOP_MACRO(min,
00664 //       MinOp,
00665 //       ;,
00666 //       std::min(expr1.val(), expr2.val()),
00667 //       expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00668 //       expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) : 
00669 //                                    expr2.fastAccessDx(i),
00670 //       expr1.val() <= expr2.val() ? value_type(0) : expr2.dx(i),
00671 //       expr1.val() <= expr2.val() ? expr1.dx(i) : value_type(0),
00672 //       expr1.val() <= expr2.val() ? value_type(0) : 
00673 //                                    expr2.fastAccessDx(i),
00674 //       expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) : 
00675 //                                    value_type(0))
00676 
00677 #undef FAD_SFINAE_BINARYOP_MACRO
00678 
00679 //-------------------------- Relational Operators -----------------------
00680 
00681 #define FAD_RELOP_MACRO(OP)           \
00682 namespace Sacado {              \
00683   namespace CacheFad {              \
00684     template <typename ExprT1, typename ExprT2>       \
00685     inline bool               \
00686     operator OP (const Expr<ExprT1>& expr1,       \
00687      const Expr<ExprT2>& expr2)       \
00688     {                 \
00689       return expr1.val() OP expr2.val();        \
00690     }                 \
00691                   \
00692     template <typename ExprT2>            \
00693     inline bool               \
00694     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00695      const Expr<ExprT2>& expr2)       \
00696     {                 \
00697       return a OP expr2.val();            \
00698     }                 \
00699                   \
00700     template <typename ExprT1>            \
00701     inline bool               \
00702     operator OP (const Expr<ExprT1>& expr1,       \
00703      const typename Expr<ExprT1>::value_type& b)    \
00704     {                 \
00705       return expr1.val() OP b;            \
00706     }                 \
00707   }                 \
00708 }
00709 
00710 FAD_RELOP_MACRO(==)
00711 FAD_RELOP_MACRO(!=)
00712 FAD_RELOP_MACRO(<)
00713 FAD_RELOP_MACRO(>)
00714 FAD_RELOP_MACRO(<=)
00715 FAD_RELOP_MACRO(>=)
00716 FAD_RELOP_MACRO(<<=)
00717 FAD_RELOP_MACRO(>>=)
00718 FAD_RELOP_MACRO(&)
00719 FAD_RELOP_MACRO(|)
00720 
00721 #undef FAD_RELOP_MACRO
00722 
00723 namespace Sacado {
00724 
00725   namespace CacheFad {
00726 
00727     template <typename ExprT>
00728     inline bool operator ! (const Expr<ExprT>& expr) 
00729     {
00730       return ! expr.val();
00731     }
00732 
00733   } // namespace CacheFad
00734 
00735 } // namespace Sacado
00736 
00737 //-------------------------- I/O Operators -----------------------
00738 
00739 namespace Sacado {
00740 
00741   namespace CacheFad {
00742 
00743     template <typename ExprT>
00744     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
00745       os << x.val() << "  [";
00746       
00747       for (int i=0; i< x.size(); i++) {
00748   os << " " << x.dx(i);
00749       }
00750 
00751       os << " ]\n";
00752       return os;
00753     }
00754 
00755   } // namespace CacheFad
00756 
00757 } // namespace Sacado
00758 
00759 
00760 #endif // SACADO_CACHEFAD_OPS_HPP