|
Sacado Package Browser (Single Doxygen Collection) Version of the Day
|
|
Sacado Package Browser (Single Doxygen Collection) Version of the Day
|
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
1.7.4
