Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_Fad_Exp_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// @HEADER
30
31#ifndef SACADO_FAD_EXP_OPS_HPP
32#define SACADO_FAD_EXP_OPS_HPP
33
34#include <type_traits>
35#include <ostream>
36
40#include "Sacado_cmath.hpp"
41
43
44#if defined(HAVE_SACADO_KOKKOSCORE)
45#include "Kokkos_Atomic.hpp"
46#include "impl/Kokkos_Error.hpp"
47#endif
48
49#define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
50namespace Sacado { \
51 namespace Fad { \
52 namespace Exp { \
53 \
54 template <typename T, typename ExprSpec> \
55 class OP {}; \
56 \
57 template <typename T> \
58 class OP< T,ExprSpecDefault > : \
59 public Expr< OP<T,ExprSpecDefault> > { \
60 public: \
61 \
62 typedef typename std::remove_cv<T>::type ExprT; \
63 typedef typename ExprT::value_type value_type; \
64 typedef typename ExprT::scalar_type scalar_type; \
65 \
66 typedef ExprSpecDefault expr_spec_type; \
67 \
68 SACADO_INLINE_FUNCTION \
69 explicit OP(const T& expr_) : expr(expr_) {} \
70 \
71 SACADO_INLINE_FUNCTION \
72 int size() const { return expr.size(); } \
73 \
74 SACADO_INLINE_FUNCTION \
75 bool hasFastAccess() const { \
76 return expr.hasFastAccess(); \
77 } \
78 \
79 SACADO_INLINE_FUNCTION \
80 value_type val() const { \
81 USING \
82 return VALUE; \
83 } \
84 \
85 SACADO_INLINE_FUNCTION \
86 value_type dx(int i) const { \
87 USING \
88 return DX; \
89 } \
90 \
91 SACADO_INLINE_FUNCTION \
92 value_type fastAccessDx(int i) const { \
93 USING \
94 return FASTACCESSDX; \
95 } \
96 \
97 protected: \
98 \
99 const T& expr; \
100 }; \
101 \
102 template <typename T> \
103 SACADO_INLINE_FUNCTION \
104 OP< typename Expr<T>::derived_type, \
105 typename T::expr_spec_type > \
106 OPNAME (const Expr<T>& expr) \
107 { \
108 typedef OP< typename Expr<T>::derived_type, \
109 typename T::expr_spec_type > expr_t; \
110 \
111 return expr_t(expr.derived()); \
112 } \
113 \
114 template <typename T, typename E> \
115 struct ExprLevel< OP< T,E > > { \
116 static const unsigned value = ExprLevel<T>::value; \
117 }; \
118 \
119 template <typename T, typename E> \
120 struct IsFadExpr< OP< T,E > > { \
121 static const unsigned value = true; \
122 }; \
123 \
124 } \
125 } \
126 \
127 template <typename T, typename E> \
128 struct IsExpr< Fad::Exp::OP< T,E > > { \
129 static const bool value = true; \
130 }; \
131 \
132 template <typename T, typename E> \
133 struct BaseExprType< Fad::Exp::OP< T,E > > { \
134 typedef typename BaseExprType<T>::type type; \
135 }; \
136 \
137 template <typename T, typename E> \
138 struct IsSimdType< Fad::Exp::OP< T,E > > { \
139 static const bool value = \
140 IsSimdType< typename Fad::Exp::OP< T,E >::scalar_type >::value; \
141 }; \
142 \
143 template <typename T, typename E> \
144 struct ValueType< Fad::Exp::OP< T,E > > { \
145 typedef typename Fad::Exp::OP< T,E >::value_type type; \
146 }; \
147 \
148 template <typename T, typename E> \
149 struct ScalarType< Fad::Exp::OP< T,E > > { \
150 typedef typename Fad::Exp::OP< T,E >::scalar_type type; \
151 }; \
152 \
153}
154
156 UnaryPlusOp,
157 ;,
158 expr.val(),
159 expr.dx(i),
160 expr.fastAccessDx(i))
161FAD_UNARYOP_MACRO(operator-,
163 ;,
164 -expr.val(),
165 -expr.dx(i),
166 -expr.fastAccessDx(i))
169 using std::exp;,
170 exp(expr.val()),
171 exp(expr.val())*expr.dx(i),
172 exp(expr.val())*expr.fastAccessDx(i))
175 using std::log;,
176 log(expr.val()),
177 expr.dx(i)/expr.val(),
178 expr.fastAccessDx(i)/expr.val())
181 using std::log10; using std::log;,
182 log10(expr.val()),
183 expr.dx(i)/( log(value_type(10))*expr.val()),
184 expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
187 using std::sqrt;,
188 sqrt(expr.val()),
189 expr.dx(i)/(value_type(2)* sqrt(expr.val())),
190 expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
193 using std::cos; using std::sin;,
194 cos(expr.val()),
195 -expr.dx(i)* sin(expr.val()),
196 -expr.fastAccessDx(i)* sin(expr.val()))
199 using std::cos; using std::sin;,
200 sin(expr.val()),
201 expr.dx(i)* cos(expr.val()),
202 expr.fastAccessDx(i)* cos(expr.val()))
205 using std::tan;,
206 tan(expr.val()),
207 expr.dx(i)*
208 (value_type(1)+ tan(expr.val())* tan(expr.val())),
209 expr.fastAccessDx(i)*
210 (value_type(1)+ tan(expr.val())* tan(expr.val())))
213 using std::acos; using std::sqrt;,
214 acos(expr.val()),
215 -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
216 -expr.fastAccessDx(i) /
217 sqrt(value_type(1)-expr.val()*expr.val()))
220 using std::asin; using std::sqrt;,
221 asin(expr.val()),
222 expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
223 expr.fastAccessDx(i) /
224 sqrt(value_type(1)-expr.val()*expr.val()))
227 using std::atan;,
228 atan(expr.val()),
229 expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
230 expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
233 using std::cosh; using std::sinh;,
234 cosh(expr.val()),
235 expr.dx(i)* sinh(expr.val()),
236 expr.fastAccessDx(i)* sinh(expr.val()))
239 using std::cosh; using std::sinh;,
240 sinh(expr.val()),
241 expr.dx(i)* cosh(expr.val()),
242 expr.fastAccessDx(i)* cosh(expr.val()))
245 using std::tanh;,
246 tanh(expr.val()),
247 expr.dx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())),
248 expr.fastAccessDx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())))
251 using std::acosh; using std::sqrt;,
252 acosh(expr.val()),
253 expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
254 (expr.val()+value_type(1))),
255 expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
256 (expr.val()+value_type(1))))
259 using std::asinh; using std::sqrt;,
260 asinh(expr.val()),
261 expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
262 expr.fastAccessDx(i)/ sqrt(value_type(1)+
263 expr.val()*expr.val()))
266 using std::atanh;,
267 atanh(expr.val()),
268 expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
269 expr.fastAccessDx(i)/(value_type(1)-
270 expr.val()*expr.val()))
273 using std::abs; using Sacado::if_then_else;,
274 abs(expr.val()),
275 if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
276 if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
279 using std::fabs; using Sacado::if_then_else;,
280 fabs(expr.val()),
281 if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
282 if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
285 using std::cbrt;,
286 cbrt(expr.val()),
287 expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
288 expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
289
290// Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
291// "simd" value types that use if_then_else(). The only reason for not using
292// if_then_else() always is to avoid evaluating the derivative if the value is
293// zero to avoid throwing FPEs.
294namespace Sacado {
295 namespace Fad {
296 namespace Exp {
297
298 template <typename T, typename ExprSpec, bool is_simd>
299 class SafeSqrtOp {};
300
301 //
302 // Implementation for simd type using if_then_else()
303 //
304 template <typename T>
305 class SafeSqrtOp< T,ExprSpecDefault,true > :
306 public Expr< SafeSqrtOp<T,ExprSpecDefault> > {
307 public:
308
309 typedef typename std::remove_cv<T>::type ExprT;
310 typedef typename ExprT::value_type value_type;
311 typedef typename ExprT::scalar_type scalar_type;
312
313 typedef ExprSpecDefault expr_spec_type;
314
316 explicit SafeSqrtOp(const T& expr_) : expr(expr_) {}
317
319 int size() const { return expr.size(); }
320
322 bool hasFastAccess() const {
323 return expr.hasFastAccess();
324 }
325
327 value_type val() const {
328 using std::sqrt;
329 return sqrt(expr.val());
330 }
331
333 value_type dx(int i) const {
334 using std::sqrt; using Sacado::if_then_else;
335 return if_then_else(
336 expr.val() == value_type(0.0), value_type(0.0),
337 value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
338 }
339
341 value_type fastAccessDx(int i) const {
342 using std::sqrt; using Sacado::if_then_else;
343 return if_then_else(
344 expr.val() == value_type(0.0), value_type(0.0),
345 value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
346 }
347
348 protected:
349
350 const T& expr;
351 };
352
353 //
354 // Specialization for scalar types using ternary operator
355 //
356 template <typename T>
357 class SafeSqrtOp< T,ExprSpecDefault,false > :
358 public Expr< SafeSqrtOp<T,ExprSpecDefault> > {
359 public:
360
361 typedef typename std::remove_cv<T>::type ExprT;
362 typedef typename ExprT::value_type value_type;
363 typedef typename ExprT::scalar_type scalar_type;
364
365 typedef ExprSpecDefault expr_spec_type;
366
368 explicit SafeSqrtOp(const T& expr_) : expr(expr_) {}
369
371 int size() const { return expr.size(); }
372
374 bool hasFastAccess() const {
375 return expr.hasFastAccess();
376 }
377
379 value_type val() const {
380 using std::sqrt;
381 return sqrt(expr.val());
382 }
383
385 value_type dx(int i) const {
386 using std::sqrt;
387 return expr.val() == value_type(0.0) ? value_type(0.0) :
388 value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
389 }
390
392 value_type fastAccessDx(int i) const {
393 using std::sqrt;
394 return expr.val() == value_type(0.0) ? value_type(0.0) :
395 value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
396 }
397
398 protected:
399
400 const T& expr;
401 };
402
403 template <typename T>
405 SafeSqrtOp< typename Expr<T>::derived_type,
406 typename T::expr_spec_type >
407 safe_sqrt (const Expr<T>& expr)
408 {
409 typedef SafeSqrtOp< typename Expr<T>::derived_type,
410 typename T::expr_spec_type > expr_t;
411
412 return expr_t(expr.derived());
413 }
414
415 template <typename T, typename E>
416 struct ExprLevel< SafeSqrtOp< T,E > > {
417 static const unsigned value = ExprLevel<T>::value;
418 };
419
420 template <typename T, typename E>
421 struct IsFadExpr< SafeSqrtOp< T,E > > {
422 static const unsigned value = true;
423 };
424
425 }
426 }
427
428 template <typename T, typename E>
429 struct IsExpr< Fad::Exp::SafeSqrtOp< T,E > > {
430 static const bool value = true;
431 };
432
433 template <typename T, typename E>
434 struct BaseExprType< Fad::Exp::SafeSqrtOp< T,E > > {
435 typedef typename BaseExprType<T>::type type;
436 };
437
438 template <typename T, typename E>
439 struct IsSimdType< Fad::Exp::SafeSqrtOp< T,E > > {
440 static const bool value =
441 IsSimdType< typename Fad::Exp::SafeSqrtOp< T,E >::scalar_type >::value;
442 };
443
444 template <typename T, typename E>
445 struct ValueType< Fad::Exp::SafeSqrtOp< T,E > > {
446 typedef typename Fad::Exp::SafeSqrtOp< T,E >::value_type type;
447 };
448
449 template <typename T, typename E>
450 struct ScalarType< Fad::Exp::SafeSqrtOp< T,E > > {
451 typedef typename Fad::Exp::SafeSqrtOp< T,E >::scalar_type type;
452 };
453
454}
455
456#undef FAD_UNARYOP_MACRO
457
458#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) \
459namespace Sacado { \
460 namespace Fad { \
461 namespace Exp { \
462 \
463 template <typename T1, typename T2, \
464 bool is_const_T1, bool is_const_T2, \
465 typename ExprSpec > \
466 class OP {}; \
467 \
468 template <typename T1, typename T2> \
469 class OP< T1, T2, false, false, ExprSpecDefault > : \
470 public Expr< OP< T1, T2, false, false, ExprSpecDefault > > { \
471 public: \
472 \
473 typedef typename std::remove_cv<T1>::type ExprT1; \
474 typedef typename std::remove_cv<T2>::type ExprT2; \
475 typedef typename ExprT1::value_type value_type_1; \
476 typedef typename ExprT2::value_type value_type_2; \
477 typedef typename Sacado::Promote<value_type_1, \
478 value_type_2>::type value_type; \
479 \
480 typedef typename ExprT1::scalar_type scalar_type_1; \
481 typedef typename ExprT2::scalar_type scalar_type_2; \
482 typedef typename Sacado::Promote<scalar_type_1, \
483 scalar_type_2>::type scalar_type; \
484 \
485 typedef ExprSpecDefault expr_spec_type; \
486 \
487 SACADO_INLINE_FUNCTION \
488 OP(const T1& expr1_, const T2& expr2_) : \
489 expr1(expr1_), expr2(expr2_) {} \
490 \
491 SACADO_INLINE_FUNCTION \
492 int size() const { \
493 const int sz1 = expr1.size(), sz2 = expr2.size(); \
494 return sz1 > sz2 ? sz1 : sz2; \
495 } \
496 \
497 SACADO_INLINE_FUNCTION \
498 bool hasFastAccess() const { \
499 return expr1.hasFastAccess() && expr2.hasFastAccess(); \
500 } \
501 \
502 SACADO_INLINE_FUNCTION \
503 value_type val() const { \
504 USING \
505 return VALUE; \
506 } \
507 \
508 SACADO_INLINE_FUNCTION \
509 value_type dx(int i) const { \
510 USING \
511 const int sz1 = expr1.size(), sz2 = expr2.size(); \
512 if (sz1 > 0 && sz2 > 0) \
513 return DX; \
514 else if (sz1 > 0) \
515 return CDX2; \
516 else \
517 return CDX1; \
518 } \
519 \
520 SACADO_INLINE_FUNCTION \
521 value_type fastAccessDx(int i) const { \
522 USING \
523 return FASTACCESSDX; \
524 } \
525 \
526 protected: \
527 \
528 const T1& expr1; \
529 const T2& expr2; \
530 \
531 }; \
532 \
533 template <typename T1, typename T2> \
534 class OP< T1, T2, false, true, ExprSpecDefault > \
535 : public Expr< OP< T1, T2, false, true, ExprSpecDefault > > { \
536 public: \
537 \
538 typedef typename std::remove_cv<T1>::type ExprT1; \
539 typedef T2 ConstT; \
540 typedef typename ExprT1::value_type value_type; \
541 typedef typename ExprT1::scalar_type scalar_type; \
542 \
543 typedef ExprSpecDefault expr_spec_type; \
544 \
545 SACADO_INLINE_FUNCTION \
546 OP(const T1& expr1_, const ConstT& c_) : \
547 expr1(expr1_), c(c_) {} \
548 \
549 SACADO_INLINE_FUNCTION \
550 int size() const { \
551 return expr1.size(); \
552 } \
553 \
554 SACADO_INLINE_FUNCTION \
555 bool hasFastAccess() const { \
556 return expr1.hasFastAccess(); \
557 } \
558 \
559 SACADO_INLINE_FUNCTION \
560 value_type val() const { \
561 USING \
562 return VAL_CONST_DX_2; \
563 } \
564 \
565 SACADO_INLINE_FUNCTION \
566 value_type dx(int i) const { \
567 USING \
568 return CONST_DX_2; \
569 } \
570 \
571 SACADO_INLINE_FUNCTION \
572 value_type fastAccessDx(int i) const { \
573 USING \
574 return CONST_FASTACCESSDX_2; \
575 } \
576 \
577 protected: \
578 \
579 const T1& expr1; \
580 const ConstT& c; \
581 }; \
582 \
583 template <typename T1, typename T2> \
584 class OP< T1, T2, true, false, ExprSpecDefault > \
585 : public Expr< OP< T1, T2, true, false, ExprSpecDefault > > { \
586 public: \
587 \
588 typedef typename std::remove_cv<T2>::type ExprT2; \
589 typedef T1 ConstT; \
590 typedef typename ExprT2::value_type value_type; \
591 typedef typename ExprT2::scalar_type scalar_type; \
592 \
593 typedef ExprSpecDefault expr_spec_type; \
594 \
595 SACADO_INLINE_FUNCTION \
596 OP(const ConstT& c_, const T2& expr2_) : \
597 c(c_), expr2(expr2_) {} \
598 \
599 SACADO_INLINE_FUNCTION \
600 int size() const { \
601 return expr2.size(); \
602 } \
603 \
604 SACADO_INLINE_FUNCTION \
605 bool hasFastAccess() const { \
606 return expr2.hasFastAccess(); \
607 } \
608 \
609 SACADO_INLINE_FUNCTION \
610 value_type val() const { \
611 USING \
612 return VAL_CONST_DX_1; \
613 } \
614 \
615 SACADO_INLINE_FUNCTION \
616 value_type dx(int i) const { \
617 USING \
618 return CONST_DX_1; \
619 } \
620 \
621 SACADO_INLINE_FUNCTION \
622 value_type fastAccessDx(int i) const { \
623 USING \
624 return CONST_FASTACCESSDX_1; \
625 } \
626 \
627 protected: \
628 \
629 const ConstT& c; \
630 const T2& expr2; \
631 }; \
632 \
633 template <typename T1, typename T2> \
634 SACADO_INLINE_FUNCTION \
635 SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP) \
636 OPNAME (const T1& expr1, const T2& expr2) \
637 { \
638 typedef OP< typename Expr<T1>::derived_type, \
639 typename Expr<T2>::derived_type, \
640 false, false, typename T1::expr_spec_type > expr_t; \
641 \
642 return expr_t(expr1.derived(), expr2.derived()); \
643 } \
644 \
645 template <typename T> \
646 SACADO_INLINE_FUNCTION \
647 OP< typename T::value_type, typename Expr<T>::derived_type, \
648 true, false, typename T::expr_spec_type > \
649 OPNAME (const typename T::value_type& c, \
650 const Expr<T>& expr) \
651 { \
652 typedef typename T::value_type ConstT; \
653 typedef OP< ConstT, typename Expr<T>::derived_type, \
654 true, false, typename T::expr_spec_type > expr_t; \
655 \
656 return expr_t(c, expr.derived()); \
657 } \
658 \
659 template <typename T> \
660 SACADO_INLINE_FUNCTION \
661 OP< typename Expr<T>::derived_type, typename T::value_type, \
662 false, true, typename T::expr_spec_type > \
663 OPNAME (const Expr<T>& expr, \
664 const typename T::value_type& c) \
665 { \
666 typedef typename T::value_type ConstT; \
667 typedef OP< typename Expr<T>::derived_type, ConstT, \
668 false, true, typename T::expr_spec_type > expr_t; \
669 \
670 return expr_t(expr.derived(), c); \
671 } \
672 \
673 template <typename T> \
674 SACADO_INLINE_FUNCTION \
675 SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP) \
676 OPNAME (const typename T::scalar_type& c, \
677 const Expr<T>& expr) \
678 { \
679 typedef typename T::scalar_type ConstT; \
680 typedef OP< ConstT, typename Expr<T>::derived_type, \
681 true, false, typename T::expr_spec_type > expr_t; \
682 \
683 return expr_t(c, expr.derived()); \
684 } \
685 \
686 template <typename T> \
687 SACADO_INLINE_FUNCTION \
688 SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP) \
689 OPNAME (const Expr<T>& expr, \
690 const typename T::scalar_type& c) \
691 { \
692 typedef typename T::scalar_type ConstT; \
693 typedef OP< typename Expr<T>::derived_type, ConstT, \
694 false, true, typename T::expr_spec_type > expr_t; \
695 \
696 return expr_t(expr.derived(), c); \
697 } \
698 \
699 template <typename T1, typename T2, bool c1, bool c2, typename E> \
700 struct ExprLevel< OP< T1, T2, c1, c2, E > > { \
701 static constexpr unsigned value_1 = ExprLevel<T1>::value; \
702 static constexpr unsigned value_2 = ExprLevel<T2>::value; \
703 static constexpr unsigned value = \
704 value_1 >= value_2 ? value_1 : value_2; \
705 }; \
706 \
707 template <typename T1, typename T2, bool c1, bool c2, typename E> \
708 struct IsFadExpr< OP< T1, T2, c1, c2, E > > { \
709 static constexpr unsigned value = true; \
710 }; \
711 \
712 } \
713 } \
714 \
715 template <typename T1, typename T2, bool c1, bool c2, typename E> \
716 struct IsExpr< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
717 static constexpr bool value = true; \
718 }; \
719 \
720 template <typename T1, typename T2, bool c1, bool c2, typename E> \
721 struct BaseExprType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
722 typedef typename BaseExprType<T1>::type base_expr_1; \
723 typedef typename BaseExprType<T2>::type base_expr_2; \
724 typedef typename Sacado::Promote<base_expr_1, \
725 base_expr_2>::type type; \
726 }; \
727 \
728 template <typename T1, typename T2, bool c1, bool c2, typename E> \
729 struct IsSimdType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
730 static const bool value = \
731 IsSimdType< typename Fad::Exp::OP< T1, T2, c1, c2, E >::value_type >::value; \
732 }; \
733 \
734 template <typename T1, typename T2, bool c1, bool c2, typename E> \
735 struct ValueType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
736 typedef typename Fad::Exp::OP< T1, T2, c1, c2, E >::value_type type;\
737 }; \
738 \
739 template <typename T1, typename T2, bool c1, bool c2, typename E> \
740 struct ScalarType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
741 typedef typename Fad::Exp::OP< T1, T2, c1, c2, E >::scalar_type type;\
742 }; \
743 \
744}
745
746
749 ;,
750 expr1.val() + expr2.val(),
751 expr1.dx(i) + expr2.dx(i),
752 expr2.dx(i),
753 expr1.dx(i),
754 expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
755 c + expr2.val(),
756 expr1.val() + c,
757 expr2.dx(i),
758 expr1.dx(i),
759 expr2.fastAccessDx(i),
760 expr1.fastAccessDx(i))
761FAD_BINARYOP_MACRO(operator-,
763 ;,
764 expr1.val() - expr2.val(),
765 expr1.dx(i) - expr2.dx(i),
766 -expr2.dx(i),
767 expr1.dx(i),
768 expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
769 c - expr2.val(),
770 expr1.val() - c,
771 -expr2.dx(i),
772 expr1.dx(i),
773 -expr2.fastAccessDx(i),
774 expr1.fastAccessDx(i))
775FAD_BINARYOP_MACRO(operator*,
777 ;,
778 expr1.val() * expr2.val(),
779 expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
780 expr1.val()*expr2.dx(i),
781 expr1.dx(i)*expr2.val(),
782 expr1.val()*expr2.fastAccessDx(i) +
783 expr1.fastAccessDx(i)*expr2.val(),
784 c * expr2.val(),
785 expr1.val() * c,
786 c*expr2.dx(i),
787 expr1.dx(i)*c,
788 c*expr2.fastAccessDx(i),
789 expr1.fastAccessDx(i)*c)
790FAD_BINARYOP_MACRO(operator/,
792 ;,
793 expr1.val() / expr2.val(),
794 (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
795 (expr2.val()*expr2.val()),
796 -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
797 expr1.dx(i)/expr2.val(),
798 (expr1.fastAccessDx(i)*expr2.val() -
799 expr2.fastAccessDx(i)*expr1.val()) /
800 (expr2.val()*expr2.val()),
801 c / expr2.val(),
802 expr1.val() / c,
803 -expr2.dx(i)*c / (expr2.val()*expr2.val()),
804 expr1.dx(i)/c,
805 -expr2.fastAccessDx(i)*c / (expr2.val()*expr2.val()),
806 expr1.fastAccessDx(i)/c)
809 using std::atan2;,
810 atan2(expr1.val(), expr2.val()),
811 (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
812 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
813 -expr1.val()*expr2.dx(i)/
814 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
815 expr2.val()*expr1.dx(i)/
816 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
817 (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
818 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
819 atan2(c, expr2.val()),
820 atan2(expr1.val(), c),
821 (-c*expr2.dx(i)) / (c*c + expr2.val()*expr2.val()),
822 (c*expr1.dx(i))/ (expr1.val()*expr1.val() + c*c),
823 (-c*expr2.fastAccessDx(i))/ (c*c + expr2.val()*expr2.val()),
824 (c*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c*c))
825// FAD_BINARYOP_MACRO(pow,
826// PowerOp,
827// using std::pow; using std::log; using Sacado::if_then_else;,
828// pow(expr1.val(), expr2.val()),
829// 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())) ),
830// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val())) ),
831// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val())) ),
832// 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())) ),
833// pow(c, expr2.val()),
834// pow(expr1.val(), c),
835// if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) ),
836// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c*expr1.dx(i)/expr1.val()*pow(expr1.val(),c)) ),
837// if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) ),
838// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c))) )
841 using Sacado::if_then_else;,
842 if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
843 if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
844 if_then_else( expr1.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
845 if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), value_type(0.0) ),
846 if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
847 if_then_else( c >= expr2.val(), value_type(c), expr2.val() ),
848 if_then_else( expr1.val() >= c, expr1.val(), value_type(c) ),
849 if_then_else( c >= expr2.val(), value_type(0.0), expr2.dx(i) ),
850 if_then_else( expr1.val() >= c, expr1.dx(i), value_type(0.0) ),
851 if_then_else( c >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
852 if_then_else( expr1.val() >= c, expr1.fastAccessDx(i), value_type(0.0) ) )
855 using Sacado::if_then_else;,
856 if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
857 if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
858 if_then_else( expr1.val() <= expr2.val(), value_type(0.0), expr2.dx(i) ),
859 if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), value_type(0.0) ),
860 if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
861 if_then_else( c <= expr2.val(), value_type(c), expr2.val() ),
862 if_then_else( expr1.val() <= c, expr1.val(), value_type(c) ),
863 if_then_else( c <= expr2.val(), value_type(0), expr2.dx(i) ),
864 if_then_else( expr1.val() <= c, expr1.dx(i), value_type(0) ),
865 if_then_else( c <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
866 if_then_else( expr1.val() <= c, expr1.fastAccessDx(i), value_type(0) ) )
867
868
869// Special handling for std::pow() to provide specializations of PowerOp for
870// "simd" value types that use if_then_else(). The only reason for not using
871// if_then_else() always is to avoid evaluating the derivative if the value is
872// zero to avoid throwing FPEs.
873namespace Sacado {
874 namespace Fad {
875 namespace Exp {
876
877 template <typename T1, typename T2,
878 bool is_const_T1, bool is_const_T2,
879 typename ExprSpec, typename Impl >
880 class PowerOp {};
881
882 //
883 // Implementation for simd type using if_then_else()
884 //
885 template <typename T1, typename T2>
886 class PowerOp< T1, T2, false, false, ExprSpecDefault,
887 PowerImpl::Simd > :
888 public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
889 PowerImpl::Simd > > {
890 public:
891
892 typedef typename std::remove_cv<T1>::type ExprT1;
893 typedef typename std::remove_cv<T2>::type ExprT2;
894 typedef typename ExprT1::value_type value_type_1;
895 typedef typename ExprT2::value_type value_type_2;
896 typedef typename Sacado::Promote<value_type_1,
897 value_type_2>::type value_type;
898
899 typedef typename ExprT1::scalar_type scalar_type_1;
900 typedef typename ExprT2::scalar_type scalar_type_2;
901 typedef typename Sacado::Promote<scalar_type_1,
902 scalar_type_2>::type scalar_type;
903
904 typedef ExprSpecDefault expr_spec_type;
905
907 PowerOp(const T1& expr1_, const T2& expr2_) :
908 expr1(expr1_), expr2(expr2_) {}
909
911 int size() const {
912 const int sz1 = expr1.size(), sz2 = expr2.size();
913 return sz1 > sz2 ? sz1 : sz2;
914 }
915
917 bool hasFastAccess() const {
918 return expr1.hasFastAccess() && expr2.hasFastAccess();
919 }
920
922 value_type val() const {
923 using std::pow;
924 return pow(expr1.val(), expr2.val());
925 }
926
928 value_type dx(int i) const {
929 using std::pow; using std::log; using Sacado::if_then_else;
930 const int sz1 = expr1.size(), sz2 = expr2.size();
931 if (sz1 > 0 && sz2 > 0)
932 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())) );
933 else if (sz1 > 0)
934 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
935 // It seems less accurate and caused convergence problems in some codes
936 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())) ));
937 else
938 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())) );
939 }
940
942 value_type fastAccessDx(int i) const {
943 using std::pow; using std::log; using Sacado::if_then_else;
944 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())));
945 }
946
947 protected:
948
949 const T1& expr1;
950 const T2& expr2;
951
952 };
953
954 template <typename T1, typename T2>
955 class PowerOp< T1, T2, false, true, ExprSpecDefault,
956 PowerImpl::Simd > :
957 public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
958 PowerImpl::Simd > > {
959 public:
960
961 typedef typename std::remove_cv<T1>::type ExprT1;
962 typedef T2 ConstT;
963 typedef typename ExprT1::value_type value_type;
964 typedef typename ExprT1::scalar_type scalar_type;
965
966 typedef ExprSpecDefault expr_spec_type;
967
969 PowerOp(const T1& expr1_, const ConstT& c_) :
970 expr1(expr1_), c(c_) {}
971
973 int size() const {
974 return expr1.size();
975 }
976
978 bool hasFastAccess() const {
979 return expr1.hasFastAccess();
980 }
981
983 value_type val() const {
984 using std::pow;
985 return pow(expr1.val(), c);
986 }
987
989 value_type dx(int i) const {
990 using std::pow; using Sacado::if_then_else;
991 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
992 // It seems less accurate and caused convergence problems in some codes
993 return if_then_else( c == scalar_type(1.0), expr1.dx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c*expr1.dx(i)/expr1.val()*pow(expr1.val(),c)) ));
994 }
995
997 value_type fastAccessDx(int i) const {
998 using std::pow; using Sacado::if_then_else;
999 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1000 // It seems less accurate and caused convergence problems in some codes
1001 return if_then_else( c == scalar_type(1.0), expr1.fastAccessDx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c)) ));
1002 }
1003
1004 protected:
1005
1006 const T1& expr1;
1007 const ConstT& c;
1008 };
1009
1010 template <typename T1, typename T2>
1011 class PowerOp< T1, T2, true, false, ExprSpecDefault,
1012 PowerImpl::Simd > :
1013 public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
1014 PowerImpl::Simd > > {
1015 public:
1016
1017 typedef typename std::remove_cv<T2>::type ExprT2;
1018 typedef T1 ConstT;
1019 typedef typename ExprT2::value_type value_type;
1020 typedef typename ExprT2::scalar_type scalar_type;
1021
1022 typedef ExprSpecDefault expr_spec_type;
1023
1025 PowerOp(const ConstT& c_, const T2& expr2_) :
1026 c(c_), expr2(expr2_) {}
1027
1029 int size() const {
1030 return expr2.size();
1031 }
1032
1034 bool hasFastAccess() const {
1035 return expr2.hasFastAccess();
1036 }
1037
1039 value_type val() const {
1040 using std::pow;
1041 return pow(c, expr2.val());
1042 }
1043
1045 value_type dx(int i) const {
1046 using std::pow; using std::log; using Sacado::if_then_else;
1047 return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) );
1048 }
1049
1051 value_type fastAccessDx(int i) const {
1052 using std::pow; using std::log; using Sacado::if_then_else;
1053 return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) );
1054 }
1055
1056 protected:
1057
1058 const ConstT& c;
1059 const T2& expr2;
1060 };
1061
1062 //
1063 // Specialization for scalar types using ternary operator
1064 //
1065
1066 template <typename T1, typename T2>
1067 class PowerOp< T1, T2, false, false, ExprSpecDefault,
1068 PowerImpl::Scalar > :
1069 public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
1070 PowerImpl::Scalar > > {
1071 public:
1072
1073 typedef typename std::remove_cv<T1>::type ExprT1;
1074 typedef typename std::remove_cv<T2>::type ExprT2;
1075 typedef typename ExprT1::value_type value_type_1;
1076 typedef typename ExprT2::value_type value_type_2;
1077 typedef typename Sacado::Promote<value_type_1,
1078 value_type_2>::type value_type;
1079
1080 typedef typename ExprT1::scalar_type scalar_type_1;
1081 typedef typename ExprT2::scalar_type scalar_type_2;
1082 typedef typename Sacado::Promote<scalar_type_1,
1083 scalar_type_2>::type scalar_type;
1084
1085 typedef ExprSpecDefault expr_spec_type;
1086
1088 PowerOp(const T1& expr1_, const T2& expr2_) :
1089 expr1(expr1_), expr2(expr2_) {}
1090
1092 int size() const {
1093 const int sz1 = expr1.size(), sz2 = expr2.size();
1094 return sz1 > sz2 ? sz1 : sz2;
1095 }
1096
1098 bool hasFastAccess() const {
1099 return expr1.hasFastAccess() && expr2.hasFastAccess();
1100 }
1101
1103 value_type val() const {
1104 using std::pow;
1105 return pow(expr1.val(), expr2.val());
1106 }
1107
1109 value_type dx(int i) const {
1110 using std::pow; using std::log;
1111 const int sz1 = expr1.size(), sz2 = expr2.size();
1112 if (sz1 > 0 && sz2 > 0)
1113 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()));
1114 else if (sz1 > 0)
1115 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1116 // It seems less accurate and caused convergence problems in some codes
1117 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()));
1118 else
1119 return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
1120 }
1121
1123 value_type fastAccessDx(int i) const {
1124 using std::pow; using std::log;
1125 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()));
1126 }
1127
1128 protected:
1129
1130 const T1& expr1;
1131 const T2& expr2;
1132
1133 };
1134
1135 template <typename T1, typename T2>
1136 class PowerOp< T1, T2, false, true, ExprSpecDefault,
1137 PowerImpl::Scalar > :
1138 public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
1139 PowerImpl::Scalar > > {
1140 public:
1141
1142 typedef typename std::remove_cv<T1>::type ExprT1;
1143 typedef T2 ConstT;
1144 typedef typename ExprT1::value_type value_type;
1145 typedef typename ExprT1::scalar_type scalar_type;
1146
1147 typedef ExprSpecDefault expr_spec_type;
1148
1150 PowerOp(const T1& expr1_, const ConstT& c_) :
1151 expr1(expr1_), c(c_) {}
1152
1154 int size() const {
1155 return expr1.size();
1156 }
1157
1159 bool hasFastAccess() const {
1160 return expr1.hasFastAccess();
1161 }
1162
1164 value_type val() const {
1165 using std::pow;
1166 return pow(expr1.val(), c);
1167 }
1168
1170 value_type dx(int i) const {
1171 using std::pow;
1172 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1173 // It seems less accurate and caused convergence problems in some codes
1174 return c == scalar_type(1.0) ? expr1.dx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)/expr1.val()*pow(expr1.val(),c));
1175 }
1176
1178 value_type fastAccessDx(int i) const {
1179 using std::pow;
1180 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1181 // It seems less accurate and caused convergence problems in some codes
1182 return c == scalar_type(1.0) ? expr1.fastAccessDx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c));
1183 }
1184
1185 protected:
1186
1187 const T1& expr1;
1188 const ConstT& c;
1189 };
1190
1191 template <typename T1, typename T2>
1192 class PowerOp< T1, T2, true, false, ExprSpecDefault,
1193 PowerImpl::Scalar > :
1194 public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
1195 PowerImpl::Scalar > > {
1196 public:
1197
1198 typedef typename std::remove_cv<T2>::type ExprT2;
1199 typedef T1 ConstT;
1200 typedef typename ExprT2::value_type value_type;
1201 typedef typename ExprT2::scalar_type scalar_type;
1202
1203 typedef ExprSpecDefault expr_spec_type;
1204
1206 PowerOp(const ConstT& c_, const T2& expr2_) :
1207 c(c_), expr2(expr2_) {}
1208
1210 int size() const {
1211 return expr2.size();
1212 }
1213
1215 bool hasFastAccess() const {
1216 return expr2.hasFastAccess();
1217 }
1218
1220 value_type val() const {
1221 using std::pow;
1222 return pow(c, expr2.val());
1223 }
1224
1226 value_type dx(int i) const {
1227 using std::pow; using std::log;
1228 return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c)*pow(c,expr2.val()));
1229 }
1230
1232 value_type fastAccessDx(int i) const {
1233 using std::pow; using std::log;
1234 return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val()));
1235 }
1236
1237 protected:
1238
1239 const ConstT& c;
1240 const T2& expr2;
1241 };
1242
1243 //
1244 // Specialization for nested derivatives. This version does not use
1245 // if_then_else/ternary-operator on the base so that nested derivatives
1246 // are correct.
1247 //
1248 template <typename T1, typename T2>
1249 class PowerOp< T1, T2, false, false, ExprSpecDefault,
1250 PowerImpl::Nested > :
1251 public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
1252 PowerImpl::Nested > > {
1253 public:
1254
1255 typedef typename std::remove_cv<T1>::type ExprT1;
1256 typedef typename std::remove_cv<T2>::type ExprT2;
1257 typedef typename ExprT1::value_type value_type_1;
1258 typedef typename ExprT2::value_type value_type_2;
1259 typedef typename Sacado::Promote<value_type_1,
1260 value_type_2>::type value_type;
1261
1262 typedef typename ExprT1::scalar_type scalar_type_1;
1263 typedef typename ExprT2::scalar_type scalar_type_2;
1264 typedef typename Sacado::Promote<scalar_type_1,
1265 scalar_type_2>::type scalar_type;
1266
1267 typedef ExprSpecDefault expr_spec_type;
1268
1270 PowerOp(const T1& expr1_, const T2& expr2_) :
1271 expr1(expr1_), expr2(expr2_) {}
1272
1274 int size() const {
1275 const int sz1 = expr1.size(), sz2 = expr2.size();
1276 return sz1 > sz2 ? sz1 : sz2;
1277 }
1278
1280 bool hasFastAccess() const {
1281 return expr1.hasFastAccess() && expr2.hasFastAccess();
1282 }
1283
1285 value_type val() const {
1286 using std::pow;
1287 return pow(expr1.val(), expr2.val());
1288 }
1289
1291 value_type dx(int i) const {
1292 using std::pow; using std::log;
1293 const int sz1 = expr1.size(), sz2 = expr2.size();
1294 if (sz1 > 0 && sz2 > 0)
1295 return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1296 else if (sz1 > 0)
1297 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)));
1298 else
1299 return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1300 }
1301
1303 value_type fastAccessDx(int i) const {
1304 using std::pow; using std::log;
1305 return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1306 }
1307
1308 protected:
1309
1310 const T1& expr1;
1311 const T2& expr2;
1312
1313 };
1314
1315 template <typename T1, typename T2>
1316 class PowerOp< T1, T2, false, true, ExprSpecDefault,
1317 PowerImpl::Nested > :
1318 public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
1319 PowerImpl::Nested > > {
1320 public:
1321
1322 typedef typename std::remove_cv<T1>::type ExprT1;
1323 typedef T2 ConstT;
1324 typedef typename ExprT1::value_type value_type;
1325 typedef typename ExprT1::scalar_type scalar_type;
1326
1327 typedef ExprSpecDefault expr_spec_type;
1328
1330 PowerOp(const T1& expr1_, const ConstT& c_) :
1331 expr1(expr1_), c(c_) {}
1332
1334 int size() const {
1335 return expr1.size();
1336 }
1337
1339 bool hasFastAccess() const {
1340 return expr1.hasFastAccess();
1341 }
1342
1344 value_type val() const {
1345 using std::pow;
1346 return pow(expr1.val(), c);
1347 }
1348
1350 value_type dx(int i) const {
1351 using std::pow;
1352 return c == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)*pow(expr1.val(),c-scalar_type(1.0)));
1353 }
1354
1356 value_type fastAccessDx(int i) const {
1357 using std::pow;
1358 return c == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)*pow(expr1.val(),c-scalar_type(1.0)));
1359 }
1360
1361 protected:
1362
1363 const T1& expr1;
1364 const ConstT& c;
1365 };
1366
1367 template <typename T1, typename T2>
1368 class PowerOp<T1, T2, true, false, ExprSpecDefault,
1369 PowerImpl::Nested > :
1370 public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
1371 PowerImpl::Nested > > {
1372 public:
1373
1374 typedef typename std::remove_cv<T2>::type ExprT2;
1375 typedef T1 ConstT;
1376 typedef typename ExprT2::value_type value_type;
1377 typedef typename ExprT2::scalar_type scalar_type;
1378
1379 typedef ExprSpecDefault expr_spec_type;
1380
1382 PowerOp(const ConstT& c_, const T2& expr2_) :
1383 c(c_), expr2(expr2_) {}
1384
1386 int size() const {
1387 return expr2.size();
1388 }
1389
1391 bool hasFastAccess() const {
1392 return expr2.hasFastAccess();
1393 }
1394
1396 value_type val() const {
1397 using std::pow;
1398 return pow(c, expr2.val());
1399 }
1400
1402 value_type dx(int i) const {
1403 using std::pow; using std::log;
1404 return expr2.dx(i)*log(c)*pow(c,expr2.val());
1405 }
1406
1408 value_type fastAccessDx(int i) const {
1409 using std::pow; using std::log;
1410 return expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val());
1411 }
1412
1413 protected:
1414
1415 const ConstT& c;
1416 const T2& expr2;
1417 };
1418
1419 //
1420 // Specialization for nested derivatives. This version does not use
1421 // if_then_else/ternary-operator on the base so that nested derivatives
1422 // are correct.
1423 //
1424 template <typename T1, typename T2>
1425 class PowerOp< T1, T2, false, false, ExprSpecDefault,
1426 PowerImpl::NestedSimd > :
1427 public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
1428 PowerImpl::NestedSimd > > {
1429 public:
1430
1431 typedef typename std::remove_cv<T1>::type ExprT1;
1432 typedef typename std::remove_cv<T2>::type ExprT2;
1433 typedef typename ExprT1::value_type value_type_1;
1434 typedef typename ExprT2::value_type value_type_2;
1435 typedef typename Sacado::Promote<value_type_1,
1436 value_type_2>::type value_type;
1437
1438 typedef typename ExprT1::scalar_type scalar_type_1;
1439 typedef typename ExprT2::scalar_type scalar_type_2;
1440 typedef typename Sacado::Promote<scalar_type_1,
1441 scalar_type_2>::type scalar_type;
1442
1443 typedef ExprSpecDefault expr_spec_type;
1444
1446 PowerOp(const T1& expr1_, const T2& expr2_) :
1447 expr1(expr1_), expr2(expr2_) {}
1448
1450 int size() const {
1451 const int sz1 = expr1.size(), sz2 = expr2.size();
1452 return sz1 > sz2 ? sz1 : sz2;
1453 }
1454
1456 bool hasFastAccess() const {
1457 return expr1.hasFastAccess() && expr2.hasFastAccess();
1458 }
1459
1461 value_type val() const {
1462 using std::pow;
1463 return pow(expr1.val(), expr2.val());
1464 }
1465
1467 value_type dx(int i) const {
1468 using std::pow; using std::log; using Sacado::if_then_else;
1469 const int sz1 = expr1.size(), sz2 = expr2.size();
1470 if (sz1 > 0 && sz2 > 0)
1471 return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1472 else if (sz1 > 0)
1473 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))));
1474 else
1475 return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1476 }
1477
1479 value_type fastAccessDx(int i) const {
1480 using std::pow; using std::log; using Sacado::if_then_else;
1481 return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1482 }
1483
1484 protected:
1485
1486 const T1& expr1;
1487 const T2& expr2;
1488
1489 };
1490
1491 template <typename T1, typename T2>
1492 class PowerOp< T1, T2, false, true, ExprSpecDefault,
1493 PowerImpl::NestedSimd > :
1494 public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
1495 PowerImpl::NestedSimd > > {
1496 public:
1497
1498 typedef typename std::remove_cv<T1>::type ExprT1;
1499 typedef T2 ConstT;
1500 typedef typename ExprT1::value_type value_type;
1501 typedef typename ExprT1::scalar_type scalar_type;
1502
1503 typedef ExprSpecDefault expr_spec_type;
1504
1506 PowerOp(const T1& 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 value_type val() const {
1521 using std::pow;
1522 return pow(expr1.val(), c);
1523 }
1524
1526 value_type dx(int i) const {
1527 using std::pow; using Sacado::if_then_else;
1528 return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(c*expr1.dx(i)*pow(expr1.val(),c-scalar_type(1.0))));
1529 }
1530
1532 value_type fastAccessDx(int i) const {
1533 using std::pow; using Sacado::if_then_else;
1534 return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(c*expr1.fastAccessDx(i)*pow(expr1.val(),c-scalar_type(1.0))));
1535 }
1536
1537 protected:
1538
1539 const T1& expr1;
1540 const ConstT& c;
1541 };
1542
1543 template <typename T1, typename T2>
1544 class PowerOp<T1, T2, true, false, ExprSpecDefault,
1545 PowerImpl::NestedSimd > :
1546 public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
1547 PowerImpl::NestedSimd > > {
1548 public:
1549
1550 typedef typename std::remove_cv<T2>::type ExprT2;
1551 typedef T1 ConstT;
1552 typedef typename ExprT2::value_type value_type;
1553 typedef typename ExprT2::scalar_type scalar_type;
1554
1555 typedef ExprSpecDefault expr_spec_type;
1556
1558 PowerOp(const ConstT& c_, const T2& expr2_) :
1559 c(c_), expr2(expr2_) {}
1560
1562 int size() const {
1563 return expr2.size();
1564 }
1565
1567 bool hasFastAccess() const {
1568 return expr2.hasFastAccess();
1569 }
1570
1572 value_type val() const {
1573 using std::pow;
1574 return pow(c, expr2.val());
1575 }
1576
1578 value_type dx(int i) const {
1579 using std::pow; using std::log;
1580 return expr2.dx(i)*log(c)*pow(c,expr2.val());
1581 }
1582
1584 value_type fastAccessDx(int i) const {
1585 using std::pow; using std::log;
1586 return expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val());
1587 }
1588
1589 protected:
1590
1591 const ConstT& c;
1592 const T2& expr2;
1593 };
1594
1595 template <typename T1, typename T2>
1598 pow (const T1& expr1, const T2& expr2)
1599 {
1600 typedef PowerOp< typename Expr<T1>::derived_type,
1601 typename Expr<T2>::derived_type,
1602 false, false, typename T1::expr_spec_type > expr_t;
1603
1604 return expr_t(expr1.derived(), expr2.derived());
1605 }
1606
1607 template <typename T>
1609 PowerOp< typename T::value_type, typename Expr<T>::derived_type,
1610 true, false, typename T::expr_spec_type >
1611 pow (const typename T::value_type& c,
1612 const Expr<T>& expr)
1613 {
1614 typedef typename T::value_type ConstT;
1615 typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1616 true, false, typename T::expr_spec_type > expr_t;
1617
1618 return expr_t(c, expr.derived());
1619 }
1620
1621 template <typename T>
1623 PowerOp< typename Expr<T>::derived_type, typename T::value_type,
1624 false, true, typename T::expr_spec_type >
1625 pow (const Expr<T>& expr,
1626 const typename T::value_type& c)
1627 {
1628 typedef typename T::value_type ConstT;
1629 typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1630 false, true, typename T::expr_spec_type > expr_t;
1631
1632 return expr_t(expr.derived(), c);
1633 }
1634
1635 template <typename T>
1638 pow (const typename T::scalar_type& c,
1639 const Expr<T>& expr)
1640 {
1641 typedef typename T::scalar_type ConstT;
1642 typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1643 true, false, typename T::expr_spec_type > expr_t;
1644
1645 return expr_t(c, expr.derived());
1646 }
1647
1648 template <typename T>
1651 pow (const Expr<T>& expr,
1652 const typename T::scalar_type& c)
1653 {
1654 typedef typename T::scalar_type ConstT;
1655 typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1656 false, true, typename T::expr_spec_type > expr_t;
1657
1658 return expr_t(expr.derived(), c);
1659 }
1660
1661 template <typename T1, typename T2, bool c1, bool c2, typename E>
1662 struct ExprLevel< PowerOp< T1, T2, c1, c2, E > > {
1663 static constexpr unsigned value_1 = ExprLevel<T1>::value;
1664 static constexpr unsigned value_2 = ExprLevel<T2>::value;
1665 static constexpr unsigned value =
1666 value_1 >= value_2 ? value_1 : value_2;
1667 };
1668
1669 template <typename T1, typename T2, bool c1, bool c2, typename E>
1670 struct IsFadExpr< PowerOp< T1, T2, c1, c2, E > > {
1671 static constexpr unsigned value = true;
1672 };
1673
1674 }
1675 }
1676
1677 template <typename T1, typename T2, bool c1, bool c2, typename E>
1678 struct IsExpr< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1679 static constexpr bool value = true;
1680 };
1681
1682 template <typename T1, typename T2, bool c1, bool c2, typename E>
1683 struct BaseExprType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1684 typedef typename BaseExprType<T1>::type base_expr_1;
1685 typedef typename BaseExprType<T2>::type base_expr_2;
1686 typedef typename Sacado::Promote<base_expr_1,
1687 base_expr_2>::type type;
1688 };
1689
1690 template <typename T1, typename T2, bool c1, bool c2, typename E>
1691 struct IsSimdType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1692 static const bool value =
1693 IsSimdType< typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::value_type >::value;
1694 };
1695
1696 template <typename T1, typename T2, bool c1, bool c2, typename E>
1697 struct ValueType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1698 typedef typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::value_type type;
1699 };
1700
1701 template <typename T1, typename T2, bool c1, bool c2, typename E>
1702 struct ScalarType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1703 typedef typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::scalar_type type;
1704 };
1705
1706}
1707
1708//--------------------------if_then_else operator -----------------------
1709// Can't use the above macros because it is a ternary operator (sort of).
1710// Also, relies on C++11
1711
1712namespace Sacado {
1713 namespace Fad {
1714 namespace Exp {
1715
1716 template <typename CondT, typename T1, typename T2,
1717 bool is_const_T1, bool is_const_T2,
1718 typename ExprSpec>
1720
1721 template <typename CondT, typename T1, typename T2>
1723 public Expr< IfThenElseOp< CondT, T1, T2, false, false, ExprSpecDefault > > {
1724
1725 public:
1726
1727 typedef typename std::remove_cv<T1>::type ExprT1;
1728 typedef typename std::remove_cv<T2>::type ExprT2;
1729 typedef typename ExprT1::value_type value_type_1;
1730 typedef typename ExprT2::value_type value_type_2;
1731 typedef typename Sacado::Promote<value_type_1,
1733
1734 typedef typename ExprT1::scalar_type scalar_type_1;
1735 typedef typename ExprT2::scalar_type scalar_type_2;
1736 typedef typename Sacado::Promote<scalar_type_1,
1738
1740
1742 IfThenElseOp(const CondT& cond_, const T1& expr1_, const T2& expr2_) :
1743 cond(cond_), expr1(expr1_), expr2(expr2_) {}
1744
1746 int size() const {
1747 int sz1 = expr1.size(), sz2 = expr2.size();
1748 return sz1 > sz2 ? sz1 : sz2;
1749 }
1750
1752 bool hasFastAccess() const {
1753 return expr1.hasFastAccess() && expr2.hasFastAccess();
1754 }
1755
1757 value_type val() const {
1759 return if_then_else( cond, expr1.val(), expr2.val() );
1760 }
1761
1763 value_type dx(int i) const {
1765 return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
1766 }
1767
1771 return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
1772 }
1773
1774 protected:
1775
1776 const CondT& cond;
1777 const T1& expr1;
1778 const T2& expr2;
1779
1780 };
1781
1782 template <typename CondT, typename T1, typename T2>
1784 public Expr< IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > > {
1785
1786 public:
1787
1788 typedef typename std::remove_cv<T1>::type ExprT1;
1789 typedef T2 ConstT;
1790 typedef typename ExprT1::value_type value_type;
1791 typedef typename ExprT1::scalar_type scalar_type;
1792
1794
1796 IfThenElseOp(const CondT& cond_, const T1& expr1_, const ConstT& c_) :
1797 cond(cond_), expr1(expr1_), c(c_) {}
1798
1800 int size() const {
1801 return expr1.size();
1802 }
1803
1805 bool hasFastAccess() const {
1806 return expr1.hasFastAccess();
1807 }
1808
1810 value_type val() const {
1812 return if_then_else( cond, expr1.val(), c );
1813 }
1814
1816 value_type dx(int i) const {
1818 return if_then_else( cond, expr1.dx(i), value_type(0.0) );
1819 }
1820
1824 return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
1825 }
1826
1827 protected:
1828
1829 const CondT& cond;
1830 const T1& expr1;
1831 const ConstT& c;
1832 };
1833
1834 template <typename CondT, typename T1, typename T2>
1836 public Expr< IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > > {
1837
1838 public:
1839
1840 typedef typename std::remove_cv<T2>::type ExprT2;
1841 typedef T1 ConstT;
1842 typedef typename ExprT2::value_type value_type;
1843 typedef typename ExprT2::scalar_type scalar_type;
1844
1846
1848 IfThenElseOp(const CondT& cond_, const ConstT& c_, const T2& expr2_) :
1849 cond(cond_), c(c_), expr2(expr2_) {}
1850
1852 int size() const {
1853 return expr2.size();
1854 }
1855
1857 bool hasFastAccess() const {
1858 return expr2.hasFastAccess();
1859 }
1860
1862 value_type val() const {
1864 return if_then_else( cond, c, expr2.val() );
1865 }
1866
1868 value_type dx(int i) const {
1870 return if_then_else( cond, value_type(0.0), expr2.dx(i) );
1871 }
1872
1876 return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
1877 }
1878
1879 protected:
1880
1881 const CondT& cond;
1882 const ConstT& c;
1883 const T2& expr2;
1884 };
1885
1886 template <typename CondT, typename T1, typename T2>
1891 typename Expr<T1>::derived_type,
1892 typename Expr<T2>::derived_type,
1893 false, false,
1894 typename T1::expr_spec_type >
1895 >::type
1896 if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
1897 {
1899 typename Expr<T2>::derived_type,
1900 false, false, typename T1::expr_spec_type > expr_t;
1901
1902 return expr_t(cond, expr1.derived(), expr2.derived());
1903 }
1904
1905 template <typename CondT, typename T>
1908 true, false, typename T::expr_spec_type >
1909 if_then_else (const CondT& cond, const typename T::value_type& c,
1910 const Expr<T>& expr)
1911 {
1912 typedef typename T::value_type ConstT;
1914 true, false, typename T::expr_spec_type > expr_t;
1915
1916 return expr_t(cond, c, expr.derived());
1917 }
1918
1919 template <typename CondT, typename T>
1922 false, true, typename T::expr_spec_type >
1923 if_then_else (const CondT& cond, const Expr<T>& expr,
1924 const typename T::value_type& c)
1925 {
1926 typedef typename T::value_type ConstT;
1928 false, true, typename T::expr_spec_type > expr_t;
1929
1930 return expr_t(cond, expr.derived(), c);
1931 }
1932
1933 template <typename CondT, typename T>
1935 typename mpl::disable_if<
1936 mpl::is_same< typename T::value_type,
1937 typename T::scalar_type >,
1938 IfThenElseOp< CondT, typename T::scalar_type,
1939 typename Expr<T>::derived_type,
1940 true, false, typename T::expr_spec_type >
1941 >::type
1942 if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
1943 const Expr<T>& expr)
1944 {
1945 typedef typename T::scalar_type ConstT;
1947 true, false, typename T::expr_spec_type > expr_t;
1948
1949 return expr_t(cond, c, expr.derived());
1950 }
1951
1952 template <typename CondT, typename T>
1954 typename mpl::disable_if<
1955 mpl::is_same< typename T::value_type,
1956 typename T::scalar_type >,
1958 typename T::scalar_type,
1959 false, true, typename T::expr_spec_type >
1960 >::type
1961 if_then_else (const CondT& cond, const Expr<T>& expr,
1962 const typename Expr<T>::scalar_type& c)
1963 {
1964 typedef typename T::scalar_type ConstT;
1966 false, true, typename T::expr_spec_type > expr_t;
1967
1968 return expr_t(cond, expr.derived(), c);
1969 }
1970
1971 template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1972 typename E>
1973 struct ExprLevel< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1974 static constexpr unsigned value_1 = ExprLevel<T1>::value;
1975 static constexpr unsigned value_2 = ExprLevel<T2>::value;
1976 static constexpr unsigned value =
1977 value_1 >= value_2 ? value_1 : value_2;
1978 };
1979
1980 template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1981 typename E>
1982 struct IsFadExpr< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1983 static constexpr unsigned value = true;
1984 };
1985
1986 }
1987 }
1988
1989 template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1990 typename E>
1991 struct IsExpr< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1992 static constexpr bool value = true;
1993 };
1994
1995 template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1996 typename E>
1997 struct BaseExprType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
2000 typedef typename Sacado::Promote<base_expr_1,
2002 };
2003
2004 template <typename CondT, typename T1, typename T2, bool c1, bool c2,
2005 typename E>
2006 struct IsSimdType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
2007 static const bool value =
2009 };
2010
2011 template <typename CondT, typename T1, typename T2, bool c1, bool c2,
2012 typename E>
2013 struct ValueType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
2015 };
2016
2017 template <typename CondT, typename T1, typename T2, bool c1, bool c2,
2018 typename E>
2019 struct ScalarType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
2021 };
2022}
2023
2024#undef FAD_BINARYOP_MACRO
2025
2026//-------------------------- Relational Operators -----------------------
2027
2028namespace Sacado {
2029 namespace Fad {
2030 namespace Exp {
2031 namespace Impl {
2032 // Helper trait to determine return type of logical comparison operations
2033 // (==, !=, ...), usually bool but maybe something else for SIMD types.
2034 // Need to use SFINAE to restrict to types that define == operator in the
2035 // conditional overloads below, otherwise instantiating ConditionaReturnType
2036 // may fail during overload resolution.
2037 template <typename T1, typename T2 = T1,
2040
2041 template <typename T1, typename T2>
2043 typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
2044 };
2045 }
2046 }
2047 }
2048}
2049
2050#define FAD_RELOP_MACRO(OP) \
2051namespace Sacado { \
2052 namespace Fad { \
2053 namespace Exp { \
2054 template <typename T1, typename T2> \
2055 SACADO_INLINE_FUNCTION \
2056 typename mpl::enable_if_c< \
2057 IsFadExpr<T1>::value && IsFadExpr<T2>::value && \
2058 ExprLevel<T1>::value == ExprLevel<T2>::value, \
2059 typename Impl::ConditionalReturnType<typename T1::value_type, \
2060 typename T2::value_type>::type \
2061 >::type \
2062 operator OP (const T1& expr1, const T2& expr2) \
2063 { \
2064 return expr1.derived().val() OP expr2.derived().val(); \
2065 } \
2066 \
2067 template <typename T2> \
2068 SACADO_INLINE_FUNCTION \
2069 typename Impl::ConditionalReturnType<typename T2::value_type>::type \
2070 operator OP (const typename T2::value_type& a, \
2071 const Expr<T2>& expr2) \
2072 { \
2073 return a OP expr2.derived().val(); \
2074 } \
2075 \
2076 template <typename T1> \
2077 SACADO_INLINE_FUNCTION \
2078 typename Impl::ConditionalReturnType<typename T1::value_type>::type \
2079 operator OP (const Expr<T1>& expr1, \
2080 const typename T1::value_type& b) \
2081 { \
2082 return expr1.derived().val() OP b; \
2083 } \
2084 } \
2085 } \
2086} \
2087
2094FAD_RELOP_MACRO(<<=)
2095FAD_RELOP_MACRO(>>=)
2098
2099#undef FAD_RELOP_MACRO
2100
2101namespace Sacado {
2102
2103 namespace Fad {
2104 namespace Exp {
2105
2106 template <typename ExprT>
2108 bool operator ! (const Expr<ExprT>& expr)
2109 {
2110 return ! expr.derived().val();
2111 }
2112
2113 } // namespace Exp
2114 } // namespace Fad
2115
2116} // namespace Sacado
2117
2118//-------------------------- Boolean Operators -----------------------
2119namespace Sacado {
2120
2121 namespace Fad {
2122 namespace Exp {
2123
2124 template <typename T>
2126 bool toBool(const Expr<T>& xx) {
2127 const typename Expr<T>::derived_type& x = xx.derived();
2128 bool is_zero = (x.val() == 0.0);
2129 for (int i=0; i<x.size(); i++)
2130 is_zero = is_zero && (x.dx(i) == 0.0);
2131 return !is_zero;
2132 }
2133
2134 } // namespace Exp
2135 } // namespace Fad
2136
2137} // namespace Sacado
2138
2139#define FAD_BOOL_MACRO(OP) \
2140namespace Sacado { \
2141 namespace Fad { \
2142 namespace Exp { \
2143 template <typename T1, typename T2> \
2144 SACADO_INLINE_FUNCTION \
2145 bool \
2146 operator OP (const Expr<T1>& expr1, \
2147 const Expr<T2>& expr2) \
2148 { \
2149 return toBool(expr1) OP toBool(expr2); \
2150 } \
2151 \
2152 template <typename T2> \
2153 SACADO_INLINE_FUNCTION \
2154 bool \
2155 operator OP (const typename Expr<T2>::value_type& a, \
2156 const Expr<T2>& expr2) \
2157 { \
2158 return a OP toBool(expr2); \
2159 } \
2160 \
2161 template <typename T1> \
2162 SACADO_INLINE_FUNCTION \
2163 bool \
2164 operator OP (const Expr<T1>& expr1, \
2165 const typename Expr<T1>::value_type& b) \
2166 { \
2167 return toBool(expr1) OP b; \
2168 } \
2169 } \
2170 } \
2171}
2172
2175
2176#undef FAD_BOOL_MACRO
2177
2178//-------------------------- I/O Operators -----------------------
2179
2180namespace Sacado {
2181
2182 namespace Fad {
2183 namespace Exp {
2184
2185 template <typename T>
2186 std::ostream& operator << (std::ostream& os, const Expr<T>& xx) {
2187 const typename Expr<T>::derived_type& x = xx.derived();
2188 os << x.val() << " [";
2189
2190 for (int i=0; i< x.size(); i++) {
2191 os << " " << x.dx(i);
2192 }
2193
2194 os << " ]";
2195 return os;
2196 }
2197
2198 } // namespace Exp
2199 } // namespace Fad
2200
2201} // namespace Sacado
2202
2203#if defined(HAVE_SACADO_KOKKOSCORE)
2204
2205//-------------------------- Atomic Operators -----------------------
2206
2207namespace Sacado {
2208
2209 namespace Fad {
2210 namespace Exp {
2211
2212 // Overload of Kokkos::atomic_add for Fad types.
2213 template <typename S>
2215 void atomic_add(GeneralFad<S>* dst, const GeneralFad<S>& x) {
2216 using Kokkos::atomic_add;
2217
2218 const int xsz = x.size();
2219 const int sz = dst->size();
2220
2221 // We currently cannot handle resizing since that would need to be
2222 // done atomically.
2223 if (xsz > sz)
2224 Kokkos::abort(
2225 "Sacado error: Fad resize within atomic_add() not supported!");
2226
2227 if (xsz != sz && sz > 0 && xsz > 0)
2228 Kokkos::abort(
2229 "Sacado error: Fad assignment of incompatiable sizes!");
2230
2231
2232 if (sz > 0 && xsz > 0) {
2234 atomic_add(&(dst->fastAccessDx(i)), x.fastAccessDx(i));
2235 }
2237 atomic_add(&(dst->val()), x.val());
2238 }
2239
2240 } // namespace Exp
2241 } // namespace Fad
2242
2243} // namespace Sacado
2244
2245#endif // HAVE_SACADO_KOKKOSCORE
2246
2247#endif // SACADO_FAD_OPS_HPP
#define SACADO_INLINE_FUNCTION
expr expr expr bar false
expr true
#define SACADO_FAD_THREAD_SINGLE
#define SACADO_FAD_DERIV_LOOP(I, SZ)
log(expr.val())
expr expr SinOp
expr expr SinhOp
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
expr expr ACoshOp
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
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())
expr expr Log10Op
#define FAD_BOOL_MACRO(OP)
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
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
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_EXP_OP_ENABLE_EXPR_SCALAR(OP)
#define SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP)
#define SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP)
#define T1(r, f)
#define T2(r, f)
Fad specializations for Teuchos::BLAS wrappers.
Wrapper for a generic expression template.
T derived_type
Typename of derived object, returned by derived()
SACADO_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const T2 &expr2_)
SACADO_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const ConstT &c_)
SACADO_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const ConstT &c_, const T2 &expr2_)
Wrapper for a generic expression template.
std::ostream & operator<<(std::ostream &os, const Expr< T > &xx)
SACADO_INLINE_FUNCTION bool toBool(const Expr< T > &xx)
SACADO_INLINE_FUNCTION mpl::enable_if_c< IsFadExpr< T1 >::value &&IsFadExpr< T2 >::value &&ExprLevel< T1 >::value==ExprLevel< T2 >::value, IfThenElseOp< CondT, typenameExpr< T1 >::derived_type, typenameExpr< T2 >::derived_type, false, false, typenameT1::expr_spec_type > >::type if_then_else(const CondT &cond, const T1 &expr1, const 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)
Get the base Fad type from a view/expression.
Meta-function for determining nesting with an expression.
decltype(std::declval< T1 >()==std::declval< T2 >()) type
Determine whether a given type is an expression.
Is a type an expression.
static const bool value
Base template specification for IsSimdType.
static const bool value
Base template specification for Promote.
Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E >::scalar_type type
Base template specification for ScalarType.
Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E >::value_type type
Base template specification for ValueType.