Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_Fad_Ops.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Sacado Package
5// Copyright (2006) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25// (etphipp@sandia.gov).
26//
27// ***********************************************************************
28//
29// The forward-mode AD classes in Sacado are a derivative work of the
30// expression template classes in the Fad package by Nicolas Di Cesare.
31// The following banner is included in the original Fad source code:
32//
33// ************ DO NOT REMOVE THIS BANNER ****************
34//
35// Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
36// http://www.ann.jussieu.fr/~dicesare
37//
38// CEMRACS 98 : C++ courses,
39// templates : new C++ techniques
40// for scientific computing
41//
42//********************************************************
43//
44// A short implementation ( not all operators and
45// functions are overloaded ) of 1st order Automatic
46// Differentiation in forward mode (FAD) using
47// EXPRESSION TEMPLATES.
48//
49//********************************************************
50// @HEADER
51
52#ifndef SACADO_FAD_OPS_HPP
53#define SACADO_FAD_OPS_HPP
54
57#include "Sacado_cmath.hpp"
58#include <ostream> // for std::ostream
59
60#define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
61namespace Sacado { \
62 namespace Fad { \
63 \
64 template <typename ExprT> \
65 class OP {}; \
66 \
67 template <typename ExprT> \
68 struct ExprSpec< OP<ExprT> > { \
69 typedef typename ExprSpec<ExprT>::type type; \
70 }; \
71 \
72 template <typename ExprT> \
73 class Expr< OP<ExprT>,ExprSpecDefault > { \
74 public: \
75 \
76 typedef typename ExprT::value_type value_type; \
77 typedef typename ExprT::scalar_type scalar_type; \
78 typedef typename ExprT::base_expr_type base_expr_type; \
79 \
80 SACADO_INLINE_FUNCTION \
81 explicit Expr(const ExprT& expr_) : expr(expr_) {} \
82 \
83 SACADO_INLINE_FUNCTION \
84 int size() const { return expr.size(); } \
85 \
86 SACADO_INLINE_FUNCTION \
87 bool hasFastAccess() const { return expr.hasFastAccess(); } \
88 \
89 SACADO_INLINE_FUNCTION \
90 bool isPassive() const { return expr.isPassive();} \
91 \
92 SACADO_INLINE_FUNCTION \
93 bool updateValue() const { return expr.updateValue(); } \
94 \
95 SACADO_INLINE_FUNCTION \
96 void cache() const {} \
97 \
98 SACADO_INLINE_FUNCTION \
99 value_type val() const { \
100 USING \
101 return VALUE; \
102 } \
103 \
104 SACADO_INLINE_FUNCTION \
105 value_type dx(int i) const { \
106 USING \
107 return DX; \
108 } \
109 \
110 SACADO_INLINE_FUNCTION \
111 value_type fastAccessDx(int i) const { \
112 USING \
113 return FASTACCESSDX; \
114 } \
115 \
116 protected: \
117 \
118 const ExprT& expr; \
119 }; \
120 \
121 template <typename T> \
122 SACADO_INLINE_FUNCTION \
123 Expr< OP< Expr<T> > > \
124 OPNAME (const Expr<T>& expr) \
125 { \
126 typedef OP< Expr<T> > expr_t; \
127 \
128 return Expr<expr_t>(expr); \
129 } \
130 } \
131 \
132}
133
135 UnaryPlusOp,
136 ;,
137 expr.val(),
138 expr.dx(i),
139 expr.fastAccessDx(i))
140FAD_UNARYOP_MACRO(operator-,
142 ;,
143 -expr.val(),
144 -expr.dx(i),
145 -expr.fastAccessDx(i))
148 using std::exp;,
149 exp(expr.val()),
150 exp(expr.val())*expr.dx(i),
151 exp(expr.val())*expr.fastAccessDx(i))
154 using std::log;,
155 log(expr.val()),
156 expr.dx(i)/expr.val(),
157 expr.fastAccessDx(i)/expr.val())
160 using std::log10; using std::log;,
161 log10(expr.val()),
162 expr.dx(i)/( log(value_type(10))*expr.val()),
163 expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
166 using std::sqrt;,
167 sqrt(expr.val()),
168 expr.dx(i)/(value_type(2)* sqrt(expr.val())),
169 expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
172 using std::cos; using std::sin;,
173 cos(expr.val()),
174 -expr.dx(i)* sin(expr.val()),
175 -expr.fastAccessDx(i)* sin(expr.val()))
178 using std::cos; using std::sin;,
179 sin(expr.val()),
180 expr.dx(i)* cos(expr.val()),
181 expr.fastAccessDx(i)* cos(expr.val()))
184 using std::tan;,
185 std::tan(expr.val()),
186 expr.dx(i)*
187 (value_type(1)+ tan(expr.val())* tan(expr.val())),
188 expr.fastAccessDx(i)*
189 (value_type(1)+ tan(expr.val())* tan(expr.val())))
192 using std::acos; using std::sqrt;,
193 acos(expr.val()),
194 -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
195 -expr.fastAccessDx(i) /
196 sqrt(value_type(1)-expr.val()*expr.val()))
199 using std::asin; using std::sqrt;,
200 asin(expr.val()),
201 expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
202 expr.fastAccessDx(i) /
203 sqrt(value_type(1)-expr.val()*expr.val()))
206 using std::atan;,
207 atan(expr.val()),
208 expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
209 expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
212 using std::cosh; using std::sinh;,
213 cosh(expr.val()),
214 expr.dx(i)* sinh(expr.val()),
215 expr.fastAccessDx(i)* sinh(expr.val()))
218 using std::cosh; using std::sinh;,
219 sinh(expr.val()),
220 expr.dx(i)* cosh(expr.val()),
221 expr.fastAccessDx(i)* cosh(expr.val()))
224 using std::tanh;,
225 tanh(expr.val()),
226 expr.dx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())),
227 expr.fastAccessDx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())))
230 using std::acosh; using std::sqrt;,
231 acosh(expr.val()),
232 expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
233 (expr.val()+value_type(1))),
234 expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
235 (expr.val()+value_type(1))))
238 using std::asinh; using std::sqrt;,
239 asinh(expr.val()),
240 expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
241 expr.fastAccessDx(i)/ sqrt(value_type(1)+
242 expr.val()*expr.val()))
245 using std::atanh;,
246 atanh(expr.val()),
247 expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
248 expr.fastAccessDx(i)/(value_type(1)-
249 expr.val()*expr.val()))
252 using std::abs; using Sacado::if_then_else;,
253 abs(expr.val()),
254 if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
255 if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
258 using std::fabs; using Sacado::if_then_else;,
259 fabs(expr.val()),
260 if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
261 if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
264 using std::cbrt;,
265 cbrt(expr.val()),
266 expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
267 expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
268
269#undef FAD_UNARYOP_MACRO
270
271// Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
272// "simd" value types that use if_then_else(). The only reason for not using
273// if_then_else() always is to avoid evaluating the derivative if the value is
274// zero to avoid throwing FPEs.
275namespace Sacado {
276 namespace Fad {
277
278 template <typename ExprT, bool is_simd>
279 class SafeSqrtOp {};
280
281 template <typename ExprT>
282 struct ExprSpec< SafeSqrtOp<ExprT> > {
283 typedef typename ExprSpec<ExprT>::type type;
284 };
285
286 //
287 // Implementation for simd type using if_then_else()
288 //
289 template <typename ExprT>
290 class Expr< SafeSqrtOp<ExprT,true>,ExprSpecDefault > {
291 public:
292
293 typedef typename ExprT::value_type value_type;
294 typedef typename ExprT::scalar_type scalar_type;
295 typedef typename ExprT::base_expr_type base_expr_type;
296
298 explicit Expr(const ExprT& expr_) : expr(expr_) {}
299
301 int size() const { return expr.size(); }
302
304 bool hasFastAccess() const { return expr.hasFastAccess(); }
305
307 bool isPassive() const { return expr.isPassive();}
308
310 bool updateValue() const { return expr.updateValue(); }
311
313 void cache() const {}
314
316 value_type val() const {
317 using std::sqrt;
318 return sqrt(expr.val());
319 }
320
322 value_type dx(int i) const {
323 using std::sqrt; using Sacado::if_then_else;
324 return if_then_else(
325 expr.val() == value_type(0.0), value_type(0.0),
326 value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
327 }
328
330 value_type fastAccessDx(int i) const {
331 using std::sqrt; using Sacado::if_then_else;
332 return if_then_else(
333 expr.val() == value_type(0.0), value_type(0.0),
334 value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
335 }
336
337 protected:
338
339 const ExprT& expr;
340 };
341
342 //
343 // Specialization for scalar types using ternary operator
344 //
345 template <typename ExprT>
346 class Expr< SafeSqrtOp<ExprT,false>,ExprSpecDefault > {
347 public:
348
349 typedef typename ExprT::value_type value_type;
350 typedef typename ExprT::scalar_type scalar_type;
351 typedef typename ExprT::base_expr_type base_expr_type;
352
354 explicit Expr(const ExprT& expr_) : expr(expr_) {}
355
357 int size() const { return expr.size(); }
358
360 bool hasFastAccess() const { return expr.hasFastAccess(); }
361
363 bool isPassive() const { return expr.isPassive();}
364
366 bool updateValue() const { return expr.updateValue(); }
367
369 void cache() const {}
370
372 value_type val() const {
373 using std::sqrt;
374 return sqrt(expr.val());
375 }
376
378 value_type dx(int i) const {
379 using std::sqrt;
380 return expr.val() == value_type(0.0) ? value_type(0.0) :
381 value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
382 }
383
385 value_type fastAccessDx(int i) const {
386 using std::sqrt;
387 return expr.val() == value_type(0.0) ? value_type(0.0) :
388 value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
389 }
390
391 protected:
392
393 const ExprT& expr;
394 };
395
396 template <typename T>
398 Expr< SafeSqrtOp< Expr<T> > >
399 safe_sqrt (const Expr<T>& expr)
400 {
401 typedef SafeSqrtOp< Expr<T> > expr_t;
402
403 return Expr<expr_t>(expr);
404 }
405 }
406
407}
408
409#define FAD_BINARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX,VAL_CONST_DX_1,VAL_CONST_DX_2,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
410namespace Sacado { \
411 namespace Fad { \
412 \
413 template <typename ExprT1, typename ExprT2> \
414 class OP {}; \
415 \
416 template <typename ExprT1, typename ExprT2> \
417 struct ExprSpec< OP< ExprT1, ExprT2 > > { \
418 typedef typename ExprSpec<ExprT1>::type type; \
419 }; \
420 \
421 template <typename ExprT1, typename ExprT2> \
422 class Expr< OP< ExprT1, ExprT2 >,ExprSpecDefault > { \
423 \
424 public: \
425 \
426 typedef typename ExprT1::value_type value_type_1; \
427 typedef typename ExprT2::value_type value_type_2; \
428 typedef typename Sacado::Promote<value_type_1, \
429 value_type_2>::type value_type; \
430 \
431 typedef typename ExprT1::scalar_type scalar_type_1; \
432 typedef typename ExprT2::scalar_type scalar_type_2; \
433 typedef typename Sacado::Promote<scalar_type_1, \
434 scalar_type_2>::type scalar_type; \
435 \
436 typedef typename ExprT1::base_expr_type base_expr_type_1; \
437 typedef typename ExprT2::base_expr_type base_expr_type_2; \
438 typedef typename Sacado::Promote<base_expr_type_1, \
439 base_expr_type_2>::type base_expr_type; \
440 \
441 SACADO_INLINE_FUNCTION \
442 Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
443 expr1(expr1_), expr2(expr2_) {} \
444 \
445 SACADO_INLINE_FUNCTION \
446 int size() const { \
447 int sz1 = expr1.size(), sz2 = expr2.size(); \
448 return sz1 > sz2 ? sz1 : sz2; \
449 } \
450 \
451 SACADO_INLINE_FUNCTION \
452 bool hasFastAccess() const { \
453 return expr1.hasFastAccess() && expr2.hasFastAccess(); \
454 } \
455 \
456 SACADO_INLINE_FUNCTION \
457 bool isPassive() const { \
458 return expr1.isPassive() && expr2.isPassive(); \
459 } \
460 \
461 SACADO_INLINE_FUNCTION \
462 bool updateValue() const { \
463 return expr1.updateValue() && expr2.updateValue(); \
464 } \
465 \
466 SACADO_INLINE_FUNCTION \
467 void cache() const {} \
468 \
469 SACADO_INLINE_FUNCTION \
470 const value_type val() const { \
471 USING \
472 return VALUE; \
473 } \
474 \
475 SACADO_INLINE_FUNCTION \
476 const value_type dx(int i) const { \
477 USING \
478 return DX; \
479 } \
480 \
481 SACADO_INLINE_FUNCTION \
482 const value_type fastAccessDx(int i) const { \
483 USING \
484 return FASTACCESSDX; \
485 } \
486 \
487 protected: \
488 \
489 const ExprT1& expr1; \
490 const ExprT2& expr2; \
491 \
492 }; \
493 \
494 template <typename ExprT1, typename T2> \
495 struct ExprSpec< OP< ExprT1, ConstExpr<T2> > > { \
496 typedef typename ExprSpec<ExprT1>::type type; \
497 }; \
498 \
499 template <typename ExprT1, typename T2> \
500 class Expr< OP< ExprT1, ConstExpr<T2> >,ExprSpecDefault > { \
501 \
502 public: \
503 \
504 typedef ConstExpr<T2> ConstT; \
505 typedef ConstExpr<T2> ExprT2; \
506 typedef typename ExprT1::value_type value_type_1; \
507 typedef typename ExprT2::value_type value_type_2; \
508 typedef typename Sacado::Promote<value_type_1, \
509 value_type_2>::type value_type; \
510 \
511 typedef typename ExprT1::scalar_type scalar_type_1; \
512 typedef typename ExprT2::scalar_type scalar_type_2; \
513 typedef typename Sacado::Promote<scalar_type_1, \
514 scalar_type_2>::type scalar_type; \
515 \
516 typedef typename ExprT1::base_expr_type base_expr_type_1; \
517 typedef typename ExprT2::base_expr_type base_expr_type_2; \
518 typedef typename Sacado::Promote<base_expr_type_1, \
519 base_expr_type_2>::type base_expr_type; \
520 \
521 SACADO_INLINE_FUNCTION \
522 Expr(const ExprT1& expr1_, const ConstT& c_) : \
523 expr1(expr1_), c(c_) {} \
524 \
525 SACADO_INLINE_FUNCTION \
526 int size() const { \
527 return expr1.size(); \
528 } \
529 \
530 SACADO_INLINE_FUNCTION \
531 bool hasFastAccess() const { \
532 return expr1.hasFastAccess(); \
533 } \
534 \
535 SACADO_INLINE_FUNCTION \
536 bool isPassive() const { \
537 return expr1.isPassive(); \
538 } \
539 \
540 SACADO_INLINE_FUNCTION \
541 bool updateValue() const { return expr1.updateValue(); } \
542 \
543 SACADO_INLINE_FUNCTION \
544 void cache() const {} \
545 \
546 SACADO_INLINE_FUNCTION \
547 const value_type val() const { \
548 USING \
549 return VAL_CONST_DX_2; \
550 } \
551 \
552 SACADO_INLINE_FUNCTION \
553 const value_type dx(int i) const { \
554 USING \
555 return CONST_DX_2; \
556 } \
557 \
558 SACADO_INLINE_FUNCTION \
559 const value_type fastAccessDx(int i) const { \
560 USING \
561 return CONST_FASTACCESSDX_2; \
562 } \
563 \
564 protected: \
565 \
566 const ExprT1& expr1; \
567 ConstT c; \
568 }; \
569 \
570 template <typename T1, typename ExprT2> \
571 struct ExprSpec< OP< ConstExpr<T1>, ExprT2 > > { \
572 typedef typename ExprSpec<ExprT2>::type type; \
573 }; \
574 \
575 template <typename T1, typename ExprT2> \
576 class Expr< OP< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > { \
577 \
578 public: \
579 \
580 typedef ConstExpr<T1> ConstT; \
581 typedef ConstExpr<T1> ExprT1; \
582 typedef typename ExprT1::value_type value_type_1; \
583 typedef typename ExprT2::value_type value_type_2; \
584 typedef typename Sacado::Promote<value_type_1, \
585 value_type_2>::type value_type; \
586 \
587 typedef typename ExprT1::scalar_type scalar_type_1; \
588 typedef typename ExprT2::scalar_type scalar_type_2; \
589 typedef typename Sacado::Promote<scalar_type_1, \
590 scalar_type_2>::type scalar_type; \
591 \
592 typedef typename ExprT1::base_expr_type base_expr_type_1; \
593 typedef typename ExprT2::base_expr_type base_expr_type_2; \
594 typedef typename Sacado::Promote<base_expr_type_1, \
595 base_expr_type_2>::type base_expr_type; \
596 \
597 \
598 SACADO_INLINE_FUNCTION \
599 Expr(const ConstT& c_, const ExprT2& expr2_) : \
600 c(c_), expr2(expr2_) {} \
601 \
602 SACADO_INLINE_FUNCTION \
603 int size() const { \
604 return expr2.size(); \
605 } \
606 \
607 SACADO_INLINE_FUNCTION \
608 bool hasFastAccess() const { \
609 return expr2.hasFastAccess(); \
610 } \
611 \
612 SACADO_INLINE_FUNCTION \
613 bool isPassive() const { \
614 return expr2.isPassive(); \
615 } \
616 \
617 SACADO_INLINE_FUNCTION \
618 bool updateValue() const { return expr2.updateValue(); } \
619 \
620 SACADO_INLINE_FUNCTION \
621 void cache() const {} \
622 \
623 SACADO_INLINE_FUNCTION \
624 const value_type val() const { \
625 USING \
626 return VAL_CONST_DX_1; \
627 } \
628 \
629 SACADO_INLINE_FUNCTION \
630 const value_type dx(int i) const { \
631 USING \
632 return CONST_DX_1; \
633 } \
634 \
635 SACADO_INLINE_FUNCTION \
636 const value_type fastAccessDx(int i) const { \
637 USING \
638 return CONST_FASTACCESSDX_1; \
639 } \
640 \
641 protected: \
642 \
643 ConstT c; \
644 const ExprT2& expr2; \
645 }; \
646 \
647 template <typename T1, typename T2> \
648 SACADO_INLINE_FUNCTION \
649 typename mpl::enable_if_c< \
650 ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value, \
651 Expr< OP< Expr<T1>, Expr<T2> > > \
652 >::type \
653 /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(OP)*/ \
654 OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2) \
655 { \
656 typedef OP< Expr<T1>, Expr<T2> > expr_t; \
657 \
658 return Expr<expr_t>(expr1, expr2); \
659 } \
660 \
661 template <typename T> \
662 SACADO_INLINE_FUNCTION \
663 Expr< OP< Expr<T>, Expr<T> > > \
664 OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
665 { \
666 typedef OP< Expr<T>, Expr<T> > expr_t; \
667 \
668 return Expr<expr_t>(expr1, expr2); \
669 } \
670 \
671 template <typename T> \
672 SACADO_INLINE_FUNCTION \
673 Expr< OP< ConstExpr<typename Expr<T>::value_type>, \
674 Expr<T> > > \
675 OPNAME (const typename Expr<T>::value_type& c, \
676 const Expr<T>& expr) \
677 { \
678 typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
679 typedef OP< ConstT, Expr<T> > expr_t; \
680 \
681 return Expr<expr_t>(ConstT(c), expr); \
682 } \
683 \
684 template <typename T> \
685 SACADO_INLINE_FUNCTION \
686 Expr< OP< Expr<T>, \
687 ConstExpr<typename Expr<T>::value_type> > > \
688 OPNAME (const Expr<T>& expr, \
689 const typename Expr<T>::value_type& c) \
690 { \
691 typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
692 typedef OP< Expr<T>, ConstT > expr_t; \
693 \
694 return Expr<expr_t>(expr, ConstT(c)); \
695 } \
696 \
697 template <typename T> \
698 SACADO_INLINE_FUNCTION \
699 SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP) \
700 OPNAME (const typename Expr<T>::scalar_type& c, \
701 const Expr<T>& expr) \
702 { \
703 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
704 typedef OP< ConstT, Expr<T> > expr_t; \
705 \
706 return Expr<expr_t>(ConstT(c), expr); \
707 } \
708 \
709 template <typename T> \
710 SACADO_INLINE_FUNCTION \
711 SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP) \
712 OPNAME (const Expr<T>& expr, \
713 const typename Expr<T>::scalar_type& c) \
714 { \
715 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
716 typedef OP< Expr<T>, ConstT > expr_t; \
717 \
718 return Expr<expr_t>(expr, ConstT(c)); \
719 } \
720 } \
721 \
722}
723
724
727 ;,
728 expr1.val() + expr2.val(),
729 expr1.dx(i) + expr2.dx(i),
730 expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
731 c.val() + expr2.val(),
732 expr1.val() + c.val(),
733 expr2.dx(i),
734 expr1.dx(i),
735 expr2.fastAccessDx(i),
736 expr1.fastAccessDx(i))
737FAD_BINARYOP_MACRO(operator-,
739 ;,
740 expr1.val() - expr2.val(),
741 expr1.dx(i) - expr2.dx(i),
742 expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
743 c.val() - expr2.val(),
744 expr1.val() - c.val(),
745 -expr2.dx(i),
746 expr1.dx(i),
747 -expr2.fastAccessDx(i),
748 expr1.fastAccessDx(i))
749// FAD_BINARYOP_MACRO(operator*,
750// MultiplicationOp,
751// ;,
752// expr1.val() * expr2.val(),
753// expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
754// expr1.val()*expr2.fastAccessDx(i) +
755// expr1.fastAccessDx(i)*expr2.val(),
756// c.val() * expr2.val(),
757// expr1.val() * c.val(),
758// c.val()*expr2.dx(i),
759// expr1.dx(i)*c.val(),
760// c.val()*expr2.fastAccessDx(i),
761// expr1.fastAccessDx(i)*c.val())
762FAD_BINARYOP_MACRO(operator/,
764 ;,
765 expr1.val() / expr2.val(),
766 (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
767 (expr2.val()*expr2.val()),
768 (expr1.fastAccessDx(i)*expr2.val() -
769 expr2.fastAccessDx(i)*expr1.val()) /
770 (expr2.val()*expr2.val()),
771 c.val() / expr2.val(),
772 expr1.val() / c.val(),
773 -expr2.dx(i)*c.val() / (expr2.val()*expr2.val()),
774 expr1.dx(i)/c.val(),
775 -expr2.fastAccessDx(i)*c.val() / (expr2.val()*expr2.val()),
776 expr1.fastAccessDx(i)/c.val())
779 using std::atan2;,
780 atan2(expr1.val(), expr2.val()),
781 (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
782 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
783 (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
784 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
785 atan2(c.val(), expr2.val()),
786 atan2(expr1.val(), c.val()),
787 (-c.val()*expr2.dx(i)) / (c.val()*c.val() + expr2.val()*expr2.val()),
788 (c.val()*expr1.dx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()),
789 (-c.val()*expr2.fastAccessDx(i))/ (c.val()*c.val() + expr2.val()*expr2.val()),
790 (c.val()*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()))
791// FAD_BINARYOP_MACRO(pow,
792// PowerOp,
793// using std::pow; using std::log; using Sacado::if_then_else;,
794// pow(expr1.val(), expr2.val()),
795// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val())) ),
796// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type((expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val())) ),
797// pow(c.val(), expr2.val()),
798// pow(expr1.val(), c.val()),
799// if_then_else( c.val() == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val())) ),
800// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val())) ),
801// if_then_else( c.val() == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val())) ),
802// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()))) )
805 using Sacado::if_then_else;,
806 if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
807 if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
808 if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
809 if_then_else( c.val() >= expr2.val(), value_type(c.val()), expr2.val() ),
810 if_then_else( expr1.val() >= c.val(), expr1.val(), value_type(c.val()) ),
811 if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
812 if_then_else( expr1.val() >= c.val(), expr1.dx(i), value_type(0.0) ),
813 if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
814 if_then_else( expr1.val() >= c.val(), expr1.fastAccessDx(i), value_type(0.0) ) )
817 using Sacado::if_then_else;,
818 if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
819 if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
820 if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
821 if_then_else( c.val() <= expr2.val(), value_type(c.val()), expr2.val() ),
822 if_then_else( expr1.val() <= c.val(), expr1.val(), value_type(c.val()) ),
823 if_then_else( c.val() <= expr2.val(), value_type(0), expr2.dx(i) ),
824 if_then_else( expr1.val() <= c.val(), expr1.dx(i), value_type(0) ),
825 if_then_else( c.val() <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
826 if_then_else( expr1.val() <= c.val(), expr1.fastAccessDx(i), value_type(0) ) )
827
828
829#undef FAD_BINARYOP_MACRO
830
831namespace Sacado {
832 namespace Fad {
833
834 template <typename ExprT1, typename ExprT2>
835 class MultiplicationOp {};
836
837 template <typename ExprT1, typename ExprT2>
838 struct ExprSpec< MultiplicationOp< ExprT1, ExprT2 > > {
839 typedef typename ExprSpec<ExprT1>::type type;
840 };
841
842 template <typename ExprT1, typename ExprT2>
843 class Expr< MultiplicationOp< ExprT1, ExprT2 >,ExprSpecDefault > {
844
845 public:
846
847 typedef typename ExprT1::value_type value_type_1;
848 typedef typename ExprT2::value_type value_type_2;
849 typedef typename Sacado::Promote<value_type_1,
850 value_type_2>::type value_type;
851
852 typedef typename ExprT1::scalar_type scalar_type_1;
853 typedef typename ExprT2::scalar_type scalar_type_2;
854 typedef typename Sacado::Promote<scalar_type_1,
855 scalar_type_2>::type scalar_type;
856
857 typedef typename ExprT1::base_expr_type base_expr_type_1;
858 typedef typename ExprT2::base_expr_type base_expr_type_2;
859 typedef typename Sacado::Promote<base_expr_type_1,
860 base_expr_type_2>::type base_expr_type;
861
863 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
864 expr1(expr1_), expr2(expr2_) {}
865
867 int size() const {
868 int sz1 = expr1.size(), sz2 = expr2.size();
869 return sz1 > sz2 ? sz1 : sz2;
870 }
871
873 bool hasFastAccess() const {
874 return expr1.hasFastAccess() && expr2.hasFastAccess();
875 }
876
878 bool isPassive() const {
879 return expr1.isPassive() && expr2.isPassive();
880 }
881
883 bool updateValue() const {
884 return expr1.updateValue() && expr2.updateValue();
885 }
886
888 void cache() const {}
889
891 const value_type val() const {
892 return expr1.val()*expr2.val();
893 }
894
896 const value_type dx(int i) const {
897 if (expr1.size() > 0 && expr2.size() > 0)
898 return expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val();
899 else if (expr1.size() > 0)
900 return expr1.dx(i)*expr2.val();
901 else
902 return expr1.val()*expr2.dx(i);
903 }
904
906 const value_type fastAccessDx(int i) const {
907 return expr1.val()*expr2.fastAccessDx(i) +
908 expr1.fastAccessDx(i)*expr2.val();
909 }
910
911 protected:
912
913 const ExprT1& expr1;
914 const ExprT2& expr2;
915
916 };
917
918 template <typename ExprT1, typename T2>
919 struct ExprSpec< MultiplicationOp< ExprT1, ConstExpr<T2> > > {
920 typedef typename ExprSpec<ExprT1>::type type;
921 };
922
923 template <typename ExprT1, typename T2>
924 class Expr< MultiplicationOp< ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
925
926 public:
927
928 typedef ConstExpr<T2> ConstT;
929 typedef ConstExpr<T2> ExprT2;
930 typedef typename ExprT1::value_type value_type_1;
931 typedef typename ExprT2::value_type value_type_2;
932 typedef typename Sacado::Promote<value_type_1,
933 value_type_2>::type value_type;
934
935 typedef typename ExprT1::scalar_type scalar_type_1;
936 typedef typename ExprT2::scalar_type scalar_type_2;
937 typedef typename Sacado::Promote<scalar_type_1,
938 scalar_type_2>::type scalar_type;
939
940 typedef typename ExprT1::base_expr_type base_expr_type_1;
941 typedef typename ExprT2::base_expr_type base_expr_type_2;
942 typedef typename Sacado::Promote<base_expr_type_1,
943 base_expr_type_2>::type base_expr_type;
944
946 Expr(const ExprT1& expr1_, const ConstT& c_) :
947 expr1(expr1_), c(c_) {}
948
950 int size() const {
951 return expr1.size();
952 }
953
955 bool hasFastAccess() const {
956 return expr1.hasFastAccess();
957 }
958
960 bool isPassive() const {
961 return expr1.isPassive();
962 }
963
965 bool updateValue() const { return expr1.updateValue(); }
966
968 void cache() const {}
969
971 const value_type val() const {
972 return expr1.val()*c.val();
973 }
974
976 const value_type dx(int i) const {
977 return expr1.dx(i)*c.val();
978 }
979
981 const value_type fastAccessDx(int i) const {
982 return expr1.fastAccessDx(i)*c.val();
983 }
984
985 protected:
986
987 const ExprT1& expr1;
988 ConstT c;
989 };
990
991 template <typename T1, typename ExprT2>
992 struct ExprSpec< MultiplicationOp< ConstExpr<T1>, ExprT2 > > {
993 typedef typename ExprSpec<ExprT2>::type type;
994 };
995
996 template <typename T1, typename ExprT2>
997 class Expr< MultiplicationOp< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
998
999 public:
1000
1001 typedef ConstExpr<T1> ConstT;
1002 typedef ConstExpr<T1> ExprT1;
1003 typedef typename ExprT1::value_type value_type_1;
1004 typedef typename ExprT2::value_type value_type_2;
1005 typedef typename Sacado::Promote<value_type_1,
1006 value_type_2>::type value_type;
1007
1008 typedef typename ExprT1::scalar_type scalar_type_1;
1009 typedef typename ExprT2::scalar_type scalar_type_2;
1010 typedef typename Sacado::Promote<scalar_type_1,
1011 scalar_type_2>::type scalar_type;
1012
1013 typedef typename ExprT1::base_expr_type base_expr_type_1;
1014 typedef typename ExprT2::base_expr_type base_expr_type_2;
1015 typedef typename Sacado::Promote<base_expr_type_1,
1016 base_expr_type_2>::type base_expr_type;
1017
1019 Expr(const ConstT& c_, const ExprT2& expr2_) :
1020 c(c_), expr2(expr2_) {}
1021
1023 int size() const {
1024 return expr2.size();
1025 }
1026
1028 bool hasFastAccess() const {
1029 return expr2.hasFastAccess();
1030 }
1031
1033 bool isPassive() const {
1034 return expr2.isPassive();
1035 }
1036
1038 bool updateValue() const { return expr2.updateValue(); }
1039
1041 void cache() const {}
1042
1044 const value_type val() const {
1045 return c.val()*expr2.val();
1046 }
1047
1049 const value_type dx(int i) const {
1050 return c.val()*expr2.dx(i);
1051 }
1052
1054 const value_type fastAccessDx(int i) const {
1055 return c.val()*expr2.fastAccessDx(i);
1056 }
1057
1058 protected:
1059
1060 ConstT c;
1061 const ExprT2& expr2;
1062 };
1063
1064 template <typename T1, typename T2>
1066 typename mpl::enable_if_c<
1067 ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value,
1068 Expr< MultiplicationOp< Expr<T1>, Expr<T2> > >
1069 >::type
1070 /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(MultiplicationOp)*/
1071 operator* (const Expr<T1>& expr1, const Expr<T2>& expr2)
1072 {
1073 typedef MultiplicationOp< Expr<T1>, Expr<T2> > expr_t;
1074
1075 return Expr<expr_t>(expr1, expr2);
1076 }
1077
1078 template <typename T>
1080 Expr< MultiplicationOp< Expr<T>, Expr<T> > >
1081 operator* (const Expr<T>& expr1, const Expr<T>& expr2)
1082 {
1083 typedef MultiplicationOp< Expr<T>, Expr<T> > expr_t;
1084
1085 return Expr<expr_t>(expr1, expr2);
1086 }
1087
1088 template <typename T>
1090 Expr< MultiplicationOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
1091 operator* (const typename Expr<T>::value_type& c,
1092 const Expr<T>& expr)
1093 {
1094 typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1095 typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1096
1097 return Expr<expr_t>(ConstT(c), expr);
1098 }
1099
1100 template <typename T>
1102 Expr< MultiplicationOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
1103 operator* (const Expr<T>& expr,
1104 const typename Expr<T>::value_type& c)
1105 {
1106 typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1107 typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1108
1109 return Expr<expr_t>(expr, ConstT(c));
1110 }
1111
1112 template <typename T>
1115 operator* (const typename Expr<T>::scalar_type& c,
1116 const Expr<T>& expr)
1117 {
1118 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1119 typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1120
1121 return Expr<expr_t>(ConstT(c), expr);
1122 }
1123
1124 template <typename T>
1127 operator* (const Expr<T>& expr,
1128 const typename Expr<T>::scalar_type& c)
1129 {
1130 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1131 typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1132
1133 return Expr<expr_t>(expr, ConstT(c));
1134 }
1135 }
1136}
1137
1138// Special handling for std::pow() to provide specializations of PowerOp for
1139// "simd" value types that use if_then_else(). The only reason for not using
1140// if_then_else() always is to avoid evaluating the derivative if the value is
1141// zero to avoid throwing FPEs.
1142namespace Sacado {
1143 namespace Fad {
1144
1145 template <typename ExprT1, typename ExprT2, typename Impl>
1146 class PowerOp {};
1147
1148 template <typename ExprT1, typename ExprT2>
1149 struct ExprSpec< PowerOp< ExprT1, ExprT2 > > {
1151 };
1152
1153 template <typename ExprT1, typename T2>
1154 struct ExprSpec<PowerOp< ExprT1, ConstExpr<T2> > > {
1156 };
1157
1158 template <typename T1, typename ExprT2>
1159 struct ExprSpec< PowerOp< ConstExpr<T1>, ExprT2 > > {
1161 };
1162
1163 //
1164 // Implementation for simd type using if_then_else()
1165 //
1166 template <typename ExprT1, typename ExprT2>
1167 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Simd >, ExprSpecDefault > {
1168
1169 public:
1170
1171 typedef typename ExprT1::value_type value_type_1;
1172 typedef typename ExprT2::value_type value_type_2;
1173 typedef typename Sacado::Promote<value_type_1,
1175
1176 typedef typename ExprT1::scalar_type scalar_type_1;
1177 typedef typename ExprT2::scalar_type scalar_type_2;
1178 typedef typename Sacado::Promote<scalar_type_1,
1180
1181 typedef typename ExprT1::base_expr_type base_expr_type_1;
1182 typedef typename ExprT2::base_expr_type base_expr_type_2;
1183 typedef typename Sacado::Promote<base_expr_type_1,
1185
1187 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1188 expr1(expr1_), expr2(expr2_) {}
1189
1191 int size() const {
1192 int sz1 = expr1.size(), sz2 = expr2.size();
1193 return sz1 > sz2 ? sz1 : sz2;
1194 }
1195
1197 bool hasFastAccess() const {
1198 return expr1.hasFastAccess() && expr2.hasFastAccess();
1199 }
1200
1202 bool isPassive() const {
1203 return expr1.isPassive() && expr2.isPassive();
1204 }
1205
1207 bool updateValue() const {
1208 return expr1.updateValue() && expr2.updateValue();
1209 }
1210
1212 void cache() const {}
1213
1215 const value_type val() const {
1216 using std::pow;
1217 return pow(expr1.val(), expr2.val());
1218 }
1219
1221 const value_type dx(int i) const {
1222 using std::pow; using std::log; using Sacado::if_then_else;
1223 const int sz1 = expr1.size(), sz2 = expr2.size();
1224 if (sz1 > 0 && sz2 > 0)
1225 return if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val())) );
1226 else if (sz1 > 0)
1227 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1228 // It seems less accurate and caused convergence problems in some codes
1229 return if_then_else( expr2.val() == scalar_type(1.0), expr1.dx(i), if_then_else(expr1.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val())) ));
1230 else
1231 return if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val())) );
1232 }
1233
1235 const value_type fastAccessDx(int i) const {
1236 using std::pow; using std::log; using Sacado::if_then_else;
1237 return if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type((expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val())) );
1238 }
1239
1240 protected:
1241
1242 const ExprT1& expr1;
1243 const ExprT2& expr2;
1244
1245 };
1246
1247 template <typename ExprT1, typename T2>
1248 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Simd >,
1249 ExprSpecDefault > {
1250
1251 public:
1252
1255 typedef typename ExprT1::value_type value_type_1;
1257 typedef typename Sacado::Promote<value_type_1,
1259
1260 typedef typename ExprT1::scalar_type scalar_type_1;
1262 typedef typename Sacado::Promote<scalar_type_1,
1264
1265 typedef typename ExprT1::base_expr_type base_expr_type_1;
1267 typedef typename Sacado::Promote<base_expr_type_1,
1269
1271 Expr(const ExprT1& expr1_, const ConstT& c_) :
1272 expr1(expr1_), c(c_) {}
1273
1275 int size() const {
1276 return expr1.size();
1277 }
1278
1280 bool hasFastAccess() const {
1281 return expr1.hasFastAccess();
1282 }
1283
1285 bool isPassive() const {
1286 return expr1.isPassive();
1287 }
1288
1290 bool updateValue() const { return expr1.updateValue(); }
1291
1293 void cache() const {}
1294
1296 const value_type val() const {
1297 using std::pow;
1298 return pow(expr1.val(), c.val());
1299 }
1300
1302 const value_type dx(int i) const {
1303 using std::pow; using Sacado::if_then_else;
1304 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1305 // It seems less accurate and caused convergence problems in some codes
1306 return if_then_else( c.val() == scalar_type(1.0), expr1.dx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val())) ));
1307 }
1308
1310 const value_type fastAccessDx(int i) const {
1311 using std::pow; using Sacado::if_then_else;
1312 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1313 // It seems less accurate and caused convergence problems in some codes
1314 return if_then_else( c.val() == scalar_type(1.0), expr1.fastAccessDx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()))));
1315 }
1316
1317 protected:
1318
1319 const ExprT1& expr1;
1321 };
1322
1323 template <typename T1, typename ExprT2>
1324 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Simd >,
1325 ExprSpecDefault > {
1326
1327 public:
1328
1332 typedef typename ExprT2::value_type value_type_2;
1333 typedef typename Sacado::Promote<value_type_1,
1335
1337 typedef typename ExprT2::scalar_type scalar_type_2;
1338 typedef typename Sacado::Promote<scalar_type_1,
1340
1342 typedef typename ExprT2::base_expr_type base_expr_type_2;
1343 typedef typename Sacado::Promote<base_expr_type_1,
1345
1346
1348 Expr(const ConstT& c_, const ExprT2& expr2_) :
1349 c(c_), expr2(expr2_) {}
1350
1352 int size() const {
1353 return expr2.size();
1354 }
1355
1357 bool hasFastAccess() const {
1358 return expr2.hasFastAccess();
1359 }
1360
1362 bool isPassive() const {
1363 return expr2.isPassive();
1364 }
1365
1367 bool updateValue() const { return expr2.updateValue(); }
1368
1370 void cache() const {}
1371
1373 const value_type val() const {
1374 using std::pow;
1375 return pow(c.val(), expr2.val());
1376 }
1377
1379 const value_type dx(int i) const {
1380 using std::pow; using std::log; using Sacado::if_then_else;
1381 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val())) );
1382 }
1383
1385 const value_type fastAccessDx(int i) const {
1386 using std::pow; using std::log; using Sacado::if_then_else;
1387 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val())) );
1388 }
1389
1390 protected:
1391
1393 const ExprT2& expr2;
1394 };
1395
1396 //
1397 // Specialization for scalar types using ternary operator
1398 //
1399
1400 template <typename ExprT1, typename ExprT2>
1401 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Scalar >,
1402 ExprSpecDefault > {
1403
1404 public:
1405
1406 typedef typename ExprT1::value_type value_type_1;
1407 typedef typename ExprT2::value_type value_type_2;
1408 typedef typename Sacado::Promote<value_type_1,
1410
1411 typedef typename ExprT1::scalar_type scalar_type_1;
1412 typedef typename ExprT2::scalar_type scalar_type_2;
1413 typedef typename Sacado::Promote<scalar_type_1,
1415
1416 typedef typename ExprT1::base_expr_type base_expr_type_1;
1417 typedef typename ExprT2::base_expr_type base_expr_type_2;
1418 typedef typename Sacado::Promote<base_expr_type_1,
1420
1422 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1423 expr1(expr1_), expr2(expr2_) {}
1424
1426 int size() const {
1427 int sz1 = expr1.size(), sz2 = expr2.size();
1428 return sz1 > sz2 ? sz1 : sz2;
1429 }
1430
1432 bool hasFastAccess() const {
1433 return expr1.hasFastAccess() && expr2.hasFastAccess();
1434 }
1435
1437 bool isPassive() const {
1438 return expr1.isPassive() && expr2.isPassive();
1439 }
1440
1442 bool updateValue() const {
1443 return expr1.updateValue() && expr2.updateValue();
1444 }
1445
1447 void cache() const {}
1448
1450 const value_type val() const {
1451 using std::pow;
1452 return pow(expr1.val(), expr2.val());
1453 }
1454
1456 const value_type dx(int i) const {
1457 using std::pow; using std::log;
1458 const int sz1 = expr1.size(), sz2 = expr2.size();
1459 if (sz1 > 0 && sz2 > 0)
1460 return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val()));
1461 else if (sz1 > 0)
1462 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1463 // It seems less accurate and caused convergence problems in some codes
1464 return expr2.val() == scalar_type(1.0) ? expr1.dx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val()));
1465 else
1466 return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
1467 }
1468
1470 const value_type fastAccessDx(int i) const {
1471 using std::pow; using std::log;
1472 return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type((expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val()));
1473 }
1474
1475 protected:
1476
1477 const ExprT1& expr1;
1478 const ExprT2& expr2;
1479
1480 };
1481
1482 template <typename ExprT1, typename T2>
1483 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Scalar >,
1484 ExprSpecDefault > {
1485
1486 public:
1487
1490 typedef typename ExprT1::value_type value_type_1;
1492 typedef typename Sacado::Promote<value_type_1,
1494
1495 typedef typename ExprT1::scalar_type scalar_type_1;
1497 typedef typename Sacado::Promote<scalar_type_1,
1499
1500 typedef typename ExprT1::base_expr_type base_expr_type_1;
1502 typedef typename Sacado::Promote<base_expr_type_1,
1504
1506 Expr(const ExprT1& expr1_, const ConstT& c_) :
1507 expr1(expr1_), c(c_) {}
1508
1510 int size() const {
1511 return expr1.size();
1512 }
1513
1515 bool hasFastAccess() const {
1516 return expr1.hasFastAccess();
1517 }
1518
1520 bool isPassive() const {
1521 return expr1.isPassive();
1522 }
1523
1525 bool updateValue() const { return expr1.updateValue(); }
1526
1528 void cache() const {}
1529
1531 const value_type val() const {
1532 using std::pow;
1533 return pow(expr1.val(), c.val());
1534 }
1535
1537 const value_type dx(int i) const {
1538 using std::pow;
1539 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1540 // It seems less accurate and caused convergence problems in some codes
1541 return c.val() == scalar_type(1.0) ? expr1.dx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val()));
1542 }
1543
1545 const value_type fastAccessDx(int i) const {
1546 using std::pow;
1547 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1548 // It seems less accurate and caused convergence problems in some codes
1549 return c.val() == scalar_type(1.0) ? expr1.fastAccessDx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()));
1550 }
1551
1552 protected:
1553
1554 const ExprT1& expr1;
1556 };
1557
1558 template <typename T1, typename ExprT2>
1559 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Scalar >,
1560 ExprSpecDefault > {
1561
1562 public:
1563
1567 typedef typename ExprT2::value_type value_type_2;
1568 typedef typename Sacado::Promote<value_type_1,
1570
1572 typedef typename ExprT2::scalar_type scalar_type_2;
1573 typedef typename Sacado::Promote<scalar_type_1,
1575
1577 typedef typename ExprT2::base_expr_type base_expr_type_2;
1578 typedef typename Sacado::Promote<base_expr_type_1,
1580
1581
1583 Expr(const ConstT& c_, const ExprT2& expr2_) :
1584 c(c_), expr2(expr2_) {}
1585
1587 int size() const {
1588 return expr2.size();
1589 }
1590
1592 bool hasFastAccess() const {
1593 return expr2.hasFastAccess();
1594 }
1595
1597 bool isPassive() const {
1598 return expr2.isPassive();
1599 }
1600
1602 bool updateValue() const { return expr2.updateValue(); }
1603
1605 void cache() const {}
1606
1608 const value_type val() const {
1609 using std::pow;
1610 return pow(c.val(), expr2.val());
1611 }
1612
1614 const value_type dx(int i) const {
1615 using std::pow; using std::log;
1616 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val()));
1617 }
1618
1620 const value_type fastAccessDx(int i) const {
1621 using std::pow; using std::log;
1622 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val()));
1623 }
1624
1625 protected:
1626
1628 const ExprT2& expr2;
1629 };
1630
1631 //
1632 // Specialization for nested derivatives. This version does not use
1633 // if_then_else/ternary-operator on the base so that nested derivatives
1634 // are correct.
1635 //
1636 template <typename ExprT1, typename ExprT2>
1637 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Nested >,
1638 ExprSpecDefault > {
1639
1640 public:
1641
1642 typedef typename ExprT1::value_type value_type_1;
1643 typedef typename ExprT2::value_type value_type_2;
1644 typedef typename Sacado::Promote<value_type_1,
1646
1647 typedef typename ExprT1::scalar_type scalar_type_1;
1648 typedef typename ExprT2::scalar_type scalar_type_2;
1649 typedef typename Sacado::Promote<scalar_type_1,
1651
1652 typedef typename ExprT1::base_expr_type base_expr_type_1;
1653 typedef typename ExprT2::base_expr_type base_expr_type_2;
1654 typedef typename Sacado::Promote<base_expr_type_1,
1656
1658 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1659 expr1(expr1_), expr2(expr2_) {}
1660
1662 int size() const {
1663 int sz1 = expr1.size(), sz2 = expr2.size();
1664 return sz1 > sz2 ? sz1 : sz2;
1665 }
1666
1668 bool hasFastAccess() const {
1669 return expr1.hasFastAccess() && expr2.hasFastAccess();
1670 }
1671
1673 bool isPassive() const {
1674 return expr1.isPassive() && expr2.isPassive();
1675 }
1676
1678 bool updateValue() const {
1679 return expr1.updateValue() && expr2.updateValue();
1680 }
1681
1683 void cache() const {}
1684
1686 const value_type val() const {
1687 using std::pow;
1688 return pow(expr1.val(), expr2.val());
1689 }
1690
1692 const value_type dx(int i) const {
1693 using std::pow; using std::log;
1694 const int sz1 = expr1.size(), sz2 = expr2.size();
1695 if (sz1 > 0 && sz2 > 0)
1696 return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1697 else if (sz1 > 0)
1698 return expr2.val() == scalar_type(0.0) ? value_type(0.0) : value_type((expr2.val()*expr1.dx(i))*pow(expr1.val(),expr2.val()-scalar_type(1.0)));
1699 else
1700 return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1701 }
1702
1704 const value_type fastAccessDx(int i) const {
1705 using std::pow; using std::log;
1706 return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1707 }
1708
1709 protected:
1710
1711 const ExprT1& expr1;
1712 const ExprT2& expr2;
1713
1714 };
1715
1716 template <typename ExprT1, typename T2>
1717 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Nested >,
1718 ExprSpecDefault > {
1719
1720 public:
1721
1724 typedef typename ExprT1::value_type value_type_1;
1726 typedef typename Sacado::Promote<value_type_1,
1728
1729 typedef typename ExprT1::scalar_type scalar_type_1;
1731 typedef typename Sacado::Promote<scalar_type_1,
1733
1734 typedef typename ExprT1::base_expr_type base_expr_type_1;
1736 typedef typename Sacado::Promote<base_expr_type_1,
1738
1740 Expr(const ExprT1& expr1_, const ConstT& c_) :
1741 expr1(expr1_), c(c_) {}
1742
1744 int size() const {
1745 return expr1.size();
1746 }
1747
1749 bool hasFastAccess() const {
1750 return expr1.hasFastAccess();
1751 }
1752
1754 bool isPassive() const {
1755 return expr1.isPassive();
1756 }
1757
1759 bool updateValue() const { return expr1.updateValue(); }
1760
1762 void cache() const {}
1763
1765 const value_type val() const {
1766 using std::pow;
1767 return pow(expr1.val(), c.val());
1768 }
1769
1771 const value_type dx(int i) const {
1772 using std::pow;
1773 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
1774 }
1775
1777 const value_type fastAccessDx(int i) const {
1778 using std::pow;
1779 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
1780 }
1781
1782 protected:
1783
1784 const ExprT1& expr1;
1786 };
1787
1788 template <typename T1, typename ExprT2>
1789 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Nested >,
1790 ExprSpecDefault > {
1791
1792 public:
1793
1797 typedef typename ExprT2::value_type value_type_2;
1798 typedef typename Sacado::Promote<value_type_1,
1800
1802 typedef typename ExprT2::scalar_type scalar_type_2;
1803 typedef typename Sacado::Promote<scalar_type_1,
1805
1807 typedef typename ExprT2::base_expr_type base_expr_type_2;
1808 typedef typename Sacado::Promote<base_expr_type_1,
1810
1811
1813 Expr(const ConstT& c_, const ExprT2& expr2_) :
1814 c(c_), expr2(expr2_) {}
1815
1817 int size() const {
1818 return expr2.size();
1819 }
1820
1822 bool hasFastAccess() const {
1823 return expr2.hasFastAccess();
1824 }
1825
1827 bool isPassive() const {
1828 return expr2.isPassive();
1829 }
1830
1832 bool updateValue() const { return expr2.updateValue(); }
1833
1835 void cache() const {}
1836
1838 const value_type val() const {
1839 using std::pow;
1840 return pow(c.val(), expr2.val());
1841 }
1842
1844 const value_type dx(int i) const {
1845 using std::pow; using std::log;
1846 return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
1847 }
1848
1850 const value_type fastAccessDx(int i) const {
1851 using std::pow; using std::log;
1852 return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
1853 }
1854
1855 protected:
1856
1858 const ExprT2& expr2;
1859 };
1860
1861 //
1862 // Specialization for nested derivatives. This version does not use
1863 // if_then_else/ternary-operator on the base so that nested derivatives
1864 // are correct.
1865 //
1866 template <typename ExprT1, typename ExprT2>
1867 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::NestedSimd >,
1868 ExprSpecDefault > {
1869
1870 public:
1871
1872 typedef typename ExprT1::value_type value_type_1;
1873 typedef typename ExprT2::value_type value_type_2;
1874 typedef typename Sacado::Promote<value_type_1,
1876
1877 typedef typename ExprT1::scalar_type scalar_type_1;
1878 typedef typename ExprT2::scalar_type scalar_type_2;
1879 typedef typename Sacado::Promote<scalar_type_1,
1881
1882 typedef typename ExprT1::base_expr_type base_expr_type_1;
1883 typedef typename ExprT2::base_expr_type base_expr_type_2;
1884 typedef typename Sacado::Promote<base_expr_type_1,
1886
1888 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1889 expr1(expr1_), expr2(expr2_) {}
1890
1892 int size() const {
1893 int sz1 = expr1.size(), sz2 = expr2.size();
1894 return sz1 > sz2 ? sz1 : sz2;
1895 }
1896
1898 bool hasFastAccess() const {
1899 return expr1.hasFastAccess() && expr2.hasFastAccess();
1900 }
1901
1903 bool isPassive() const {
1904 return expr1.isPassive() && expr2.isPassive();
1905 }
1906
1908 bool updateValue() const {
1909 return expr1.updateValue() && expr2.updateValue();
1910 }
1911
1913 void cache() const {}
1914
1916 const value_type val() const {
1917 using std::pow;
1918 return pow(expr1.val(), expr2.val());
1919 }
1920
1922 const value_type dx(int i) const {
1923 using std::pow; using std::log; using Sacado::if_then_else;
1924 const int sz1 = expr1.size(), sz2 = expr2.size();
1925 if (sz1 > 0 && sz2 > 0)
1926 return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1927 else if (sz1 > 0)
1928 return if_then_else( expr2.val() == scalar_type(0.0), value_type(0.0), value_type((expr2.val()*expr1.dx(i))*pow(expr1.val(),expr2.val()-scalar_type(1.0))));
1929 else
1930 return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1931 }
1932
1934 const value_type fastAccessDx(int i) const {
1935 using std::pow; using std::log;
1936 return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1937 }
1938
1939 protected:
1940
1941 const ExprT1& expr1;
1942 const ExprT2& expr2;
1943
1944 };
1945
1946 template <typename ExprT1, typename T2>
1947 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::NestedSimd >,
1948 ExprSpecDefault > {
1949
1950 public:
1951
1954 typedef typename ExprT1::value_type value_type_1;
1956 typedef typename Sacado::Promote<value_type_1,
1958
1959 typedef typename ExprT1::scalar_type scalar_type_1;
1961 typedef typename Sacado::Promote<scalar_type_1,
1963
1964 typedef typename ExprT1::base_expr_type base_expr_type_1;
1966 typedef typename Sacado::Promote<base_expr_type_1,
1968
1970 Expr(const ExprT1& expr1_, const ConstT& c_) :
1971 expr1(expr1_), c(c_) {}
1972
1974 int size() const {
1975 return expr1.size();
1976 }
1977
1979 bool hasFastAccess() const {
1980 return expr1.hasFastAccess();
1981 }
1982
1984 bool isPassive() const {
1985 return expr1.isPassive();
1986 }
1987
1989 bool updateValue() const { return expr1.updateValue(); }
1990
1992 void cache() const {}
1993
1995 const value_type val() const {
1996 using std::pow;
1997 return pow(expr1.val(), c.val());
1998 }
1999
2001 const value_type dx(int i) const {
2002 using std::pow; using Sacado::if_then_else;
2003 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))));
2004 }
2005
2007 const value_type fastAccessDx(int i) const {
2008 using std::pow; using Sacado::if_then_else;
2009 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))));
2010 }
2011
2012 protected:
2013
2014 const ExprT1& expr1;
2016 };
2017
2018 template <typename T1, typename ExprT2>
2019 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::NestedSimd >,
2020 ExprSpecDefault > {
2021
2022 public:
2023
2027 typedef typename ExprT2::value_type value_type_2;
2028 typedef typename Sacado::Promote<value_type_1,
2030
2032 typedef typename ExprT2::scalar_type scalar_type_2;
2033 typedef typename Sacado::Promote<scalar_type_1,
2035
2037 typedef typename ExprT2::base_expr_type base_expr_type_2;
2038 typedef typename Sacado::Promote<base_expr_type_1,
2040
2041
2043 Expr(const ConstT& c_, const ExprT2& expr2_) :
2044 c(c_), expr2(expr2_) {}
2045
2047 int size() const {
2048 return expr2.size();
2049 }
2050
2052 bool hasFastAccess() const {
2053 return expr2.hasFastAccess();
2054 }
2055
2057 bool isPassive() const {
2058 return expr2.isPassive();
2059 }
2060
2062 bool updateValue() const { return expr2.updateValue(); }
2063
2065 void cache() const {}
2066
2068 const value_type val() const {
2069 using std::pow;
2070 return pow(c.val(), expr2.val());
2071 }
2072
2074 const value_type dx(int i) const {
2075 using std::pow; using std::log;
2076 return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
2077 }
2078
2080 const value_type fastAccessDx(int i) const {
2081 using std::pow; using std::log;
2082 return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
2083 }
2084
2085 protected:
2086
2088 const ExprT2& expr2;
2089 };
2090
2091 template <typename T1, typename T2>
2093 typename mpl::enable_if_c<
2094 ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value,
2096 >::type
2097 /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(PowerOp)*/
2098 pow (const Expr<T1>& expr1, const Expr<T2>& expr2)
2099 {
2100 typedef PowerOp< Expr<T1>, Expr<T2> > expr_t;
2101
2102 return Expr<expr_t>(expr1, expr2);
2103 }
2104
2105 template <typename T>
2107 Expr< PowerOp< Expr<T>, Expr<T> > >
2108 pow (const Expr<T>& expr1, const Expr<T>& expr2)
2109 {
2110 typedef PowerOp< Expr<T>, Expr<T> > expr_t;
2111
2112 return Expr<expr_t>(expr1, expr2);
2113 }
2114
2115 template <typename T>
2117 Expr< PowerOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
2118 pow (const typename Expr<T>::value_type& c,
2119 const Expr<T>& expr)
2120 {
2122 typedef PowerOp< ConstT, Expr<T> > expr_t;
2123
2124 return Expr<expr_t>(ConstT(c), expr);
2125 }
2126
2127 template <typename T>
2129 Expr< PowerOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
2130 pow (const Expr<T>& expr,
2131 const typename Expr<T>::value_type& c)
2132 {
2134 typedef PowerOp< Expr<T>, ConstT > expr_t;
2135
2136 return Expr<expr_t>(expr, ConstT(c));
2137 }
2138
2139 template <typename T>
2142 pow (const typename Expr<T>::scalar_type& c,
2143 const Expr<T>& expr)
2144 {
2146 typedef PowerOp< ConstT, Expr<T> > expr_t;
2147
2148 return Expr<expr_t>(ConstT(c), expr);
2149 }
2150
2151 template <typename T>
2154 pow (const Expr<T>& expr,
2155 const typename Expr<T>::scalar_type& c)
2156 {
2158 typedef PowerOp< Expr<T>, ConstT > expr_t;
2159
2160 return Expr<expr_t>(expr, ConstT(c));
2161 }
2162 }
2163}
2164
2165//--------------------------if_then_else operator -----------------------
2166// Can't use the above macros because it is a ternary operator (sort of).
2167// Also, relies on C++11
2168
2169
2170namespace Sacado {
2171 namespace Fad {
2172
2173 template <typename CondT, typename ExprT1, typename ExprT2>
2175
2176 template <typename CondT, typename ExprT1, typename ExprT2>
2177 struct ExprSpec< IfThenElseOp< CondT, ExprT1, ExprT2 > > {
2179 };
2180
2181 template <typename CondT, typename ExprT1, typename ExprT2>
2182 class Expr< IfThenElseOp< CondT, ExprT1, ExprT2 >,ExprSpecDefault > {
2183
2184 public:
2185
2186 typedef typename ExprT1::value_type value_type_1;
2187 typedef typename ExprT2::value_type value_type_2;
2188 typedef typename Sacado::Promote<value_type_1,
2190
2191 typedef typename ExprT1::scalar_type scalar_type_1;
2192 typedef typename ExprT2::scalar_type scalar_type_2;
2193 typedef typename Sacado::Promote<scalar_type_1,
2195
2196 typedef typename ExprT1::base_expr_type base_expr_type_1;
2197 typedef typename ExprT2::base_expr_type base_expr_type_2;
2198 typedef typename Sacado::Promote<base_expr_type_1,
2200
2202 Expr(const CondT& cond_, const ExprT1& expr1_, const ExprT2& expr2_) :
2203 cond(cond_), expr1(expr1_), expr2(expr2_) {}
2204
2206 int size() const {
2207 int sz1 = expr1.size(), sz2 = expr2.size();
2208 return sz1 > sz2 ? sz1 : sz2;
2209 }
2210
2212 bool hasFastAccess() const {
2213 return expr1.hasFastAccess() && expr2.hasFastAccess();
2214 }
2215
2217 bool isPassive() const {
2218 return expr1.isPassive() && expr2.isPassive();
2219 }
2220
2222 bool updateValue() const {
2223 return expr1.updateValue() && expr2.updateValue();
2224 }
2225
2227 void cache() const {}
2228
2230 const value_type val() const {
2232 return if_then_else( cond, expr1.val(), expr2.val() );
2233 }
2234
2236 const value_type dx(int i) const {
2238 return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
2239 }
2240
2242 const value_type fastAccessDx(int i) const {
2244 return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
2245 }
2246
2247 protected:
2248
2249 const CondT& cond;
2250 const ExprT1& expr1;
2251 const ExprT2& expr2;
2252
2253 };
2254
2255 template <typename CondT, typename ExprT1, typename T2>
2256 struct ExprSpec< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> > > {
2258 };
2259
2260 template <typename CondT, typename ExprT1, typename T2>
2261 class Expr< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
2262
2263 public:
2264
2267 typedef typename ExprT1::value_type value_type_1;
2269 typedef typename Sacado::Promote<value_type_1,
2271
2272 typedef typename ExprT1::scalar_type scalar_type_1;
2274 typedef typename Sacado::Promote<scalar_type_1,
2276
2277 typedef typename ExprT1::base_expr_type base_expr_type_1;
2279 typedef typename Sacado::Promote<base_expr_type_1,
2281
2283 Expr(const CondT& cond_, const ExprT1& expr1_, const ConstT& c_) :
2284 cond(cond_), expr1(expr1_), c(c_) {}
2285
2287 int size() const {
2288 return expr1.size();
2289 }
2290
2292 bool hasFastAccess() const {
2293 return expr1.hasFastAccess();
2294 }
2295
2297 bool isPassive() const {
2298 return expr1.isPassive();
2299 }
2300
2302 bool updateValue() const { return expr1.updateValue(); }
2303
2305 void cache() const {}
2306
2308 const value_type val() const {
2310 return if_then_else( cond, expr1.val(), c.val() );
2311 }
2312
2314 const value_type dx(int i) const {
2316 return if_then_else( cond, expr1.dx(i), value_type(0.0) );
2317 }
2318
2320 const value_type fastAccessDx(int i) const {
2322 return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
2323 }
2324
2325 protected:
2326
2327 const CondT& cond;
2328 const ExprT1& expr1;
2330 };
2331
2332 template <typename CondT, typename T1, typename ExprT2>
2333 struct ExprSpec< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 > > {
2335 };
2336
2337 template <typename CondT, typename T1, typename ExprT2>
2338 class Expr< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
2339
2340 public:
2341
2345 typedef typename ExprT2::value_type value_type_2;
2346 typedef typename Sacado::Promote<value_type_1,
2348
2350 typedef typename ExprT2::scalar_type scalar_type_2;
2351 typedef typename Sacado::Promote<scalar_type_1,
2353
2355 typedef typename ExprT2::base_expr_type base_expr_type_2;
2356 typedef typename Sacado::Promote<base_expr_type_1,
2358
2360 Expr(const CondT& cond_, const ConstT& c_, const ExprT2& expr2_) :
2361 cond(cond_), c(c_), expr2(expr2_) {}
2362
2364 int size() const {
2365 return expr2.size();
2366 }
2367
2369 bool hasFastAccess() const {
2370 return expr2.hasFastAccess();
2371 }
2372
2374 bool isPassive() const {
2375 return expr2.isPassive();
2376 }
2377
2379 bool updateValue() const { return expr2.updateValue(); }
2380
2382 void cache() const {}
2383
2385 const value_type val() const {
2387 return if_then_else( cond, c.val(), expr2.val() );
2388 }
2389
2391 const value_type dx(int i) const {
2393 return if_then_else( cond, value_type(0.0), expr2.dx(i) );
2394 }
2395
2397 const value_type fastAccessDx(int i) const {
2399 return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
2400 }
2401
2402 protected:
2403
2404 const CondT& cond;
2406 const ExprT2& expr2;
2407 };
2408
2409 template <typename CondT, typename T1, typename T2>
2414 >::type
2415 if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
2416 {
2417 typedef IfThenElseOp< CondT, T1, T2 > expr_t;
2418
2419 return Expr<expr_t>(cond, expr1, expr2);
2420 }
2421
2422 template <typename CondT, typename T>
2424 Expr< IfThenElseOp< CondT, Expr<T>, Expr<T> > >
2425 if_then_else (const CondT& cond, const Expr<T>& expr1, const Expr<T>& expr2)
2426 {
2427 typedef IfThenElseOp< CondT, Expr<T>, Expr<T> > expr_t;
2428
2429 return Expr<expr_t>(cond, expr1, expr2);
2430 }
2431
2432 template <typename CondT, typename T>
2434 Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::value_type>,
2435 Expr<T> > >
2436 if_then_else (const CondT& cond, const typename Expr<T>::value_type& c,
2437 const Expr<T>& expr)
2438 {
2440 typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2441
2442 return Expr<expr_t>(cond, ConstT(c), expr);
2443 }
2444
2445 template <typename CondT, typename T>
2447 Expr< IfThenElseOp< CondT, Expr<T>,
2448 ConstExpr<typename Expr<T>::value_type> > >
2449 if_then_else (const CondT& cond, const Expr<T>& expr,
2450 const typename Expr<T>::value_type& c)
2451 {
2453 typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2454
2455 return Expr<expr_t>(cond, expr, ConstT(c));
2456 }
2457
2458 template <typename CondT, typename T>
2460 typename mpl::disable_if<
2461 std::is_same< typename Expr<T>::value_type,
2462 typename Expr<T>::scalar_type>,
2463 Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::scalar_type>,
2464 Expr<T> > >
2465 >::type
2466 if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
2467 const Expr<T>& expr)
2468 {
2470 typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2471
2472 return Expr<expr_t>(cond, ConstT(c), expr);
2473 }
2474
2475 template <typename CondT, typename T>
2477 typename mpl::disable_if<
2478 std::is_same< typename Expr<T>::value_type,
2479 typename Expr<T>::scalar_type>,
2480 Expr< IfThenElseOp< CondT, Expr<T>,
2481 ConstExpr<typename Expr<T>::scalar_type> > >
2482 >::type
2483 if_then_else (const CondT& cond, const Expr<T>& expr,
2484 const typename Expr<T>::scalar_type& c)
2485 {
2487 typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2488
2489 return Expr<expr_t>(cond, expr, ConstT(c));
2490 }
2491 }
2492}
2493
2494//-------------------------- Relational Operators -----------------------
2495
2496namespace Sacado {
2497 namespace Fad {
2498 template <typename T1, typename T2 = T1>
2500 typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
2501 };
2502 }
2503}
2504
2505#define FAD_RELOP_MACRO(OP) \
2506namespace Sacado { \
2507 namespace Fad { \
2508 template <typename ExprT1, typename ExprT2> \
2509 SACADO_INLINE_FUNCTION \
2510 typename ConditionalReturnType<typename Expr<ExprT1>::value_type, \
2511 typename Expr<ExprT2>::value_type>::type \
2512 operator OP (const Expr<ExprT1>& expr1, \
2513 const Expr<ExprT2>& expr2) \
2514 { \
2515 return expr1.val() OP expr2.val(); \
2516 } \
2517 \
2518 template <typename ExprT2> \
2519 SACADO_INLINE_FUNCTION \
2520 typename ConditionalReturnType<typename Expr<ExprT2>::value_type>::type \
2521 operator OP (const typename Expr<ExprT2>::value_type& a, \
2522 const Expr<ExprT2>& expr2) \
2523 { \
2524 return a OP expr2.val(); \
2525 } \
2526 \
2527 template <typename ExprT1> \
2528 SACADO_INLINE_FUNCTION \
2529 typename ConditionalReturnType<typename Expr<ExprT1>::value_type>::type \
2530 operator OP (const Expr<ExprT1>& expr1, \
2531 const typename Expr<ExprT1>::value_type& b) \
2532 { \
2533 return expr1.val() OP b; \
2534 } \
2535 } \
2536}
2537
2544FAD_RELOP_MACRO(<<=)
2545FAD_RELOP_MACRO(>>=)
2548
2549#undef FAD_RELOP_MACRO
2550
2551namespace Sacado {
2552
2553 namespace Fad {
2554
2555 template <typename ExprT>
2557 bool operator ! (const Expr<ExprT>& expr)
2558 {
2559 return ! expr.val();
2560 }
2561
2562 } // namespace Fad
2563
2564} // namespace Sacado
2565
2566//-------------------------- Boolean Operators -----------------------
2567namespace Sacado {
2568
2569 namespace Fad {
2570
2571 template <typename ExprT>
2573 bool toBool(const Expr<ExprT>& x) {
2574 bool is_zero = (x.val() == 0.0);
2575 for (int i=0; i<x.size(); i++)
2576 is_zero = is_zero && (x.dx(i) == 0.0);
2577 return !is_zero;
2578 }
2579
2580 } // namespace Fad
2581
2582} // namespace Sacado
2583
2584#define FAD_BOOL_MACRO(OP) \
2585namespace Sacado { \
2586 namespace Fad { \
2587 template <typename ExprT1, typename ExprT2> \
2588 SACADO_INLINE_FUNCTION \
2589 bool \
2590 operator OP (const Expr<ExprT1>& expr1, \
2591 const Expr<ExprT2>& expr2) \
2592 { \
2593 return toBool(expr1) OP toBool(expr2); \
2594 } \
2595 \
2596 template <typename ExprT2> \
2597 SACADO_INLINE_FUNCTION \
2598 bool \
2599 operator OP (const typename Expr<ExprT2>::value_type& a, \
2600 const Expr<ExprT2>& expr2) \
2601 { \
2602 return a OP toBool(expr2); \
2603 } \
2604 \
2605 template <typename ExprT1> \
2606 SACADO_INLINE_FUNCTION \
2607 bool \
2608 operator OP (const Expr<ExprT1>& expr1, \
2609 const typename Expr<ExprT1>::value_type& b) \
2610 { \
2611 return toBool(expr1) OP b; \
2612 } \
2613 } \
2614}
2615
2618
2619#undef FAD_BOOL_MACRO
2620
2621//-------------------------- I/O Operators -----------------------
2622
2623namespace Sacado {
2624
2625 namespace Fad {
2626
2627 template <typename ExprT>
2628 std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
2629 os << x.val() << " [";
2630
2631 for (int i=0; i< x.size(); i++) {
2632 os << " " << x.dx(i);
2633 }
2634
2635 os << " ]";
2636 return os;
2637 }
2638
2639 } // namespace Fad
2640
2641} // namespace Sacado
2642
2643#endif // SACADO_FAD_OPS_HPP
#define SACADO_INLINE_FUNCTION
expr expr expr bar false
expr true
fabs(expr.val())
asinh(expr.val())
log(expr.val())
tan(expr.val())
abs(expr.val())
expr expr SinOp
expr expr SinhOp
cos(expr.val())
expr expr SqrtOp
expr expr dx(i)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
expr expr ACosOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
expr expr ATanOp
cosh(expr.val())
expr expr ACoshOp
acos(expr.val())
sin(expr.val())
sinh(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 DivisionOp
expr expr ASinOp
cbrt(expr.val())
log10(expr.val())
expr expr TanhOp
expr expr TanOp
expr expr CoshOp
#define FAD_BINARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, CDX1, CDX2, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
#define FAD_RELOP_MACRO(OP)
expr expr ASinhOp
if_then_else(expr.val() >=0, expr.dx(i), value_type(-expr.dx(i)))
exp(expr.val())
atan(expr.val())
expr expr Log10Op
#define FAD_BOOL_MACRO(OP)
acosh(expr.val())
expr expr expr ExpOp
expr expr AbsOp
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
atan2(expr1.val(), expr2.val())
sqrt(expr.val())
expr expr ATanhOp
tanh(expr.val())
atanh(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 Atan2Op
expr expr CosOp
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
asin(expr.val())
log(expr.val())
#define FAD_BINARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
log10(expr.val())
#define FAD_RELOP_MACRO(OP)
exp(expr.val())
#define FAD_BOOL_MACRO(OP)
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 PowerOp
#define SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP)
#define SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP)
#define T
#define T1(r, f)
#define T2(r, f)
adouble max(const adouble &a, const adouble &b)
adouble min(const adouble &a, const adouble &b)
Constant expression template.
ScalarType< value_type >::type scalar_type
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const CondT &cond_, const ConstT &c_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const CondT &cond_, const ExprT1 &expr1_, const ConstT &c_)
SACADO_INLINE_FUNCTION const value_type fastAccessDx(int i) const
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
SACADO_INLINE_FUNCTION Expr(const CondT &cond_, const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
Wrapper for a generic expression template.
Namespace for forward-mode AD classes.
SimpleFad< ValueT > log(const SimpleFad< ValueT > &a)
std::ostream & operator<<(std::ostream &os, const Expr< ExprT > &x)
SimpleFad< ValueT > sqrt(const SimpleFad< ValueT > &a)
SimpleFad< ValueT > operator*(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
SACADO_INLINE_FUNCTION Expr< SafeSqrtOp< Expr< T > > > safe_sqrt(const Expr< T > &)
SACADO_INLINE_FUNCTION bool toBool(const Expr< ExprT > &x)
SACADO_INLINE_FUNCTION mpl::enable_if_c< IsFadExpr< T1 >::value &&IsFadExpr< T2 >::value &&ExprLevel< T1 >::value==ExprLevel< T2 >::value, Expr< IfThenElseOp< CondT, T1, T2 > > >::type if_then_else(const CondT &cond, const T1 &expr1, const T2 &expr2)
SACADO_INLINE_FUNCTION mpl::enable_if_c< ExprLevel< Expr< T1 > >::value==ExprLevel< Expr< T2 > >::value, Expr< PowerOp< Expr< T1 >, Expr< T2 > > > >::type pow(const Expr< T1 > &expr1, const Expr< T2 > &expr2)
SACADO_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
SACADO_INLINE_FUNCTION T if_then_else(const Cond cond, const T &a, const T &b)
PowExprType< Expr< T1 >, Expr< T2 > >::expr_type pow(const Expr< T1 > &expr1, const Expr< T2 > &expr2)
Log10ExprType< T >::expr_type log10(const Expr< T > &expr)
decltype(std::declval< T1 >()==std::declval< T2 >()) type
Meta-function for determining nesting with an expression.
Base template specification for Promote.