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