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