Sacado Package Browser (Single Doxygen Collection) Version of the Day
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 namespace Sacado {
00063   namespace CacheFad {
00064 
00065     //
00066     // UnaryPlusOp
00067     //
00068 
00069     template <typename ExprT>
00070     class UnaryPlusOp {};
00071 
00072     template <typename ExprT>
00073     class Expr< UnaryPlusOp<ExprT> > {
00074     public:
00075 
00076       typedef typename ExprT::value_type value_type;
00077 
00078       typedef typename ExprT::scalar_type scalar_type;
00079 
00080       Expr(const ExprT& expr_) : expr(expr_)  {}
00081 
00082       int size() const { return expr.size(); }
00083 
00084       bool updateValue() const { return expr.updateValue(); }
00085 
00086       void cache() const {
00087   expr.cache();
00088       }
00089 
00090       value_type val() const {
00091   return expr.val();
00092       }
00093 
00094       bool isLinear() const {
00095         return expr.isLinear();
00096       }
00097 
00098       bool hasFastAccess() const {
00099         return expr.hasFastAccess();
00100       }
00101 
00102       const value_type dx(int i) const {
00103         return expr.dx(i);
00104       }
00105  
00106       const value_type fastAccessDx(int i) const {
00107         return expr.fastAccessDx(i);
00108       }
00109 
00110     protected:
00111 
00112       const ExprT& expr;
00113     };
00114 
00115     template <typename T>
00116     inline Expr< UnaryPlusOp< Expr<T> > >
00117     operator+ (const Expr<T>& expr)
00118     {
00119       typedef UnaryPlusOp< Expr<T> > expr_t;
00120       
00121       return Expr<expr_t>(expr);
00122     }
00123 
00124     //
00125     // UnaryMinusOp
00126     //
00127     template <typename ExprT>
00128     class UnaryMinusOp {};
00129 
00130     template <typename ExprT>
00131     class Expr< UnaryMinusOp<ExprT> > {
00132     public:
00133 
00134       typedef typename ExprT::value_type value_type;
00135 
00136       typedef typename ExprT::scalar_type scalar_type;
00137 
00138       Expr(const ExprT& expr_) : expr(expr_)  {}
00139 
00140       int size() const { return expr.size(); }
00141 
00142       bool updateValue() const { return expr.updateValue(); }
00143 
00144       void cache() const {
00145   expr.cache();
00146       }
00147 
00148       value_type val() const {
00149   return -expr.val();
00150       }
00151 
00152       bool isLinear() const {
00153         return expr.isLinear();
00154       }
00155 
00156       bool hasFastAccess() const {
00157         return expr.hasFastAccess();
00158       }
00159 
00160       const value_type dx(int i) const {
00161         return -expr.dx(i);
00162       }
00163  
00164       const value_type fastAccessDx(int i) const {
00165         return -expr.fastAccessDx(i);
00166       }
00167 
00168     protected:
00169 
00170       const ExprT& expr;
00171     };
00172 
00173     template <typename T>
00174     inline Expr< UnaryMinusOp< Expr<T> > >
00175     operator- (const Expr<T>& expr)
00176     {
00177       typedef UnaryMinusOp< Expr<T> > expr_t;
00178       
00179       return Expr<expr_t>(expr);
00180     }
00181 
00182     //
00183     // AbsOp
00184     //
00185 
00186     template <typename ExprT>
00187     class AbsOp {};
00188 
00189     template <typename ExprT>
00190     class Expr< AbsOp<ExprT> > {
00191     public:
00192 
00193       typedef typename ExprT::value_type value_type;
00194 
00195       typedef typename ExprT::scalar_type scalar_type;
00196 
00197       Expr(const ExprT& expr_) : expr(expr_)  {}
00198 
00199       int size() const { return expr.size(); }
00200 
00201       bool updateValue() const { return expr.updateValue(); }
00202 
00203       void cache() const {
00204   expr.cache();
00205   v = expr.val();
00206   v_pos = (v >= 0);
00207       }
00208 
00209       value_type val() const {
00210   return std::abs(v);
00211       }
00212 
00213       bool isLinear() const {
00214         return false;
00215       }
00216 
00217       bool hasFastAccess() const {
00218         return expr.hasFastAccess();
00219       }
00220 
00221       const value_type dx(int i) const {
00222         return v_pos ? expr.dx(i) : value_type(-expr.dx(i));
00223       }
00224  
00225       const value_type fastAccessDx(int i) const {
00226         return v_pos ? expr.fastAccessDx(i) : value_type(-expr.fastAccessDx(i));
00227       }
00228 
00229     protected:
00230 
00231       const ExprT& expr;
00232       mutable value_type v;
00233       mutable bool v_pos;
00234     };
00235 
00236     template <typename T>
00237     inline Expr< AbsOp< Expr<T> > >
00238     abs (const Expr<T>& expr)
00239     {
00240       typedef AbsOp< Expr<T> > expr_t;
00241       
00242       return Expr<expr_t>(expr);
00243     }
00244 
00245     //
00246     // FAbsOp
00247     //
00248 
00249     template <typename ExprT>
00250     class FAbsOp {};
00251 
00252     template <typename ExprT>
00253     class Expr< FAbsOp<ExprT> > {
00254     public:
00255 
00256       typedef typename ExprT::value_type value_type;
00257 
00258       typedef typename ExprT::scalar_type scalar_type;
00259 
00260       Expr(const ExprT& expr_) : expr(expr_)  {}
00261 
00262       int size() const { return expr.size(); }
00263 
00264       bool updateValue() const { return expr.updateValue(); }
00265 
00266       void cache() const {
00267   expr.cache();
00268   v = expr.val();
00269   v_pos = (v >= 0);
00270       }
00271 
00272       value_type val() const {
00273   return std::fabs(v);
00274       }
00275 
00276       bool isLinear() const {
00277         return false;
00278       }
00279 
00280       bool hasFastAccess() const {
00281         return expr.hasFastAccess();
00282       }
00283 
00284       const value_type dx(int i) const {
00285         return v_pos ? expr.dx(i) : value_type(-expr.dx(i));
00286       }
00287  
00288       const value_type fastAccessDx(int i) const {
00289         return v_pos ? expr.fastAccessDx(i) : value_type(-expr.fastAccessDx(i));
00290       }
00291 
00292     protected:
00293 
00294       const ExprT& expr;
00295       mutable value_type v;
00296       mutable bool v_pos;
00297     };
00298 
00299     template <typename T>
00300     inline Expr< FAbsOp< Expr<T> > >
00301     fabs (const Expr<T>& expr)
00302     {
00303       typedef FAbsOp< Expr<T> > expr_t;
00304       
00305       return Expr<expr_t>(expr);
00306     }
00307 
00308   }
00309 }
00310 
00311 #define FAD_UNARYOP_MACRO(OPNAME,OP,PARTIAL,VALUE)      \
00312 namespace Sacado {              \
00313   namespace CacheFad {              \
00314                   \
00315     template <typename ExprT>           \
00316     class OP {};              \
00317                   \
00318     template <typename ExprT>           \
00319     class Expr< OP<ExprT> > {           \
00320     public:               \
00321                   \
00322       typedef typename ExprT::value_type value_type;      \
00323       typedef typename ExprT::scalar_type scalar_type;      \
00324                   \
00325       Expr(const ExprT& expr_) : expr(expr_)  {}      \
00326                   \
00327       int size() const { return expr.size(); }        \
00328                   \
00329       bool hasFastAccess() const { return expr.hasFastAccess(); } \
00330                   \
00331       bool isPassive() const { return expr.isPassive();}    \
00332                   \
00333       bool updateValue() const { return expr.updateValue(); }   \
00334                   \
00335       void cache() const {            \
00336   expr.cache();             \
00337         v = expr.val();             \
00338   PARTIAL;              \
00339       }                 \
00340                   \
00341       value_type val() const {            \
00342   return VALUE;             \
00343       }                 \
00344                   \
00345       value_type dx(int i) const {          \
00346   return expr.dx(i)*a;            \
00347       }                 \
00348                   \
00349       value_type fastAccessDx(int i) const {        \
00350   return expr.fastAccessDx(i)*a;          \
00351       }                 \
00352                   \
00353     protected:                \
00354                   \
00355       const ExprT& expr;            \
00356       mutable value_type v;           \
00357       mutable value_type a;           \
00358     };                  \
00359                   \
00360     template <typename T>           \
00361     inline Expr< OP< Expr<T> > >          \
00362     OPNAME (const Expr<T>& expr)          \
00363     {                 \
00364       typedef OP< Expr<T> > expr_t;         \
00365                         \
00366       return Expr<expr_t>(expr);          \
00367     }                 \
00368   }                 \
00369 }
00370 
00371 FAD_UNARYOP_MACRO(exp,
00372       ExpOp, 
00373       a = std::exp(v),
00374       a)
00375 FAD_UNARYOP_MACRO(log,
00376       LogOp, 
00377       a=value_type(1)/v,
00378       std::log(v))
00379 FAD_UNARYOP_MACRO(log10,
00380       Log10Op, 
00381       a = value_type(1)/(std::log(value_type(10))*v),
00382       std::log10(v))
00383 FAD_UNARYOP_MACRO(sqrt,
00384       SqrtOp, 
00385       a = value_type(1)/(value_type(2)*std::sqrt(v)),
00386       std::sqrt(v))
00387 FAD_UNARYOP_MACRO(cos,
00388       CosOp, 
00389       a = -std::sin(v),
00390       std::cos(v))
00391 FAD_UNARYOP_MACRO(sin,
00392       SinOp, 
00393       a = std::cos(v),
00394       std::sin(v))
00395 FAD_UNARYOP_MACRO(tan,
00396       TanOp, 
00397       a = value_type(1)+std::tan(v)*std::tan(v),
00398       std::tan(v))
00399 FAD_UNARYOP_MACRO(acos,
00400       ACosOp, 
00401       a = value_type(-1)/std::sqrt(value_type(1)-v*v),
00402       std::acos(v))
00403 FAD_UNARYOP_MACRO(asin,
00404       ASinOp, 
00405       a = value_type(1)/std::sqrt(value_type(1)-v*v),
00406       std::asin(v))
00407 FAD_UNARYOP_MACRO(atan,
00408       ATanOp, 
00409       a = value_type(1)/(value_type(1)+v*v),
00410       std::atan(v))
00411 FAD_UNARYOP_MACRO(cosh,
00412       CoshOp, 
00413       a = std::sinh(v),
00414       std::cosh(v))
00415 FAD_UNARYOP_MACRO(sinh,
00416       SinhOp, 
00417       a = std::cosh(v),
00418       std::sinh(v))
00419 FAD_UNARYOP_MACRO(tanh,
00420       TanhOp, 
00421       a = value_type(1)/(std::cosh(v)*std::cosh(v)),
00422       std::tanh(v))
00423 FAD_UNARYOP_MACRO(acosh,
00424       ACoshOp, 
00425       a = value_type(1)/std::sqrt((v-value_type(1))*(v+value_type(1))),
00426       std::acosh(v))
00427 FAD_UNARYOP_MACRO(asinh,
00428       ASinhOp, 
00429       a = value_type(1)/std::sqrt(value_type(1)+v*v),
00430       std::asinh(v))
00431 FAD_UNARYOP_MACRO(atanh,
00432       ATanhOp, 
00433       a = value_type(1)/(value_type(1)-v*v),
00434       std::atanh(v))
00435 
00436 #undef FAD_UNARYOP_MACRO
00437 
00438 //
00439 // Binary operators
00440 //
00441 namespace Sacado {
00442   namespace CacheFad {
00443 
00444     //
00445     // AdditionOp
00446     //
00447 
00448     template <typename ExprT1, typename ExprT2>
00449     class AdditionOp {};
00450 
00451     template <typename ExprT1, typename ExprT2>
00452     class Expr< AdditionOp<ExprT1,ExprT2> > {
00453 
00454     public:
00455 
00456       typedef typename ExprT1::value_type value_type_1;
00457       typedef typename ExprT2::value_type value_type_2;
00458       typedef typename Sacado::Promote<value_type_1,
00459                value_type_2>::type value_type;
00460 
00461       typedef typename ExprT1::scalar_type scalar_type_1;
00462       typedef typename ExprT2::scalar_type scalar_type_2;
00463       typedef typename Sacado::Promote<scalar_type_1,
00464            scalar_type_2>::type scalar_type;
00465 
00466       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00467   expr1(expr1_), expr2(expr2_) {}
00468 
00469       int size() const {
00470   int sz1 = expr1.size(), sz2 = expr2.size();
00471   return sz1 > sz2 ? sz1 : sz2;
00472       }
00473 
00474       bool updateValue() const {
00475   return expr1.updateValue() && expr2.updateValue();
00476       }
00477 
00478       void cache() const {
00479   expr1.cache();
00480   expr2.cache();
00481       }
00482 
00483       value_type val() const {
00484   return expr1.val()+expr2.val();
00485       }
00486 
00487       bool isLinear() const {
00488         return expr1.isLinear() && expr2.isLinear();
00489       }
00490 
00491       bool hasFastAccess() const {
00492         return expr1.hasFastAccess() && expr2.hasFastAccess();
00493       }
00494 
00495       const value_type dx(int i) const {
00496         return expr1.dx(i) + expr2.dx(i);
00497       }
00498 
00499       const value_type fastAccessDx(int i) const {
00500         return expr1.fastAccessDx(i) + expr2.fastAccessDx(i);
00501       }
00502 
00503     protected:
00504 
00505       const ExprT1& expr1;
00506       const ExprT2& expr2;
00507 
00508     };
00509 
00510     template <typename ExprT1, typename T2>
00511     class Expr< AdditionOp<ExprT1, ConstExpr<T2> > > {
00512 
00513     public:
00514 
00515       typedef ConstExpr<T2> ExprT2;
00516       typedef typename ExprT1::value_type value_type_1;
00517       typedef typename ExprT2::value_type value_type_2;
00518       typedef typename Sacado::Promote<value_type_1,
00519                value_type_2>::type value_type;
00520 
00521       typedef typename ExprT1::scalar_type scalar_type_1;
00522       typedef typename ExprT2::scalar_type scalar_type_2;
00523       typedef typename Sacado::Promote<scalar_type_1,
00524            scalar_type_2>::type scalar_type;
00525 
00526       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00527   expr1(expr1_), expr2(expr2_) {}
00528 
00529       int size() const {
00530   return expr1.size();
00531       }
00532 
00533       bool updateValue() const {
00534   return expr1.updateValue();
00535       }
00536 
00537       void cache() const {
00538   expr1.cache();
00539       }
00540 
00541       value_type val() const {
00542   return expr1.val() + expr2.val();
00543       }
00544 
00545       bool isLinear() const {
00546         return expr1.isLinear();
00547       }
00548 
00549       bool hasFastAccess() const {
00550         return expr1.hasFastAccess();
00551       }
00552 
00553       const value_type dx(int i) const {
00554         return expr1.dx(i);
00555       }
00556 
00557       const value_type fastAccessDx(int i) const {
00558         return expr1.fastAccessDx(i);
00559       }
00560 
00561     protected:
00562 
00563       const ExprT1& expr1;
00564       ExprT2 expr2;
00565 
00566     };
00567 
00568     template <typename T1, typename ExprT2>
00569     class Expr< AdditionOp< ConstExpr<T1>,ExprT2> > {
00570 
00571     public:
00572 
00573       typedef ConstExpr<T1> ExprT1;
00574       typedef typename ExprT1::value_type value_type_1;
00575       typedef typename ExprT2::value_type value_type_2;
00576       typedef typename Sacado::Promote<value_type_1,
00577                value_type_2>::type value_type;
00578 
00579       typedef typename ExprT1::scalar_type scalar_type_1;
00580       typedef typename ExprT2::scalar_type scalar_type_2;
00581       typedef typename Sacado::Promote<scalar_type_1,
00582            scalar_type_2>::type scalar_type;
00583 
00584       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00585   expr1(expr1_), expr2(expr2_) {}
00586 
00587       int size() const {
00588   return expr2.size();
00589       }
00590 
00591       bool updateValue() const {
00592   return expr2.updateValue();
00593       }
00594 
00595       void cache() const {
00596   expr2.cache();
00597       }
00598 
00599       value_type val() const {
00600   return expr1.val() + expr2.val();
00601       }
00602 
00603       bool isLinear() const {
00604         return expr2.isLinear();
00605       }
00606 
00607       bool hasFastAccess() const {
00608         return expr2.hasFastAccess();
00609       }
00610 
00611       const value_type dx(int i) const {
00612         return expr2.dx(i);
00613       }
00614 
00615       const value_type fastAccessDx(int i) const {
00616         return expr2.fastAccessDx(i);
00617       }
00618 
00619     protected:
00620 
00621       ExprT1 expr1;
00622       const ExprT2& expr2;
00623 
00624     };
00625 
00626     //
00627     // SubtractionOp
00628     //
00629 
00630     template <typename ExprT1, typename ExprT2>
00631     class SubtractionOp {};
00632 
00633     template <typename ExprT1, typename ExprT2>
00634     class Expr< SubtractionOp<ExprT1,ExprT2> > {
00635 
00636     public:
00637 
00638       typedef typename ExprT1::value_type value_type_1;
00639       typedef typename ExprT2::value_type value_type_2;
00640       typedef typename Sacado::Promote<value_type_1,
00641                value_type_2>::type value_type;
00642 
00643       typedef typename ExprT1::scalar_type scalar_type_1;
00644       typedef typename ExprT2::scalar_type scalar_type_2;
00645       typedef typename Sacado::Promote<scalar_type_1,
00646            scalar_type_2>::type scalar_type;
00647 
00648       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00649   expr1(expr1_), expr2(expr2_) {}
00650 
00651       int size() const {
00652   int sz1 = expr1.size(), sz2 = expr2.size();
00653   return sz1 > sz2 ? sz1 : sz2;
00654       }
00655 
00656       bool updateValue() const {
00657   return expr1.updateValue() && expr2.updateValue();
00658       }
00659 
00660       void cache() const {
00661   expr1.cache();
00662   expr2.cache();
00663       }
00664 
00665       value_type val() const {
00666   return expr1.val()-expr2.val();
00667       }
00668 
00669       bool isLinear() const {
00670         return expr1.isLinear() && expr2.isLinear();
00671       }
00672 
00673       bool hasFastAccess() const {
00674         return expr1.hasFastAccess() && expr2.hasFastAccess();
00675       }
00676 
00677       const value_type dx(int i) const {
00678         return expr1.dx(i) - expr2.dx(i);
00679       }
00680 
00681       const value_type fastAccessDx(int i) const {
00682         return expr1.fastAccessDx(i) - expr2.fastAccessDx(i);
00683       }
00684 
00685     protected:
00686 
00687       const ExprT1& expr1;
00688       const ExprT2& expr2;
00689 
00690     };
00691 
00692     template <typename ExprT1, typename T2>
00693     class Expr< SubtractionOp<ExprT1, ConstExpr<T2> > > {
00694 
00695     public:
00696 
00697       typedef ConstExpr<T2> ExprT2;
00698       typedef typename ExprT1::value_type value_type_1;
00699       typedef typename ExprT2::value_type value_type_2;
00700       typedef typename Sacado::Promote<value_type_1,
00701                value_type_2>::type value_type;
00702 
00703       typedef typename ExprT1::scalar_type scalar_type_1;
00704       typedef typename ExprT2::scalar_type scalar_type_2;
00705       typedef typename Sacado::Promote<scalar_type_1,
00706            scalar_type_2>::type scalar_type;
00707 
00708       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00709   expr1(expr1_), expr2(expr2_) {}
00710 
00711       int size() const {
00712   return expr1.size();
00713       }
00714 
00715       bool updateValue() const {
00716   return expr1.updateValue();
00717       }
00718 
00719       void cache() const {
00720   expr1.cache();
00721       }
00722 
00723       value_type val() const {
00724   return expr1.val() - expr2.val();
00725       }
00726 
00727       bool isLinear() const {
00728         return expr1.isLinear();
00729       }
00730 
00731       bool hasFastAccess() const {
00732         return expr1.hasFastAccess();
00733       }
00734 
00735       const value_type dx(int i) const {
00736         return expr1.dx(i);
00737       }
00738 
00739       const value_type fastAccessDx(int i) const {
00740         return expr1.fastAccessDx(i);
00741       }
00742 
00743     protected:
00744 
00745       const ExprT1& expr1;
00746       ExprT2 expr2;
00747 
00748     };
00749 
00750     template <typename T1, typename ExprT2>
00751     class Expr< SubtractionOp< ConstExpr<T1>,ExprT2> > {
00752 
00753     public:
00754 
00755       typedef ConstExpr<T1> ExprT1;
00756       typedef typename ExprT1::value_type value_type_1;
00757       typedef typename ExprT2::value_type value_type_2;
00758       typedef typename Sacado::Promote<value_type_1,
00759                value_type_2>::type value_type;
00760 
00761       typedef typename ExprT1::scalar_type scalar_type_1;
00762       typedef typename ExprT2::scalar_type scalar_type_2;
00763       typedef typename Sacado::Promote<scalar_type_1,
00764            scalar_type_2>::type scalar_type;
00765 
00766       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00767   expr1(expr1_), expr2(expr2_) {}
00768 
00769       int size() const {
00770   return expr2.size();
00771       }
00772 
00773       bool updateValue() const {
00774   return expr2.updateValue();
00775       }
00776 
00777       void cache() const {
00778   expr2.cache();
00779       }
00780 
00781       value_type val() const {
00782   return expr1.val() - expr2.val();
00783       }
00784 
00785       bool isLinear() const {
00786         return expr2.isLinear();
00787       }
00788 
00789       bool hasFastAccess() const {
00790         return expr2.hasFastAccess();
00791       }
00792 
00793       const value_type dx(int i) const {
00794         return -expr2.dx(i);
00795       }
00796 
00797       const value_type fastAccessDx(int i) const {
00798         return -expr2.fastAccessDx(i);
00799       }
00800 
00801     protected:
00802 
00803       ExprT1 expr1;
00804       const ExprT2& expr2;
00805 
00806     };
00807 
00808     //
00809     // MultiplicationOp
00810     //
00811 
00812     template <typename ExprT1, typename ExprT2>
00813     class MultiplicationOp {};
00814 
00815     template <typename ExprT1, typename ExprT2>
00816     class Expr< MultiplicationOp<ExprT1,ExprT2> > {
00817 
00818     public:
00819 
00820       typedef typename ExprT1::value_type value_type_1;
00821       typedef typename ExprT2::value_type value_type_2;
00822       typedef typename Sacado::Promote<value_type_1,
00823                value_type_2>::type value_type;
00824 
00825       typedef typename ExprT1::scalar_type scalar_type_1;
00826       typedef typename ExprT2::scalar_type scalar_type_2;
00827       typedef typename Sacado::Promote<scalar_type_1,
00828            scalar_type_2>::type scalar_type;
00829 
00830       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00831   expr1(expr1_), expr2(expr2_) {}
00832 
00833       int size() const {
00834   int sz1 = expr1.size(), sz2 = expr2.size();
00835   return sz1 > sz2 ? sz1 : sz2;
00836       }
00837 
00838       bool updateValue() const {
00839   return expr1.updateValue() && expr2.updateValue();
00840       }
00841 
00842       void cache() const {
00843   expr1.cache();
00844   expr2.cache();
00845         v1 = expr1.val();
00846   v2 = expr2.val();
00847       }
00848 
00849       value_type val() const {
00850   return v1*v2;
00851       }
00852 
00853       bool isLinear() const {
00854         return false;
00855       }
00856 
00857       bool hasFastAccess() const {
00858         return expr1.hasFastAccess() && expr2.hasFastAccess();
00859       }
00860 
00861       const value_type dx(int i) const {
00862         if (expr1.size() > 0 && expr2.size() > 0)
00863     return v1*expr2.dx(i) + expr1.dx(i)*v2;
00864   else if (expr1.size() > 0)
00865     return expr1.dx(i)*v2;
00866   else
00867     return v1*expr2.dx(i);
00868       }
00869 
00870       const value_type fastAccessDx(int i) const {
00871         return v1*expr2.fastAccessDx(i) + expr1.fastAccessDx(i)*v2;
00872       }
00873 
00874     protected:
00875 
00876       const ExprT1& expr1;
00877       const ExprT2& expr2;
00878       mutable value_type_1 v1;
00879       mutable value_type_2 v2;
00880 
00881     };
00882 
00883     template <typename ExprT1, typename T2>
00884     class Expr< MultiplicationOp<ExprT1, ConstExpr<T2> > > {
00885 
00886     public:
00887 
00888       typedef ConstExpr<T2> ExprT2;
00889       typedef typename ExprT1::value_type value_type_1;
00890       typedef typename ExprT2::value_type value_type_2;
00891       typedef typename Sacado::Promote<value_type_1,
00892                value_type_2>::type value_type;
00893 
00894       typedef typename ExprT1::scalar_type scalar_type_1;
00895       typedef typename ExprT2::scalar_type scalar_type_2;
00896       typedef typename Sacado::Promote<scalar_type_1,
00897            scalar_type_2>::type scalar_type;
00898 
00899       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00900   expr1(expr1_), expr2(expr2_) {}
00901 
00902       int size() const {
00903   return expr1.size();
00904       }
00905 
00906       bool updateValue() const {
00907   return expr1.updateValue();
00908       }
00909 
00910       void cache() const {
00911   expr1.cache();
00912       }
00913 
00914       value_type val() const {
00915   return expr1.val()*expr2.val();
00916       }
00917 
00918       bool isLinear() const {
00919         return expr1.isLinear();
00920       }
00921 
00922       bool hasFastAccess() const {
00923         return expr1.hasFastAccess();
00924       }
00925 
00926       const value_type dx(int i) const {
00927         return expr1.dx(i)*expr2.val();
00928       }
00929 
00930       const value_type fastAccessDx(int i) const {
00931         return expr1.fastAccessDx(i)*expr2.val();
00932       }
00933 
00934     protected:
00935 
00936       const ExprT1& expr1;
00937       ExprT2 expr2;
00938 
00939     };
00940 
00941     template <typename T1, typename ExprT2>
00942     class Expr< MultiplicationOp< ConstExpr<T1>,ExprT2> > {
00943 
00944     public:
00945 
00946       typedef ConstExpr<T1> ExprT1;
00947       typedef typename ExprT1::value_type value_type_1;
00948       typedef typename ExprT2::value_type value_type_2;
00949       typedef typename Sacado::Promote<value_type_1,
00950                value_type_2>::type value_type;
00951 
00952       typedef typename ExprT1::scalar_type scalar_type_1;
00953       typedef typename ExprT2::scalar_type scalar_type_2;
00954       typedef typename Sacado::Promote<scalar_type_1,
00955            scalar_type_2>::type scalar_type;
00956 
00957       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
00958   expr1(expr1_), expr2(expr2_) {}
00959 
00960       int size() const {
00961   return expr2.size();
00962       }
00963 
00964       bool updateValue() const {
00965   return expr2.updateValue();
00966       }
00967 
00968       void cache() const {
00969   expr2.cache();
00970       }
00971 
00972       value_type val() const {
00973   return expr1.val()*expr2.val();
00974       }
00975 
00976       bool isLinear() const {
00977         return expr2.isLinear();
00978       }
00979 
00980       bool hasFastAccess() const {
00981         return expr2.hasFastAccess();
00982       }
00983 
00984       const value_type dx(int i) const {
00985         return expr1.val()*expr2.dx(i);
00986       }
00987 
00988       const value_type fastAccessDx(int i) const {
00989         return expr1.val()*expr2.fastAccessDx(i);
00990       }
00991 
00992     protected:
00993 
00994       ExprT1 expr1;
00995       const ExprT2& expr2;
00996 
00997     };
00998 
00999     //
01000     // DivisionOp
01001     //
01002 
01003     template <typename ExprT1, typename ExprT2>
01004     class DivisionOp {};
01005 
01006     template <typename ExprT1, typename ExprT2>
01007     class Expr< DivisionOp<ExprT1,ExprT2> > {
01008 
01009     public:
01010 
01011       typedef typename ExprT1::value_type value_type_1;
01012       typedef typename ExprT2::value_type value_type_2;
01013       typedef typename Sacado::Promote<value_type_1,
01014                value_type_2>::type value_type;
01015 
01016       typedef typename ExprT1::scalar_type scalar_type_1;
01017       typedef typename ExprT2::scalar_type scalar_type_2;
01018       typedef typename Sacado::Promote<scalar_type_1,
01019            scalar_type_2>::type scalar_type;
01020 
01021       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01022   expr1(expr1_), expr2(expr2_) {}
01023 
01024       int size() const {
01025   int sz1 = expr1.size(), sz2 = expr2.size();
01026   return sz1 > sz2 ? sz1 : sz2;
01027       }
01028 
01029       bool updateValue() const {
01030   return expr1.updateValue() && expr2.updateValue();
01031       }
01032 
01033       void cache() const {
01034   expr1.cache();
01035   expr2.cache();
01036         const value_type_1 v1 = expr1.val();
01037   const value_type_2 v2 = expr2.val();
01038   a = value_type(1)/v2;
01039   v = v1*a;
01040   b = -v/v2;
01041       }
01042 
01043       value_type val() const {
01044   return v;
01045       }
01046 
01047       bool isLinear() const {
01048         return false;
01049       }
01050 
01051       bool hasFastAccess() const {
01052         return expr1.hasFastAccess() && expr2.hasFastAccess();
01053       }
01054 
01055       const value_type dx(int i) const {
01056   if (expr1.size() > 0 && expr2.size() > 0)
01057     return expr1.dx(i)*a + expr2.dx(i)*b;
01058   else if (expr1.size() > 0)
01059     return expr1.dx(i)*a;
01060   else
01061     return expr1.val()*b;
01062       }
01063 
01064       const value_type fastAccessDx(int i) const {
01065         return  expr1.fastAccessDx(i)*a + expr2.fastAccessDx(i)*b;
01066       }
01067 
01068     protected:
01069 
01070       const ExprT1& expr1;
01071       const ExprT2& expr2;
01072       mutable value_type v;
01073       mutable value_type a;
01074       mutable value_type b;
01075 
01076     };
01077 
01078     template <typename ExprT1, typename T2>
01079     class Expr< DivisionOp<ExprT1, ConstExpr<T2> > > {
01080 
01081     public:
01082 
01083       typedef ConstExpr<T2> ExprT2;
01084       typedef typename ExprT1::value_type value_type_1;
01085       typedef typename ExprT2::value_type value_type_2;
01086       typedef typename Sacado::Promote<value_type_1,
01087                value_type_2>::type value_type;
01088 
01089       typedef typename ExprT1::scalar_type scalar_type_1;
01090       typedef typename ExprT2::scalar_type scalar_type_2;
01091       typedef typename Sacado::Promote<scalar_type_1,
01092            scalar_type_2>::type scalar_type;
01093 
01094       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01095   expr1(expr1_), expr2(expr2_) {}
01096 
01097       int size() const {
01098   return expr1.size();
01099       }
01100 
01101       bool updateValue() const {
01102   return expr1.updateValue();
01103       }
01104 
01105       void cache() const {
01106   expr1.cache();
01107   const value_type_1 v1 = expr1.val();
01108   a = value_type(1)/expr2.val();
01109   v = v1*a;
01110       }
01111 
01112       value_type val() const {
01113   return v;
01114       }
01115 
01116       bool isLinear() const {
01117         return expr1.isLinear();
01118       }
01119 
01120       bool hasFastAccess() const {
01121         return expr1.hasFastAccess();
01122       }
01123 
01124       const value_type dx(int i) const {
01125         return expr1.dx(i)*a;
01126       }
01127 
01128       const value_type fastAccessDx(int i) const {
01129         return expr1.fastAccessDx(i)*a;
01130       }
01131 
01132     protected:
01133 
01134       const ExprT1& expr1;
01135       ExprT2 expr2;
01136       mutable value_type v;
01137       mutable value_type a;
01138 
01139     };
01140 
01141     template <typename T1, typename ExprT2>
01142     class Expr< DivisionOp< ConstExpr<T1>,ExprT2> > {
01143 
01144     public:
01145 
01146       typedef ConstExpr<T1> ExprT1;
01147       typedef typename ExprT1::value_type value_type_1;
01148       typedef typename ExprT2::value_type value_type_2;
01149       typedef typename Sacado::Promote<value_type_1,
01150                value_type_2>::type value_type;
01151 
01152       typedef typename ExprT1::scalar_type scalar_type_1;
01153       typedef typename ExprT2::scalar_type scalar_type_2;
01154       typedef typename Sacado::Promote<scalar_type_1,
01155            scalar_type_2>::type scalar_type;
01156 
01157       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01158   expr1(expr1_), expr2(expr2_) {}
01159 
01160       int size() const {
01161   return expr2.size();
01162       }
01163 
01164       bool updateValue() const {
01165   return expr2.updateValue();
01166       }
01167 
01168       void cache() const {
01169   expr2.cache();
01170   const value_type_2 v2 = expr2.val();
01171   v = expr1.val()/v2;
01172   b = -v/v2;
01173       }
01174 
01175       value_type val() const {
01176   return v;
01177       }
01178 
01179       bool isLinear() const {
01180         return false;
01181       }
01182 
01183       bool hasFastAccess() const {
01184         return expr2.hasFastAccess();
01185       }
01186 
01187       const value_type dx(int i) const {
01188         return expr2.dx(i)*b;
01189       }
01190 
01191       const value_type fastAccessDx(int i) const {
01192         return expr2.fastAccessDx(i)*b;
01193       }
01194 
01195     protected:
01196 
01197       ExprT1 expr1;
01198       const ExprT2& expr2;
01199       mutable value_type v;
01200       mutable value_type b;
01201 
01202     };
01203 
01204     //
01205     // Atan2Op
01206     //
01207 
01208     template <typename ExprT1, typename ExprT2>
01209     class Atan2Op {};
01210 
01211     template <typename ExprT1, typename ExprT2>
01212     class Expr< Atan2Op<ExprT1,ExprT2> > {
01213 
01214     public:
01215 
01216       typedef typename ExprT1::value_type value_type_1;
01217       typedef typename ExprT2::value_type value_type_2;
01218       typedef typename Sacado::Promote<value_type_1,
01219                value_type_2>::type value_type;
01220 
01221       typedef typename ExprT1::scalar_type scalar_type_1;
01222       typedef typename ExprT2::scalar_type scalar_type_2;
01223       typedef typename Sacado::Promote<scalar_type_1,
01224            scalar_type_2>::type scalar_type;
01225 
01226       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01227   expr1(expr1_), expr2(expr2_) {}
01228 
01229       int size() const {
01230   int sz1 = expr1.size(), sz2 = expr2.size();
01231   return sz1 > sz2 ? sz1 : sz2;
01232       }
01233 
01234       bool updateValue() const {
01235   return expr1.updateValue() && expr2.updateValue();
01236       }
01237 
01238       void cache() const {
01239   expr1.cache();
01240   expr2.cache();
01241         const value_type_1 v1 = expr1.val();
01242   const value_type_2 v2 = expr2.val();
01243   a = value_type(1)/(v1*v1 + v2*v2);
01244   b = -v1*a;
01245   a = v2*a;
01246   v = std::atan2(v1,v2);
01247       }
01248 
01249       value_type val() const {
01250   return v;
01251       }
01252 
01253       bool isLinear() const {
01254         return false;
01255       }
01256 
01257       bool hasFastAccess() const {
01258         return expr1.hasFastAccess() && expr2.hasFastAccess();
01259       }
01260 
01261       const value_type dx(int i) const {
01262         if (expr1.size() > 0 && expr2.size() > 0)
01263     return expr1.dx(i)*a + expr2.dx(i)*b;
01264   else if (expr1.size() > 0)
01265     return expr1.dx(i)*a;
01266   else
01267     return expr1.val()*b;
01268       }
01269 
01270       const value_type fastAccessDx(int i) const {
01271         return expr1.fastAccessDx(i)*a + expr2.fastAccessDx(i)*b;
01272       }
01273 
01274     protected:
01275 
01276       const ExprT1& expr1;
01277       const ExprT2& expr2;
01278       mutable value_type v;
01279       mutable value_type a;
01280       mutable value_type b;
01281 
01282     };
01283 
01284     template <typename ExprT1, typename T2>
01285     class Expr< Atan2Op<ExprT1, ConstExpr<T2> > > {
01286 
01287     public:
01288 
01289       typedef ConstExpr<T2> ExprT2;
01290       typedef typename ExprT1::value_type value_type_1;
01291       typedef typename ExprT2::value_type value_type_2;
01292       typedef typename Sacado::Promote<value_type_1,
01293                value_type_2>::type value_type;
01294 
01295       typedef typename ExprT1::scalar_type scalar_type_1;
01296       typedef typename ExprT2::scalar_type scalar_type_2;
01297       typedef typename Sacado::Promote<scalar_type_1,
01298            scalar_type_2>::type scalar_type;
01299 
01300       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01301   expr1(expr1_), expr2(expr2_) {}
01302 
01303       int size() const {
01304   return expr1.size();
01305       }
01306 
01307       bool updateValue() const {
01308   return expr1.updateValue();
01309       }
01310 
01311       void cache() const {
01312   expr1.cache();
01313   const value_type_1 v1 = expr1.val();
01314   const value_type_2 v2 = expr2.val();
01315   a = v2/(v1*v1 + v2*v2);
01316   v = std::atan2(v1,v2);
01317       }
01318 
01319       value_type val() const {
01320   return v;
01321       }
01322 
01323       bool isLinear() const {
01324         return false;
01325       }
01326 
01327       bool hasFastAccess() const {
01328         return expr1.hasFastAccess();
01329       }
01330 
01331       const value_type dx(int i) const {
01332         return expr1.dx(i)*a;
01333       }
01334 
01335       const value_type fastAccessDx(int i) const {
01336         return expr1.fastAccessDx(i)*a;
01337       }
01338 
01339     protected:
01340 
01341       const ExprT1& expr1;
01342       ExprT2 expr2;
01343       mutable value_type v;
01344       mutable value_type a;
01345 
01346     };
01347 
01348     template <typename T1, typename ExprT2>
01349     class Expr< Atan2Op< ConstExpr<T1>,ExprT2> > {
01350 
01351     public:
01352 
01353       typedef ConstExpr<T1> ExprT1;
01354       typedef typename ExprT1::value_type value_type_1;
01355       typedef typename ExprT2::value_type value_type_2;
01356       typedef typename Sacado::Promote<value_type_1,
01357                value_type_2>::type value_type;
01358 
01359       typedef typename ExprT1::scalar_type scalar_type_1;
01360       typedef typename ExprT2::scalar_type scalar_type_2;
01361       typedef typename Sacado::Promote<scalar_type_1,
01362            scalar_type_2>::type scalar_type;
01363 
01364       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01365   expr1(expr1_), expr2(expr2_) {}
01366 
01367       int size() const {
01368   return expr2.size();
01369       }
01370 
01371       bool updateValue() const {
01372   return expr2.updateValue();
01373       }
01374 
01375       void cache() const {
01376   expr2.cache();
01377   const value_type_1 v1 = expr1.val();
01378   const value_type_2 v2 = expr2.val();
01379   b = -v1/(v1*v1 + v2*v2);
01380   v = std::atan2(v1,v2);
01381       }
01382 
01383       value_type val() const {
01384   return v;
01385       }
01386 
01387       bool isLinear() const {
01388         return false;
01389       }
01390 
01391       bool hasFastAccess() const {
01392         return expr2.hasFastAccess();
01393       }
01394 
01395       const value_type dx(int i) const {
01396         return expr2.dx(i)*b;
01397       }
01398 
01399       const value_type fastAccessDx(int i) const {
01400         return expr2.fastAccessDx(i)*b;
01401       }
01402 
01403     protected:
01404 
01405       ExprT1 expr1;
01406       const ExprT2& expr2;
01407       mutable value_type v;
01408       mutable value_type b;
01409 
01410     };
01411 
01412     //
01413     // PowerOp
01414     //
01415 
01416     template <typename ExprT1, typename ExprT2>
01417     class PowerOp {};
01418 
01419     template <typename ExprT1, typename ExprT2>
01420     class Expr< PowerOp<ExprT1,ExprT2> > {
01421 
01422     public:
01423 
01424       typedef typename ExprT1::value_type value_type_1;
01425       typedef typename ExprT2::value_type value_type_2;
01426       typedef typename Sacado::Promote<value_type_1,
01427                value_type_2>::type value_type;
01428 
01429       typedef typename ExprT1::scalar_type scalar_type_1;
01430       typedef typename ExprT2::scalar_type scalar_type_2;
01431       typedef typename Sacado::Promote<scalar_type_1,
01432            scalar_type_2>::type scalar_type;
01433 
01434       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01435   expr1(expr1_), expr2(expr2_) {}
01436 
01437       int size() const {
01438   int sz1 = expr1.size(), sz2 = expr2.size();
01439   return sz1 > sz2 ? sz1 : sz2;
01440       }
01441 
01442       bool updateValue() const {
01443   return expr1.updateValue() && expr2.updateValue();
01444       }
01445 
01446       void cache() const {
01447   expr1.cache();
01448   expr2.cache();
01449         const value_type_1 v1 = expr1.val();
01450   const value_type_2 v2 = expr2.val();
01451   v = std::pow(v1,v2);
01452   if (v1 == value_type(0)) {
01453     a = value_type(0);
01454     b = value_type(0);
01455   }
01456   else {
01457     a = v*v2/v1;
01458     b = v*std::log(v1);
01459   }
01460       }
01461 
01462       value_type val() const {
01463   return v;
01464       }
01465 
01466       bool isLinear() const {
01467         return false;
01468       }
01469 
01470       bool hasFastAccess() const {
01471         return expr1.hasFastAccess() && expr2.hasFastAccess();
01472       }
01473 
01474       const value_type dx(int i) const {
01475         if (expr1.size() > 0 && expr2.size() > 0)
01476     return expr1.dx(i)*a + expr2.dx(i)*b;
01477   else if (expr1.size() > 0)
01478     return expr1.dx(i)*a;
01479   else
01480     return expr1.val()*b;
01481       }
01482 
01483       const value_type fastAccessDx(int i) const {
01484         return expr1.fastAccessDx(i)*a + expr2.fastAccessDx(i)*b;
01485       }
01486 
01487     protected:
01488 
01489       const ExprT1& expr1;
01490       const ExprT2& expr2;
01491       mutable value_type v;
01492       mutable value_type a;
01493       mutable value_type b;
01494 
01495     };
01496 
01497     template <typename ExprT1, typename T2>
01498     class Expr< PowerOp<ExprT1, ConstExpr<T2> > > {
01499 
01500     public:
01501 
01502       typedef ConstExpr<T2> ExprT2;
01503       typedef typename ExprT1::value_type value_type_1;
01504       typedef typename ExprT2::value_type value_type_2;
01505       typedef typename Sacado::Promote<value_type_1,
01506                value_type_2>::type value_type;
01507 
01508       typedef typename ExprT1::scalar_type scalar_type_1;
01509       typedef typename ExprT2::scalar_type scalar_type_2;
01510       typedef typename Sacado::Promote<scalar_type_1,
01511            scalar_type_2>::type scalar_type;
01512 
01513       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01514   expr1(expr1_), expr2(expr2_) {}
01515 
01516       int size() const {
01517   return expr1.size();
01518       }
01519 
01520       bool updateValue() const {
01521   return expr1.updateValue();
01522       }
01523 
01524       void cache() const {
01525   expr1.cache();
01526   const value_type_1 v1 = expr1.val();
01527   const value_type_2 v2 = expr2.val();
01528   v = std::pow(v1,v2);
01529   if (v1 == value_type(0)) {
01530     a = value_type(0);
01531   }
01532   else {
01533     a = v*v2/v1;
01534   }
01535       }
01536 
01537       value_type val() const {
01538   return v;
01539       }
01540 
01541       bool isLinear() const {
01542         return false;
01543       }
01544 
01545       bool hasFastAccess() const {
01546         return expr1.hasFastAccess();
01547       }
01548 
01549       const value_type dx(int i) const {
01550         return expr1.dx(i)*a;
01551       }
01552 
01553       const value_type fastAccessDx(int i) const {
01554         return expr1.fastAccessDx(i)*a;
01555       }
01556 
01557     protected:
01558 
01559       const ExprT1& expr1;
01560       ExprT2 expr2;
01561       mutable value_type v;
01562       mutable value_type a;
01563 
01564     };
01565 
01566     template <typename T1, typename ExprT2>
01567     class Expr< PowerOp< ConstExpr<T1>,ExprT2> > {
01568 
01569     public:
01570 
01571       typedef ConstExpr<T1> ExprT1;
01572       typedef typename ExprT1::value_type value_type_1;
01573       typedef typename ExprT2::value_type value_type_2;
01574       typedef typename Sacado::Promote<value_type_1,
01575                value_type_2>::type value_type;
01576 
01577       typedef typename ExprT1::scalar_type scalar_type_1;
01578       typedef typename ExprT2::scalar_type scalar_type_2;
01579       typedef typename Sacado::Promote<scalar_type_1,
01580            scalar_type_2>::type scalar_type;
01581 
01582       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01583   expr1(expr1_), expr2(expr2_) {}
01584 
01585       int size() const {
01586   return expr2.size();
01587       }
01588 
01589       bool updateValue() const {
01590   return expr2.updateValue();
01591       }
01592 
01593       void cache() const {
01594   expr2.cache();
01595   const value_type_1 v1 = expr1.val();
01596   const value_type_2 v2 = expr2.val();
01597   v = std::pow(v1,v2);
01598   if (v1 == value_type(0)) {
01599     b = value_type(0);
01600   }
01601   else {
01602     b = v*std::log(v1);
01603   }
01604       }
01605 
01606       value_type val() const {
01607   return v;
01608       }
01609 
01610       bool isLinear() const {
01611         return false;
01612       }
01613 
01614       bool hasFastAccess() const {
01615         return expr2.hasFastAccess();
01616       }
01617 
01618       const value_type dx(int i) const {
01619         return expr2.dx(i)*b;
01620       }
01621 
01622       const value_type fastAccessDx(int i) const {
01623         return expr2.fastAccessDx(i)*b;
01624       }
01625 
01626     protected:
01627 
01628       ExprT1 expr1;
01629       const ExprT2& expr2;
01630       mutable value_type v;
01631       mutable value_type b;
01632 
01633     };
01634 
01635     //
01636     // MaxOp
01637     //
01638 
01639     template <typename ExprT1, typename ExprT2>
01640     class MaxOp {};
01641 
01642     template <typename ExprT1, typename ExprT2>
01643     class Expr< MaxOp<ExprT1,ExprT2> > {
01644 
01645     public:
01646 
01647       typedef typename ExprT1::value_type value_type_1;
01648       typedef typename ExprT2::value_type value_type_2;
01649       typedef typename Sacado::Promote<value_type_1,
01650                value_type_2>::type value_type;
01651 
01652       typedef typename ExprT1::scalar_type scalar_type_1;
01653       typedef typename ExprT2::scalar_type scalar_type_2;
01654       typedef typename Sacado::Promote<scalar_type_1,
01655            scalar_type_2>::type scalar_type;
01656 
01657       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01658   expr1(expr1_), expr2(expr2_) {}
01659 
01660       int size() const {
01661   int sz1 = expr1.size(), sz2 = expr2.size();
01662   return sz1 > sz2 ? sz1 : sz2;
01663       }
01664 
01665       bool updateValue() const {
01666   return expr1.updateValue() && expr2.updateValue();
01667       }
01668 
01669       void cache() const {
01670   expr1.cache();
01671   expr2.cache();
01672         const value_type_1 v1 = expr1.val();
01673   const value_type_2 v2 = expr2.val();
01674   v = std::max(v1,v2);
01675   max_v1 = (v1 >= v2);
01676       }
01677 
01678       value_type val() const {
01679   return v;
01680       }
01681 
01682       bool isLinear() const {
01683         return false;
01684       }
01685 
01686       bool hasFastAccess() const {
01687         return expr1.hasFastAccess() && expr2.hasFastAccess();
01688       }
01689 
01690       const value_type dx(int i) const {
01691         return max_v1 ? expr1.dx(i) : expr2.dx(i);
01692       }
01693 
01694       const value_type fastAccessDx(int i) const {
01695         return max_v1 ? expr1.fastAccessDx(i) : expr2.fastAccessDx(i);
01696       }
01697 
01698     protected:
01699 
01700       const ExprT1& expr1;
01701       const ExprT2& expr2;
01702       mutable value_type v;
01703       mutable bool max_v1;
01704 
01705     };
01706 
01707     template <typename ExprT1, typename T2>
01708     class Expr< MaxOp<ExprT1, ConstExpr<T2> > > {
01709 
01710     public:
01711 
01712       typedef ConstExpr<T2> ExprT2;
01713       typedef typename ExprT1::value_type value_type_1;
01714       typedef typename ExprT2::value_type value_type_2;
01715       typedef typename Sacado::Promote<value_type_1,
01716                value_type_2>::type value_type;
01717 
01718       typedef typename ExprT1::scalar_type scalar_type_1;
01719       typedef typename ExprT2::scalar_type scalar_type_2;
01720       typedef typename Sacado::Promote<scalar_type_1,
01721            scalar_type_2>::type scalar_type;
01722 
01723       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01724   expr1(expr1_), expr2(expr2_) {}
01725 
01726       int size() const {
01727   return expr1.size();
01728       }
01729 
01730       bool updateValue() const {
01731   return expr1.updateValue();
01732       }
01733 
01734       void cache() const {
01735   expr1.cache();
01736   const value_type_1 v1 = expr1.val();
01737   const value_type_2 v2 = expr2.val();
01738   v = std::max(v1,v2);
01739   max_v1 = (v1 >= v2);
01740       }
01741 
01742       value_type val() const {
01743   return v;
01744       }
01745 
01746       bool isLinear() const {
01747         return false;
01748       }
01749 
01750       bool hasFastAccess() const {
01751         return expr1.hasFastAccess();
01752       }
01753 
01754       const value_type dx(int i) const {
01755         return max_v1 ? expr1.dx(i) : value_type(0);
01756       }
01757 
01758       const value_type fastAccessDx(int i) const {
01759         return max_v1 ? expr1.fastAccessDx(i) : value_type(0);
01760       }
01761 
01762     protected:
01763 
01764       const ExprT1& expr1;
01765       ExprT2 expr2;
01766       mutable value_type v;
01767       mutable bool max_v1;
01768 
01769     };
01770 
01771     template <typename T1, typename ExprT2>
01772     class Expr< MaxOp< ConstExpr<T1>,ExprT2> > {
01773 
01774     public:
01775 
01776       typedef ConstExpr<T1> ExprT1;
01777       typedef typename ExprT1::value_type value_type_1;
01778       typedef typename ExprT2::value_type value_type_2;
01779       typedef typename Sacado::Promote<value_type_1,
01780                value_type_2>::type value_type;
01781 
01782       typedef typename ExprT1::scalar_type scalar_type_1;
01783       typedef typename ExprT2::scalar_type scalar_type_2;
01784       typedef typename Sacado::Promote<scalar_type_1,
01785            scalar_type_2>::type scalar_type;
01786 
01787       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01788   expr1(expr1_), expr2(expr2_) {}
01789 
01790       int size() const {
01791   return expr2.size();
01792       }
01793 
01794       bool updateValue() const {
01795   return expr2.updateValue();
01796       }
01797 
01798       void cache() const {
01799   expr2.cache();
01800   const value_type_1 v1 = expr1.val();
01801   const value_type_2 v2 = expr2.val();
01802   v = std::max(v1,v2);
01803   max_v1 = (v1 >= v2);
01804       }
01805 
01806       value_type val() const {
01807   return v;
01808       }
01809 
01810       bool isLinear() const {
01811         return false;
01812       }
01813 
01814       bool hasFastAccess() const {
01815         return expr2.hasFastAccess();
01816       }
01817 
01818       const value_type dx(int i) const {
01819         return max_v1 ? value_type(0) : expr2.dx(i);
01820       }
01821 
01822       const value_type fastAccessDx(int i) const {
01823         return max_v1 ? value_type(0) : expr2.fastAccessDx(i);
01824       }
01825 
01826     protected:
01827 
01828       ExprT1 expr1;
01829       const ExprT2& expr2;
01830       mutable value_type v;
01831       mutable bool max_v1;
01832 
01833     };
01834 
01835     //
01836     // MinOp
01837     //
01838 
01839     template <typename ExprT1, typename ExprT2>
01840     class MinOp {};
01841 
01842     template <typename ExprT1, typename ExprT2>
01843     class Expr< MinOp<ExprT1,ExprT2> > {
01844 
01845     public:
01846 
01847       typedef typename ExprT1::value_type value_type_1;
01848       typedef typename ExprT2::value_type value_type_2;
01849       typedef typename Sacado::Promote<value_type_1,
01850                value_type_2>::type value_type;
01851 
01852       typedef typename ExprT1::scalar_type scalar_type_1;
01853       typedef typename ExprT2::scalar_type scalar_type_2;
01854       typedef typename Sacado::Promote<scalar_type_1,
01855            scalar_type_2>::type scalar_type;
01856 
01857       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01858   expr1(expr1_), expr2(expr2_) {}
01859 
01860       int size() const {
01861   int sz1 = expr1.size(), sz2 = expr2.size();
01862   return sz1 > sz2 ? sz1 : sz2;
01863       }
01864 
01865       bool updateValue() const {
01866   return expr1.updateValue() && expr2.updateValue();
01867       }
01868 
01869       void cache() const {
01870   expr1.cache();
01871   expr2.cache();
01872         const value_type_1 v1 = expr1.val();
01873   const value_type_2 v2 = expr2.val();
01874   v = std::min(v1,v2);
01875   min_v1 = (v1 <= v2);
01876       }
01877 
01878       value_type val() const {
01879   return v;
01880       }
01881 
01882       bool isLinear() const {
01883         return false;
01884       }
01885 
01886       bool hasFastAccess() const {
01887         return expr1.hasFastAccess() && expr2.hasFastAccess();
01888       }
01889 
01890       const value_type dx(int i) const {
01891         return min_v1 ? expr1.dx(i) : expr2.dx(i);
01892       }
01893 
01894       const value_type fastAccessDx(int i) const {
01895         return min_v1 ? expr1.fastAccessDx(i) : expr2.fastAccessDx(i);
01896       }
01897 
01898     protected:
01899 
01900       const ExprT1& expr1;
01901       const ExprT2& expr2;
01902       mutable value_type v;
01903       mutable bool min_v1;
01904 
01905     };
01906 
01907     template <typename ExprT1, typename T2>
01908     class Expr< MinOp<ExprT1, ConstExpr<T2> > > {
01909 
01910     public:
01911 
01912       typedef ConstExpr<T2> ExprT2;
01913       typedef typename ExprT1::value_type value_type_1;
01914       typedef typename ExprT2::value_type value_type_2;
01915       typedef typename Sacado::Promote<value_type_1,
01916                value_type_2>::type value_type;
01917 
01918       typedef typename ExprT1::scalar_type scalar_type_1;
01919       typedef typename ExprT2::scalar_type scalar_type_2;
01920       typedef typename Sacado::Promote<scalar_type_1,
01921            scalar_type_2>::type scalar_type;
01922 
01923       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01924   expr1(expr1_), expr2(expr2_) {}
01925 
01926       int size() const {
01927   return expr1.size();
01928       }
01929 
01930       bool updateValue() const {
01931   return expr1.updateValue();
01932       }
01933 
01934       void cache() const {
01935   expr1.cache();
01936   const value_type_1 v1 = expr1.val();
01937   const value_type_2 v2 = expr2.val();
01938   v = std::min(v1,v2);
01939   min_v1 = (v1 <= v2);
01940       }
01941 
01942       value_type val() const {
01943   return v;
01944       }
01945 
01946       bool isLinear() const {
01947         return false;
01948       }
01949 
01950       bool hasFastAccess() const {
01951         return expr1.hasFastAccess();
01952       }
01953 
01954       const value_type dx(int i) const {
01955         return min_v1 ? expr1.dx(i) : value_type(0);
01956       }
01957 
01958       const value_type fastAccessDx(int i) const {
01959         return min_v1 ? expr1.fastAccessDx(i) : value_type(0);
01960       }
01961 
01962     protected:
01963 
01964       const ExprT1& expr1;
01965       ExprT2 expr2;
01966       mutable value_type v;
01967       mutable bool min_v1;
01968 
01969     };
01970 
01971     template <typename T1, typename ExprT2>
01972     class Expr< MinOp< ConstExpr<T1>,ExprT2> > {
01973 
01974     public:
01975 
01976       typedef ConstExpr<T1> ExprT1;
01977       typedef typename ExprT1::value_type value_type_1;
01978       typedef typename ExprT2::value_type value_type_2;
01979       typedef typename Sacado::Promote<value_type_1,
01980                value_type_2>::type value_type;
01981 
01982       typedef typename ExprT1::scalar_type scalar_type_1;
01983       typedef typename ExprT2::scalar_type scalar_type_2;
01984       typedef typename Sacado::Promote<scalar_type_1,
01985            scalar_type_2>::type scalar_type;
01986 
01987       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
01988   expr1(expr1_), expr2(expr2_) {}
01989 
01990       int size() const {
01991   return expr2.size();
01992       }
01993 
01994       bool updateValue() const {
01995   return expr2.updateValue();
01996       }
01997 
01998       void cache() const {
01999   expr2.cache();
02000   const value_type_1 v1 = expr1.val();
02001   const value_type_2 v2 = expr2.val();
02002   v = std::min(v1,v2);
02003   min_v1 = (v1 <= v2);
02004       }
02005 
02006       value_type val() const {
02007   return v;
02008       }
02009 
02010       bool isLinear() const {
02011         return false;
02012       }
02013 
02014       bool hasFastAccess() const {
02015         return expr2.hasFastAccess();
02016       }
02017 
02018       const value_type dx(int i) const {
02019         return min_v1 ? value_type(0) : expr2.dx(i);
02020       }
02021 
02022       const value_type fastAccessDx(int i) const {
02023         return min_v1 ? value_type(0) : expr2.fastAccessDx(i);
02024       }
02025 
02026     protected:
02027 
02028       ExprT1 expr1;
02029       const ExprT2& expr2;
02030       mutable value_type v;
02031       mutable bool min_v1;
02032 
02033     };
02034 
02035   }
02036 
02037 }
02038 
02039 #define FAD_BINARYOP_MACRO(OPNAME,OP)         \
02040 namespace Sacado {              \
02041   namespace CacheFad {              \
02042                   \
02043     template <typename T1, typename T2>         \
02044     inline Expr< OP< Expr<T1>, Expr<T2> > >       \
02045     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
02046     {                 \
02047       typedef OP< Expr<T1>, Expr<T2> > expr_t;        \
02048                       \
02049       return Expr<expr_t>(expr1, expr2);        \
02050     }                 \
02051                   \
02052     template <typename T>           \
02053     inline Expr< OP< Expr<T>, Expr<T> > >       \
02054     OPNAME (const Expr<T>& expr1, const Expr<T>& expr2)     \
02055     {                 \
02056       typedef OP< Expr<T>, Expr<T> > expr_t;        \
02057                       \
02058       return Expr<expr_t>(expr1, expr2);        \
02059     }                 \
02060                   \
02061     template <typename T>           \
02062     inline Expr< OP< ConstExpr<typename Expr<T>::value_type>,   \
02063          Expr<T> > >          \
02064     OPNAME (const typename Expr<T>::value_type& c,      \
02065       const Expr<T>& expr)          \
02066     {                 \
02067       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
02068       typedef OP< ConstT, Expr<T> > expr_t;       \
02069                   \
02070       return Expr<expr_t>(ConstT(c), expr);       \
02071     }                 \
02072                   \
02073     template <typename T>           \
02074     inline Expr< OP< Expr<T>,           \
02075          ConstExpr<typename Expr<T>::value_type> > >  \
02076     OPNAME (const Expr<T>& expr,          \
02077       const typename Expr<T>::value_type& c)      \
02078     {                 \
02079       typedef ConstExpr<typename Expr<T>::value_type> ConstT;   \
02080       typedef OP< Expr<T>, ConstT > expr_t;       \
02081                   \
02082       return Expr<expr_t>(expr, ConstT(c));       \
02083     }                 \
02084   }                 \
02085 }
02086 
02087 
02088 FAD_BINARYOP_MACRO(operator+, AdditionOp)
02089 FAD_BINARYOP_MACRO(operator-, SubtractionOp)
02090 FAD_BINARYOP_MACRO(operator*, MultiplicationOp)
02091 FAD_BINARYOP_MACRO(operator/, DivisionOp)
02092 FAD_BINARYOP_MACRO(atan2, Atan2Op)
02093 FAD_BINARYOP_MACRO(pow, PowerOp)
02094 FAD_BINARYOP_MACRO(max, MaxOp)
02095 FAD_BINARYOP_MACRO(min, MinOp)
02096 
02097 #undef FAD_BINARYOP_MACRO
02098 
02099 //-------------------------- Relational Operators -----------------------
02100 
02101 #define FAD_RELOP_MACRO(OP)           \
02102 namespace Sacado {              \
02103   namespace CacheFad {              \
02104     template <typename ExprT1, typename ExprT2>       \
02105     inline bool               \
02106     operator OP (const Expr<ExprT1>& expr1,       \
02107      const Expr<ExprT2>& expr2)       \
02108     {                 \
02109       return expr1.val() OP expr2.val();        \
02110     }                 \
02111                   \
02112     template <typename ExprT2>            \
02113     inline bool               \
02114     operator OP (const typename Expr<ExprT2>::value_type& a,    \
02115      const Expr<ExprT2>& expr2)       \
02116     {                 \
02117       return a OP expr2.val();            \
02118     }                 \
02119                   \
02120     template <typename ExprT1>            \
02121     inline bool               \
02122     operator OP (const Expr<ExprT1>& expr1,       \
02123      const typename Expr<ExprT1>::value_type& b)    \
02124     {                 \
02125       return expr1.val() OP b;            \
02126     }                 \
02127   }                 \
02128 }
02129 
02130 FAD_RELOP_MACRO(==)
02131 FAD_RELOP_MACRO(!=)
02132 FAD_RELOP_MACRO(<)
02133 FAD_RELOP_MACRO(>)
02134 FAD_RELOP_MACRO(<=)
02135 FAD_RELOP_MACRO(>=)
02136 FAD_RELOP_MACRO(<<=)
02137 FAD_RELOP_MACRO(>>=)
02138 FAD_RELOP_MACRO(&)
02139 FAD_RELOP_MACRO(|)
02140 
02141 #undef FAD_RELOP_MACRO
02142 
02143 namespace Sacado {
02144 
02145   namespace CacheFad {
02146 
02147     template <typename ExprT>
02148     inline bool operator ! (const Expr<ExprT>& expr) 
02149     {
02150       return ! expr.val();
02151     }
02152 
02153   } // namespace CacheFad
02154 
02155 } // namespace Sacado
02156 
02157 //-------------------------- Boolean Operators -----------------------
02158 namespace Sacado {
02159 
02160   namespace CacheFad {
02161 
02162     template <typename ExprT>
02163     bool toBool(const Expr<ExprT>& x) {
02164       bool is_zero = (x.val() == 0.0);
02165       for (int i=0; i<x.size(); i++)
02166   is_zero = is_zero && (x.dx(i) == 0.0);
02167       return !is_zero;
02168     }
02169 
02170   } // namespace Fad
02171 
02172 } // namespace Sacado
02173 
02174 #define FAD_BOOL_MACRO(OP)            \
02175 namespace Sacado {              \
02176   namespace CacheFad {              \
02177     template <typename ExprT1, typename ExprT2>       \
02178     inline bool               \
02179     operator OP (const Expr<ExprT1>& expr1,       \
02180      const Expr<ExprT2>& expr2)       \
02181     {                 \
02182       return toBool(expr1) OP toBool(expr2);        \
02183     }                 \
02184                   \
02185     template <typename ExprT2>            \
02186     inline bool               \
02187     operator OP (const typename Expr<ExprT2>::value_type& a,    \
02188      const Expr<ExprT2>& expr2)       \
02189     {                 \
02190       return a OP toBool(expr2);          \
02191     }                 \
02192                   \
02193     template <typename ExprT1>            \
02194     inline bool               \
02195     operator OP (const Expr<ExprT1>& expr1,       \
02196      const typename Expr<ExprT1>::value_type& b)    \
02197     {                 \
02198       return toBool(expr1) OP b;          \
02199     }                 \
02200   }                 \
02201 }
02202 
02203 FAD_BOOL_MACRO(&&)
02204 FAD_BOOL_MACRO(||)
02205 
02206 #undef FAD_BOOL_MACRO
02207 
02208 //-------------------------- I/O Operators -----------------------
02209 
02210 namespace Sacado {
02211 
02212   namespace CacheFad {
02213 
02214     template <typename ExprT>
02215     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
02216       os << x.val() << " [";
02217       
02218       for (int i=0; i< x.size(); i++) {
02219         os << " " << x.dx(i);
02220       }
02221 
02222       os << " ]";
02223       return os;
02224     }
02225 
02226   } // namespace CacheFad
02227 
02228 } // namespace Sacado
02229 
02230 
02231 #endif // SACADO_CACHEFAD_OPS_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines