Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_LFad_LogicalSparseOps.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_LFAD_LOGICALSPARSEOPS_HPP
00055 #define SACADO_LFAD_LOGICALSPARSEOPS_HPP
00056 
00057 #include <cmath>
00058 #include <algorithm>  // for std::min and std::max
00059 #include <ostream>  // for std::ostream
00060 
00061 #define FAD_UNARYOP_MACRO(OPNAME,OP,VALUE)        \
00062 namespace Sacado {              \
00063   namespace LFad {              \
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       typedef typename ExprT::logical_type logical_type;    \
00074                   \
00075       Expr(const ExprT& expr_) : expr(expr_)  {}      \
00076                   \
00077       int size() const { return expr.size(); }        \
00078                   \
00079       bool hasFastAccess() const { return expr.hasFastAccess(); } \
00080                   \
00081       bool isPassive() const { return expr.isPassive();}    \
00082                         \
00083       value_type val() const {            \
00084   return VALUE;             \
00085       }                 \
00086                   \
00087       logical_type dx(int i) const {          \
00088   return expr.dx(i);            \
00089       }                 \
00090                   \
00091       logical_type fastAccessDx(int i) const {        \
00092   return expr.fastAccessDx(i);          \
00093       }                 \
00094                   \
00095     protected:                \
00096                   \
00097       const ExprT& expr;            \
00098     };                  \
00099                   \
00100     template <typename T>           \
00101     inline Expr< OP< Expr<T> > >          \
00102     OPNAME (const Expr<T>& expr)          \
00103     {                 \
00104       typedef OP< Expr<T> > expr_t;         \
00105                         \
00106       return Expr<expr_t>(expr);          \
00107     }                 \
00108   }                 \
00109 }
00110 
00111 FAD_UNARYOP_MACRO(operator+,
00112       UnaryPlusOp, 
00113       expr.val())
00114 FAD_UNARYOP_MACRO(operator-,
00115       UnaryMinusOp, 
00116       -expr.val())
00117 FAD_UNARYOP_MACRO(exp,
00118       ExpOp, 
00119       std::exp(expr.val()))
00120 FAD_UNARYOP_MACRO(log,
00121       LogOp, 
00122       std::log(expr.val()))
00123 FAD_UNARYOP_MACRO(log10,
00124       Log10Op, 
00125       std::log10(expr.val()))
00126 FAD_UNARYOP_MACRO(sqrt,
00127       SqrtOp, 
00128       std::sqrt(expr.val()))
00129 FAD_UNARYOP_MACRO(cos,
00130       CosOp, 
00131       std::cos(expr.val()))
00132 FAD_UNARYOP_MACRO(sin,
00133       SinOp, 
00134       std::sin(expr.val()))
00135 FAD_UNARYOP_MACRO(tan,
00136       TanOp, 
00137       std::tan(expr.val()))
00138 FAD_UNARYOP_MACRO(acos,
00139       ACosOp, 
00140       std::acos(expr.val()))
00141 FAD_UNARYOP_MACRO(asin,
00142       ASinOp, 
00143       std::asin(expr.val()))
00144 FAD_UNARYOP_MACRO(atan,
00145       ATanOp, 
00146       std::atan(expr.val()))
00147 FAD_UNARYOP_MACRO(cosh,
00148       CoshOp, 
00149       std::cosh(expr.val()))
00150 FAD_UNARYOP_MACRO(sinh,
00151       SinhOp, 
00152       std::sinh(expr.val()))
00153 FAD_UNARYOP_MACRO(tanh,
00154       TanhOp, 
00155       std::tanh(expr.val()))
00156 FAD_UNARYOP_MACRO(acosh,
00157       ACoshOp, 
00158       acosh(expr.val()))
00159 FAD_UNARYOP_MACRO(asinh,
00160       ASinhOp, 
00161       asinh(expr.val()))
00162 FAD_UNARYOP_MACRO(atanh,
00163       ATanhOp, 
00164       atanh(expr.val()))
00165 FAD_UNARYOP_MACRO(abs,
00166       AbsOp, 
00167       std::abs(expr.val()))
00168 FAD_UNARYOP_MACRO(fabs,
00169       FAbsOp, 
00170       std::fabs(expr.val()))
00171 
00172 #undef FAD_UNARYOP_MACRO
00173 
00174 #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) \
00175 namespace Sacado {              \
00176   namespace LFad {              \
00177                   \
00178     template <typename ExprT1, typename ExprT2>       \
00179     class OP {};              \
00180                   \
00181     template <typename T1, typename T2>         \
00182     class Expr< OP< Expr<T1>, Expr<T2> > > {        \
00183                   \
00184     public:               \
00185                   \
00186       typedef Expr<T1> ExprT1;            \
00187       typedef Expr<T2> ExprT2;            \
00188       typedef typename ExprT1::value_type value_type_1;     \
00189       typedef typename ExprT2::value_type value_type_2;     \
00190       typedef typename Sacado::Promote<value_type_1,      \
00191                value_type_2>::type value_type;  \
00192       typedef typename ExprT1::logical_type logical_type;   \
00193                   \
00194       Expr(const ExprT1& expr1_, const ExprT2& expr2_) :    \
00195   expr1(expr1_), expr2(expr2_) {}         \
00196                   \
00197       int size() const {            \
00198   int sz1 = expr1.size(), sz2 = expr2.size();     \
00199   return sz1 > sz2 ? sz1 : sz2;         \
00200       }                 \
00201                   \
00202       bool hasFastAccess() const {          \
00203   return expr1.hasFastAccess() && expr2.hasFastAccess();    \
00204       }                 \
00205                   \
00206       bool isPassive() const {            \
00207   return expr1.isPassive() && expr2.isPassive();      \
00208       }                 \
00209                   \
00210       const value_type val() const {          \
00211   return VALUE;             \
00212       }                 \
00213                   \
00214       const logical_type dx(int i) const {        \
00215   return DX;              \
00216       }                 \
00217                   \
00218       const logical_type fastAccessDx(int i) const {      \
00219   return FASTACCESSDX;            \
00220       }                 \
00221                         \
00222     protected:                \
00223                   \
00224       const ExprT1& expr1;            \
00225       const ExprT2& expr2;            \
00226                   \
00227     };                  \
00228                   \
00229     template <typename T1>            \
00230     class Expr< OP< Expr<T1>, typename Expr<T1>::value_type> > {  \
00231                   \
00232     public:               \
00233                   \
00234       typedef Expr<T1> ExprT1;            \
00235       typedef typename ExprT1::value_type value_type;     \
00236       typedef typename ExprT1::value_type ConstT;     \
00237       typedef typename ExprT1::logical_type logical_type;   \
00238                   \
00239       Expr(const ExprT1& expr1_, const ConstT& c_) :      \
00240   expr1(expr1_), c(c_) {}           \
00241                   \
00242       int size() const {            \
00243   return expr1.size();            \
00244       }                 \
00245                   \
00246       bool hasFastAccess() const {          \
00247   return expr1.hasFastAccess();         \
00248       }                 \
00249                   \
00250       bool isPassive() const {            \
00251   return expr1.isPassive();         \
00252       }                 \
00253                   \
00254       const value_type val() const {          \
00255   return VAL_CONST_DX_2;            \
00256       }                 \
00257                   \
00258       const logical_type dx(int i) const {        \
00259   return CONST_DX_2;            \
00260       }                 \
00261                   \
00262       const logical_type fastAccessDx(int i) const {      \
00263   return CONST_FASTACCESSDX_2;          \
00264       }                 \
00265                   \
00266     protected:                \
00267                   \
00268       const ExprT1& expr1;            \
00269       const ConstT& c;              \
00270     };                  \
00271                   \
00272     template <typename T2>            \
00273     class Expr< OP< typename Expr<T2>::value_type, Expr<T2> > > { \
00274                   \
00275     public:               \
00276                   \
00277       typedef Expr<T2> ExprT2;            \
00278       typedef typename ExprT2::value_type value_type;     \
00279       typedef typename ExprT2::value_type ConstT;     \
00280       typedef typename ExprT2::logical_type logical_type;   \
00281                   \
00282       Expr(const ConstT& c_, const ExprT2& expr2_) :      \
00283   c(c_), expr2(expr2_) {}           \
00284                   \
00285       int size() const {            \
00286   return expr2.size();            \
00287       }                 \
00288                   \
00289       bool hasFastAccess() const {          \
00290   return expr2.hasFastAccess();         \
00291       }                 \
00292                   \
00293       bool isPassive() const {            \
00294   return expr2.isPassive();         \
00295       }                 \
00296                   \
00297       const value_type val() const {          \
00298   return VAL_CONST_DX_1;            \
00299       }                 \
00300                   \
00301       const logical_type dx(int i) const {        \
00302   return CONST_DX_1;            \
00303       }                 \
00304                   \
00305       const logical_type fastAccessDx(int i) const {      \
00306   return CONST_FASTACCESSDX_1;          \
00307       }                 \
00308                         \
00309     protected:                \
00310                   \
00311       const ConstT& c;              \
00312       const ExprT2& expr2;            \
00313     };                  \
00314                   \
00315     template <typename T1, typename T2>         \
00316     inline Expr< OP< Expr<T1>, Expr<T2> > >       \
00317     OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2)   \
00318     {                 \
00319       typedef OP< Expr<T1>, Expr<T2> > expr_t;        \
00320                       \
00321       return Expr<expr_t>(expr1, expr2);        \
00322     }                 \
00323                   \
00324     template <typename T>           \
00325     inline Expr< OP< Expr<T>, Expr<T> > >       \
00326     OPNAME (const Expr<T>& expr1, const Expr<T>& expr2)     \
00327     {                 \
00328       typedef OP< Expr<T>, Expr<T> > expr_t;        \
00329                       \
00330       return Expr<expr_t>(expr1, expr2);        \
00331     }                 \
00332                   \
00333     template <typename T>           \
00334     inline Expr< OP< typename Expr<T>::value_type, Expr<T> > >    \
00335     OPNAME (const typename Expr<T>::value_type& c,      \
00336       const Expr<T>& expr)          \
00337     {                 \
00338       typedef typename Expr<T>::value_type ConstT;      \
00339       typedef OP< ConstT, Expr<T> > expr_t;       \
00340                   \
00341       return Expr<expr_t>(c, expr);         \
00342     }                 \
00343                   \
00344     template <typename T>           \
00345     inline Expr< OP< Expr<T>, typename Expr<T>::value_type > >    \
00346     OPNAME (const Expr<T>& expr,          \
00347       const typename Expr<T>::value_type& c)      \
00348     {                 \
00349       typedef typename Expr<T>::value_type ConstT;      \
00350       typedef OP< Expr<T>, ConstT > expr_t;       \
00351                   \
00352       return Expr<expr_t>(expr, c);         \
00353     }                 \
00354   }                 \
00355 }
00356 
00357 
00358 FAD_BINARYOP_MACRO(operator+,
00359        AdditionOp, 
00360        expr1.val() + expr2.val(),
00361        expr1.dx(i) || expr2.dx(i),
00362        expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
00363        c + expr2.val(),
00364        expr1.val() + c,
00365        expr2.dx(i),
00366        expr1.dx(i),
00367        expr2.fastAccessDx(i),
00368        expr1.fastAccessDx(i))
00369 FAD_BINARYOP_MACRO(operator-,
00370        SubtractionOp, 
00371        expr1.val() - expr2.val(),
00372        expr1.dx(i) || expr2.dx(i),
00373        expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
00374        c - expr2.val(),
00375        expr1.val() - c,
00376        expr2.dx(i),
00377        expr1.dx(i),
00378        expr2.fastAccessDx(i),
00379        expr1.fastAccessDx(i))
00380 FAD_BINARYOP_MACRO(operator*,
00381        MultiplicationOp, 
00382        expr1.val() * expr2.val(),
00383        expr1.dx(i) || expr2.dx(i),
00384        expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
00385        c * expr2.val(),
00386        expr1.val() * c,
00387        expr2.dx(i),
00388        expr1.dx(i),
00389        expr2.fastAccessDx(i),
00390        expr1.fastAccessDx(i))
00391 FAD_BINARYOP_MACRO(operator/,
00392        DivisionOp, 
00393        expr1.val() / expr2.val(),
00394        expr1.dx(i) || expr2.dx(i),
00395        expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
00396        c / expr2.val(),
00397        expr1.val() / c,
00398        expr2.dx(i),
00399        expr1.dx(i),
00400        expr2.fastAccessDx(i),
00401        expr1.fastAccessDx(i))
00402 FAD_BINARYOP_MACRO(atan2,
00403        Atan2Op,
00404        std::atan2(expr1.val(), expr2.val()),
00405        expr1.dx(i) || expr2.dx(i),
00406        expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
00407        std::atan2(c, expr2.val()),
00408        std::atan2(expr1.val(), c),
00409        expr2.dx(i),
00410        expr1.dx(i),
00411        expr2.fastAccessDx(i),
00412        expr1.fastAccessDx(i))
00413 FAD_BINARYOP_MACRO(pow,
00414        PowerOp,
00415        std::pow(expr1.val(), expr2.val()),
00416        expr1.dx(i) || expr2.dx(i),
00417        expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
00418        std::pow(c, expr2.val()),
00419        std::pow(expr1.val(), c),
00420        expr2.dx(i),
00421        expr1.dx(i),
00422        expr2.fastAccessDx(i),
00423        expr1.fastAccessDx(i))
00424 FAD_BINARYOP_MACRO(max,
00425                    MaxOp,
00426                    std::max(expr1.val(), expr2.val()),
00427        expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00428                    expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) : 
00429                                                 expr2.fastAccessDx(i),
00430                    std::max(c, expr2.val()),
00431        std::max(expr1.val(), c),
00432                    c >= expr2.val() ? logical_type(0) : expr2.dx(i),
00433                    expr1.val() >= c ? expr1.dx(i) : logical_type(0),
00434                    c >= expr2.val() ? logical_type(0) : expr2.fastAccessDx(i),
00435                    expr1.val() >= c ? expr1.fastAccessDx(i) : logical_type(0))
00436 FAD_BINARYOP_MACRO(min,
00437                    MinOp,
00438                    std::min(expr1.val(), expr2.val()),
00439        expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
00440                    expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) : 
00441                                                 expr2.fastAccessDx(i),
00442                    std::min(c, expr2.val()),
00443        std::min(expr1.val(), c),
00444        c <= expr2.val() ? logical_type(0) : expr2.dx(i),
00445                    expr1.val() <= c ? expr1.dx(i) : logical_type(0),
00446                    c <= expr2.val() ? logical_type(0) : expr2.fastAccessDx(i),
00447                    expr1.val() <= c ? expr1.fastAccessDx(i) : logical_type(0))
00448 
00449 #undef FAD_BINARYOP_MACRO
00450 
00451 //-------------------------- Relational Operators -----------------------
00452 
00453 #define FAD_RELOP_MACRO(OP)           \
00454 namespace Sacado {              \
00455   namespace LFad {              \
00456     template <typename ExprT1, typename ExprT2>       \
00457     inline bool               \
00458     operator OP (const Expr<ExprT1>& expr1,       \
00459      const Expr<ExprT2>& expr2)       \
00460     {                 \
00461       return expr1.val() OP expr2.val();        \
00462     }                 \
00463                   \
00464     template <typename ExprT2>            \
00465     inline bool               \
00466     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00467      const Expr<ExprT2>& expr2)       \
00468     {                 \
00469       return a OP expr2.val();            \
00470     }                 \
00471                   \
00472     template <typename ExprT1>            \
00473     inline bool               \
00474     operator OP (const Expr<ExprT1>& expr1,       \
00475      const typename Expr<ExprT1>::value_type& b)    \
00476     {                 \
00477       return expr1.val() OP b;            \
00478     }                 \
00479   }                 \
00480 }
00481 
00482 FAD_RELOP_MACRO(==)
00483 FAD_RELOP_MACRO(!=)
00484 FAD_RELOP_MACRO(<)
00485 FAD_RELOP_MACRO(>)
00486 FAD_RELOP_MACRO(<=)
00487 FAD_RELOP_MACRO(>=)
00488 FAD_RELOP_MACRO(<<=)
00489 FAD_RELOP_MACRO(>>=)
00490 FAD_RELOP_MACRO(&)
00491 FAD_RELOP_MACRO(|)
00492 
00493 #undef FAD_RELOP_MACRO
00494 
00495 namespace Sacado {
00496 
00497   namespace LFad {
00498 
00499     template <typename ExprT>
00500     inline bool operator ! (const Expr<ExprT>& expr) 
00501     {
00502       return ! expr.val();
00503     }
00504 
00505   } // namespace LFad
00506 
00507 } // namespace Sacado
00508 
00509 //-------------------------- Boolean Operators -----------------------
00510 namespace Sacado {
00511 
00512   namespace LFad {
00513 
00514     template <typename ExprT>
00515     bool toBool(const Expr<ExprT>& x) {
00516       bool is_zero = (x.val() == 0.0);
00517       for (int i=0; i<x.size(); i++)
00518   is_zero = is_zero && (x.dx(i) == 0.0);
00519       return !is_zero;
00520     }
00521 
00522   } // namespace Fad
00523 
00524 } // namespace Sacado
00525 
00526 #define FAD_BOOL_MACRO(OP)            \
00527 namespace Sacado {              \
00528   namespace LFad {              \
00529     template <typename ExprT1, typename ExprT2>       \
00530     inline bool               \
00531     operator OP (const Expr<ExprT1>& expr1,       \
00532      const Expr<ExprT2>& expr2)       \
00533     {                 \
00534       return toBool(expr1) OP toBool(expr2);        \
00535     }                 \
00536                   \
00537     template <typename ExprT2>            \
00538     inline bool               \
00539     operator OP (const typename Expr<ExprT2>::value_type& a,    \
00540      const Expr<ExprT2>& expr2)       \
00541     {                 \
00542       return a OP toBool(expr2);          \
00543     }                 \
00544                   \
00545     template <typename ExprT1>            \
00546     inline bool               \
00547     operator OP (const Expr<ExprT1>& expr1,       \
00548      const typename Expr<ExprT1>::value_type& b)    \
00549     {                 \
00550       return toBool(expr1) OP b;          \
00551     }                 \
00552   }                 \
00553 }
00554 
00555 FAD_BOOL_MACRO(&&)
00556 FAD_BOOL_MACRO(||)
00557 
00558 #undef FAD_BOOL_MACRO
00559 
00560 //-------------------------- I/O Operators -----------------------
00561 
00562 namespace Sacado {
00563 
00564   namespace LFad {
00565 
00566     template <typename ExprT>
00567     std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
00568       os << x.val() << " [";
00569       
00570       for (int i=0; i< x.size(); i++) {
00571         os << " " << x.dx(i);
00572       }
00573 
00574       os << " ]";
00575       return os;
00576     }
00577 
00578   } // namespace LFad
00579 
00580 } // namespace Sacado
00581 
00582 
00583 #endif // SACADO_LFAD_LOGICALSPARSEOPS_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines