Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_Fad_Exp_Ops.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Sacado Package
4 //
5 // Copyright 2006 NTESS and the Sacado contributors.
6 // SPDX-License-Identifier: LGPL-2.1-or-later
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef SACADO_FAD_EXP_OPS_HPP
11 #define SACADO_FAD_EXP_OPS_HPP
12 
13 #include <type_traits>
14 #include <ostream>
15 
19 #include "Sacado_cmath.hpp"
20 
22 
23 #if defined(HAVE_SACADO_KOKKOS)
24 #include "Kokkos_Atomic.hpp"
25 #include "impl/Kokkos_Error.hpp"
26 #endif
27 
28 #define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
29 namespace Sacado { \
30  namespace Fad { \
31  namespace Exp { \
32  \
33  template <typename T, typename ExprSpec> \
34  class OP {}; \
35  \
36  template <typename T> \
37  class OP< T,ExprSpecDefault > : \
38  public Expr< OP<T,ExprSpecDefault> > { \
39  public: \
40  \
41  typedef typename std::remove_cv<T>::type ExprT; \
42  typedef typename ExprT::value_type value_type; \
43  typedef typename ExprT::scalar_type scalar_type; \
44  \
45  typedef ExprSpecDefault expr_spec_type; \
46  \
47  SACADO_INLINE_FUNCTION \
48  explicit OP(const T& expr_) : expr(expr_) {} \
49  \
50  SACADO_INLINE_FUNCTION \
51  int size() const { return expr.size(); } \
52  \
53  SACADO_INLINE_FUNCTION \
54  bool hasFastAccess() const { \
55  return expr.hasFastAccess(); \
56  } \
57  \
58  SACADO_INLINE_FUNCTION \
59  value_type val() const { \
60  USING \
61  return VALUE; \
62  } \
63  \
64  SACADO_INLINE_FUNCTION \
65  value_type dx(int i) const { \
66  USING \
67  return DX; \
68  } \
69  \
70  SACADO_INLINE_FUNCTION \
71  value_type fastAccessDx(int i) const { \
72  USING \
73  return FASTACCESSDX; \
74  } \
75  \
76  protected: \
77  \
78  const T& expr; \
79  }; \
80  \
81  template <typename T> \
82  SACADO_INLINE_FUNCTION \
83  OP< typename Expr<T>::derived_type, \
84  typename T::expr_spec_type > \
85  OPNAME (const Expr<T>& expr) \
86  { \
87  typedef OP< typename Expr<T>::derived_type, \
88  typename T::expr_spec_type > expr_t; \
89  \
90  return expr_t(expr.derived()); \
91  } \
92  \
93  template <typename T, typename E> \
94  struct ExprLevel< OP< T,E > > { \
95  static const unsigned value = ExprLevel<T>::value; \
96  }; \
97  \
98  template <typename T, typename E> \
99  struct IsFadExpr< OP< T,E > > { \
100  static const unsigned value = true; \
101  }; \
102  \
103  } \
104  } \
105  \
106  template <typename T, typename E> \
107  struct IsExpr< Fad::Exp::OP< T,E > > { \
108  static const bool value = true; \
109  }; \
110  \
111  template <typename T, typename E> \
112  struct BaseExprType< Fad::Exp::OP< T,E > > { \
113  typedef typename BaseExprType<T>::type type; \
114  }; \
115  \
116  template <typename T, typename E> \
117  struct IsSimdType< Fad::Exp::OP< T,E > > { \
118  static const bool value = \
119  IsSimdType< typename Fad::Exp::OP< T,E >::scalar_type >::value; \
120  }; \
121  \
122  template <typename T, typename E> \
123  struct ValueType< Fad::Exp::OP< T,E > > { \
124  typedef typename Fad::Exp::OP< T,E >::value_type type; \
125  }; \
126  \
127  template <typename T, typename E> \
128  struct ScalarType< Fad::Exp::OP< T,E > > { \
129  typedef typename Fad::Exp::OP< T,E >::scalar_type type; \
130  }; \
131  \
132 }
133 
134 FAD_UNARYOP_MACRO(operator+,
135  UnaryPlusOp,
136  ;,
137  expr.val(),
138  expr.dx(i),
139  expr.fastAccessDx(i))
140 FAD_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  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()))
216 FAD_UNARYOP_MACRO(sinh,
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 // Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
270 // "simd" value types that use if_then_else(). The only reason for not using
271 // if_then_else() always is to avoid evaluating the derivative if the value is
272 // zero to avoid throwing FPEs.
273 namespace Sacado {
274  namespace Fad {
275  namespace Exp {
276 
277  template <typename T, typename ExprSpec, bool is_simd>
278  class SafeSqrtOp {};
279 
280  //
281  // Implementation for simd type using if_then_else()
282  //
283  template <typename T>
284  class SafeSqrtOp< T,ExprSpecDefault,true > :
285  public Expr< SafeSqrtOp<T,ExprSpecDefault> > {
286  public:
287 
288  typedef typename std::remove_cv<T>::type ExprT;
289  typedef typename ExprT::value_type value_type;
290  typedef typename ExprT::scalar_type scalar_type;
291 
292  typedef ExprSpecDefault expr_spec_type;
293 
295  explicit SafeSqrtOp(const T& expr_) : expr(expr_) {}
296 
298  int size() const { return expr.size(); }
299 
301  bool hasFastAccess() const {
302  return expr.hasFastAccess();
303  }
304 
306  value_type val() const {
307  using std::sqrt;
308  return sqrt(expr.val());
309  }
310 
312  value_type dx(int i) const {
313  using std::sqrt; using Sacado::if_then_else;
314  return if_then_else(
315  expr.val() == value_type(0.0), value_type(0.0),
316  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
317  }
318 
320  value_type fastAccessDx(int i) const {
321  using std::sqrt; using Sacado::if_then_else;
322  return if_then_else(
323  expr.val() == value_type(0.0), value_type(0.0),
324  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
325  }
326 
327  protected:
328 
329  const T& expr;
330  };
331 
332  //
333  // Specialization for scalar types using ternary operator
334  //
335  template <typename T>
336  class SafeSqrtOp< T,ExprSpecDefault,false > :
337  public Expr< SafeSqrtOp<T,ExprSpecDefault> > {
338  public:
339 
340  typedef typename std::remove_cv<T>::type ExprT;
341  typedef typename ExprT::value_type value_type;
342  typedef typename ExprT::scalar_type scalar_type;
343 
344  typedef ExprSpecDefault expr_spec_type;
345 
347  explicit SafeSqrtOp(const T& expr_) : expr(expr_) {}
348 
350  int size() const { return expr.size(); }
351 
353  bool hasFastAccess() const {
354  return expr.hasFastAccess();
355  }
356 
358  value_type val() const {
359  using std::sqrt;
360  return sqrt(expr.val());
361  }
362 
364  value_type dx(int i) const {
365  using std::sqrt;
366  return expr.val() == value_type(0.0) ? value_type(0.0) :
367  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
368  }
369 
371  value_type fastAccessDx(int i) const {
372  using std::sqrt;
373  return expr.val() == value_type(0.0) ? value_type(0.0) :
374  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
375  }
376 
377  protected:
378 
379  const T& expr;
380  };
381 
382  template <typename T>
384  SafeSqrtOp< typename Expr<T>::derived_type,
385  typename T::expr_spec_type >
386  safe_sqrt (const Expr<T>& expr)
387  {
388  typedef SafeSqrtOp< typename Expr<T>::derived_type,
389  typename T::expr_spec_type > expr_t;
390 
391  return expr_t(expr.derived());
392  }
393 
394  template <typename T, typename E>
395  struct ExprLevel< SafeSqrtOp< T,E > > {
396  static const unsigned value = ExprLevel<T>::value;
397  };
398 
399  template <typename T, typename E>
400  struct IsFadExpr< SafeSqrtOp< T,E > > {
401  static const unsigned value = true;
402  };
403 
404  }
405  }
406 
407  template <typename T, typename E>
408  struct IsExpr< Fad::Exp::SafeSqrtOp< T,E > > {
409  static const bool value = true;
410  };
411 
412  template <typename T, typename E>
413  struct BaseExprType< Fad::Exp::SafeSqrtOp< T,E > > {
414  typedef typename BaseExprType<T>::type type;
415  };
416 
417  template <typename T, typename E>
418  struct IsSimdType< Fad::Exp::SafeSqrtOp< T,E > > {
419  static const bool value =
420  IsSimdType< typename Fad::Exp::SafeSqrtOp< T,E >::scalar_type >::value;
421  };
422 
423  template <typename T, typename E>
424  struct ValueType< Fad::Exp::SafeSqrtOp< T,E > > {
425  typedef typename Fad::Exp::SafeSqrtOp< T,E >::value_type type;
426  };
427 
428  template <typename T, typename E>
429  struct ScalarType< Fad::Exp::SafeSqrtOp< T,E > > {
430  typedef typename Fad::Exp::SafeSqrtOp< T,E >::scalar_type type;
431  };
432 
433 }
434 
435 #undef FAD_UNARYOP_MACRO
436 
437 #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) \
438 namespace Sacado { \
439  namespace Fad { \
440  namespace Exp { \
441  \
442  template <typename T1, typename T2, \
443  bool is_const_T1, bool is_const_T2, \
444  typename ExprSpec > \
445  class OP {}; \
446  \
447  template <typename T1, typename T2> \
448  class OP< T1, T2, false, false, ExprSpecDefault > : \
449  public Expr< OP< T1, T2, false, false, ExprSpecDefault > > { \
450  public: \
451  \
452  typedef typename std::remove_cv<T1>::type ExprT1; \
453  typedef typename std::remove_cv<T2>::type ExprT2; \
454  typedef typename ExprT1::value_type value_type_1; \
455  typedef typename ExprT2::value_type value_type_2; \
456  typedef typename Sacado::Promote<value_type_1, \
457  value_type_2>::type value_type; \
458  \
459  typedef typename ExprT1::scalar_type scalar_type_1; \
460  typedef typename ExprT2::scalar_type scalar_type_2; \
461  typedef typename Sacado::Promote<scalar_type_1, \
462  scalar_type_2>::type scalar_type; \
463  \
464  typedef ExprSpecDefault expr_spec_type; \
465  \
466  SACADO_INLINE_FUNCTION \
467  OP(const T1& expr1_, const T2& expr2_) : \
468  expr1(expr1_), expr2(expr2_) {} \
469  \
470  SACADO_INLINE_FUNCTION \
471  int size() const { \
472  const int sz1 = expr1.size(), sz2 = expr2.size(); \
473  return sz1 > sz2 ? sz1 : sz2; \
474  } \
475  \
476  SACADO_INLINE_FUNCTION \
477  bool hasFastAccess() const { \
478  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
479  } \
480  \
481  SACADO_INLINE_FUNCTION \
482  value_type val() const { \
483  USING \
484  return VALUE; \
485  } \
486  \
487  SACADO_INLINE_FUNCTION \
488  value_type dx(int i) const { \
489  USING \
490  const int sz1 = expr1.size(), sz2 = expr2.size(); \
491  if (sz1 > 0 && sz2 > 0) \
492  return DX; \
493  else if (sz1 > 0) \
494  return CDX2; \
495  else \
496  return CDX1; \
497  } \
498  \
499  SACADO_INLINE_FUNCTION \
500  value_type fastAccessDx(int i) const { \
501  USING \
502  return FASTACCESSDX; \
503  } \
504  \
505  protected: \
506  \
507  const T1& expr1; \
508  const T2& expr2; \
509  \
510  }; \
511  \
512  template <typename T1, typename T2> \
513  class OP< T1, T2, false, true, ExprSpecDefault > \
514  : public Expr< OP< T1, T2, false, true, ExprSpecDefault > > { \
515  public: \
516  \
517  typedef typename std::remove_cv<T1>::type ExprT1; \
518  typedef T2 ConstT; \
519  typedef typename ExprT1::value_type value_type; \
520  typedef typename ExprT1::scalar_type scalar_type; \
521  \
522  typedef ExprSpecDefault expr_spec_type; \
523  \
524  SACADO_INLINE_FUNCTION \
525  OP(const T1& expr1_, const ConstT& c_) : \
526  expr1(expr1_), c(c_) {} \
527  \
528  SACADO_INLINE_FUNCTION \
529  int size() const { \
530  return expr1.size(); \
531  } \
532  \
533  SACADO_INLINE_FUNCTION \
534  bool hasFastAccess() const { \
535  return expr1.hasFastAccess(); \
536  } \
537  \
538  SACADO_INLINE_FUNCTION \
539  value_type val() const { \
540  USING \
541  return VAL_CONST_DX_2; \
542  } \
543  \
544  SACADO_INLINE_FUNCTION \
545  value_type dx(int i) const { \
546  USING \
547  return CONST_DX_2; \
548  } \
549  \
550  SACADO_INLINE_FUNCTION \
551  value_type fastAccessDx(int i) const { \
552  USING \
553  return CONST_FASTACCESSDX_2; \
554  } \
555  \
556  protected: \
557  \
558  const T1& expr1; \
559  const ConstT& c; \
560  }; \
561  \
562  template <typename T1, typename T2> \
563  class OP< T1, T2, true, false, ExprSpecDefault > \
564  : public Expr< OP< T1, T2, true, false, ExprSpecDefault > > { \
565  public: \
566  \
567  typedef typename std::remove_cv<T2>::type ExprT2; \
568  typedef T1 ConstT; \
569  typedef typename ExprT2::value_type value_type; \
570  typedef typename ExprT2::scalar_type scalar_type; \
571  \
572  typedef ExprSpecDefault expr_spec_type; \
573  \
574  SACADO_INLINE_FUNCTION \
575  OP(const ConstT& c_, const T2& expr2_) : \
576  c(c_), expr2(expr2_) {} \
577  \
578  SACADO_INLINE_FUNCTION \
579  int size() const { \
580  return expr2.size(); \
581  } \
582  \
583  SACADO_INLINE_FUNCTION \
584  bool hasFastAccess() const { \
585  return expr2.hasFastAccess(); \
586  } \
587  \
588  SACADO_INLINE_FUNCTION \
589  value_type val() const { \
590  USING \
591  return VAL_CONST_DX_1; \
592  } \
593  \
594  SACADO_INLINE_FUNCTION \
595  value_type dx(int i) const { \
596  USING \
597  return CONST_DX_1; \
598  } \
599  \
600  SACADO_INLINE_FUNCTION \
601  value_type fastAccessDx(int i) const { \
602  USING \
603  return CONST_FASTACCESSDX_1; \
604  } \
605  \
606  protected: \
607  \
608  const ConstT& c; \
609  const T2& expr2; \
610  }; \
611  \
612  template <typename T1, typename T2> \
613  SACADO_INLINE_FUNCTION \
614  SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP) \
615  OPNAME (const T1& expr1, const T2& expr2) \
616  { \
617  typedef OP< typename Expr<T1>::derived_type, \
618  typename Expr<T2>::derived_type, \
619  false, false, typename T1::expr_spec_type > expr_t; \
620  \
621  return expr_t(expr1.derived(), expr2.derived()); \
622  } \
623  \
624  template <typename T> \
625  SACADO_INLINE_FUNCTION \
626  OP< typename T::value_type, typename Expr<T>::derived_type, \
627  true, false, typename T::expr_spec_type > \
628  OPNAME (const typename T::value_type& c, \
629  const Expr<T>& expr) \
630  { \
631  typedef typename T::value_type ConstT; \
632  typedef OP< ConstT, typename Expr<T>::derived_type, \
633  true, false, typename T::expr_spec_type > expr_t; \
634  \
635  return expr_t(c, expr.derived()); \
636  } \
637  \
638  template <typename T> \
639  SACADO_INLINE_FUNCTION \
640  OP< typename Expr<T>::derived_type, typename T::value_type, \
641  false, true, typename T::expr_spec_type > \
642  OPNAME (const Expr<T>& expr, \
643  const typename T::value_type& c) \
644  { \
645  typedef typename T::value_type ConstT; \
646  typedef OP< typename Expr<T>::derived_type, ConstT, \
647  false, true, typename T::expr_spec_type > expr_t; \
648  \
649  return expr_t(expr.derived(), c); \
650  } \
651  \
652  template <typename T> \
653  SACADO_INLINE_FUNCTION \
654  SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP) \
655  OPNAME (const typename T::scalar_type& c, \
656  const Expr<T>& expr) \
657  { \
658  typedef typename T::scalar_type ConstT; \
659  typedef OP< ConstT, typename Expr<T>::derived_type, \
660  true, false, typename T::expr_spec_type > expr_t; \
661  \
662  return expr_t(c, expr.derived()); \
663  } \
664  \
665  template <typename T> \
666  SACADO_INLINE_FUNCTION \
667  SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP) \
668  OPNAME (const Expr<T>& expr, \
669  const typename T::scalar_type& c) \
670  { \
671  typedef typename T::scalar_type ConstT; \
672  typedef OP< typename Expr<T>::derived_type, ConstT, \
673  false, true, typename T::expr_spec_type > expr_t; \
674  \
675  return expr_t(expr.derived(), c); \
676  } \
677  \
678  template <typename T1, typename T2, bool c1, bool c2, typename E> \
679  struct ExprLevel< OP< T1, T2, c1, c2, E > > { \
680  static constexpr unsigned value_1 = ExprLevel<T1>::value; \
681  static constexpr unsigned value_2 = ExprLevel<T2>::value; \
682  static constexpr unsigned value = \
683  value_1 >= value_2 ? value_1 : value_2; \
684  }; \
685  \
686  template <typename T1, typename T2, bool c1, bool c2, typename E> \
687  struct IsFadExpr< OP< T1, T2, c1, c2, E > > { \
688  static constexpr unsigned value = true; \
689  }; \
690  \
691  } \
692  } \
693  \
694  template <typename T1, typename T2, bool c1, bool c2, typename E> \
695  struct IsExpr< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
696  static constexpr bool value = true; \
697  }; \
698  \
699  template <typename T1, typename T2, bool c1, bool c2, typename E> \
700  struct BaseExprType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
701  typedef typename BaseExprType<T1>::type base_expr_1; \
702  typedef typename BaseExprType<T2>::type base_expr_2; \
703  typedef typename Sacado::Promote<base_expr_1, \
704  base_expr_2>::type type; \
705  }; \
706  \
707  template <typename T1, typename T2, bool c1, bool c2, typename E> \
708  struct IsSimdType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
709  static const bool value = \
710  IsSimdType< typename Fad::Exp::OP< T1, T2, c1, c2, E >::value_type >::value; \
711  }; \
712  \
713  template <typename T1, typename T2, bool c1, bool c2, typename E> \
714  struct ValueType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
715  typedef typename Fad::Exp::OP< T1, T2, c1, c2, E >::value_type type;\
716  }; \
717  \
718  template <typename T1, typename T2, bool c1, bool c2, typename E> \
719  struct ScalarType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
720  typedef typename Fad::Exp::OP< T1, T2, c1, c2, E >::scalar_type type;\
721  }; \
722  \
723 }
724 
725 
726 FAD_BINARYOP_MACRO(operator+,
727  AdditionOp,
728  ;,
729  expr1.val() + expr2.val(),
730  expr1.dx(i) + expr2.dx(i),
731  expr2.dx(i),
732  expr1.dx(i),
733  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
734  c + expr2.val(),
735  expr1.val() + c,
736  expr2.dx(i),
737  expr1.dx(i),
738  expr2.fastAccessDx(i),
739  expr1.fastAccessDx(i))
740 FAD_BINARYOP_MACRO(operator-,
742  ;,
743  expr1.val() - expr2.val(),
744  expr1.dx(i) - expr2.dx(i),
745  -expr2.dx(i),
746  expr1.dx(i),
747  expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
748  c - expr2.val(),
749  expr1.val() - c,
750  -expr2.dx(i),
751  expr1.dx(i),
752  -expr2.fastAccessDx(i),
753  expr1.fastAccessDx(i))
754 FAD_BINARYOP_MACRO(operator*,
756  ;,
757  expr1.val() * expr2.val(),
758  expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
759  expr1.val()*expr2.dx(i),
760  expr1.dx(i)*expr2.val(),
761  expr1.val()*expr2.fastAccessDx(i) +
762  expr1.fastAccessDx(i)*expr2.val(),
763  c * expr2.val(),
764  expr1.val() * c,
765  c*expr2.dx(i),
766  expr1.dx(i)*c,
767  c*expr2.fastAccessDx(i),
768  expr1.fastAccessDx(i)*c)
769 FAD_BINARYOP_MACRO(operator/,
771  ;,
772  expr1.val() / expr2.val(),
773  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
774  (expr2.val()*expr2.val()),
775  -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
776  expr1.dx(i)/expr2.val(),
777  (expr1.fastAccessDx(i)*expr2.val() -
778  expr2.fastAccessDx(i)*expr1.val()) /
779  (expr2.val()*expr2.val()),
780  c / expr2.val(),
781  expr1.val() / c,
782  -expr2.dx(i)*c / (expr2.val()*expr2.val()),
783  expr1.dx(i)/c,
784  -expr2.fastAccessDx(i)*c / (expr2.val()*expr2.val()),
785  expr1.fastAccessDx(i)/c)
788  using std::atan2;,
789  atan2(expr1.val(), expr2.val()),
790  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
791  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
792  -expr1.val()*expr2.dx(i)/
793  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
794  expr2.val()*expr1.dx(i)/
795  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
796  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
797  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
798  atan2(c, expr2.val()),
799  atan2(expr1.val(), c),
800  (-c*expr2.dx(i)) / (c*c + expr2.val()*expr2.val()),
801  (c*expr1.dx(i))/ (expr1.val()*expr1.val() + c*c),
802  (-c*expr2.fastAccessDx(i))/ (c*c + expr2.val()*expr2.val()),
803  (c*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c*c))
804 // FAD_BINARYOP_MACRO(pow,
805 // PowerOp,
806 // using std::pow; using std::log; using Sacado::if_then_else;,
807 // pow(expr1.val(), expr2.val()),
808 // 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())) ),
809 // 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())) ),
810 // 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())) ),
811 // 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())) ),
812 // pow(c, expr2.val()),
813 // pow(expr1.val(), c),
814 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) ),
815 // 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)) ),
816 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) ),
817 // 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))) )
820  using Sacado::if_then_else;,
821  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
822  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
823  if_then_else( expr1.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
824  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), value_type(0.0) ),
825  if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
826  if_then_else( c >= expr2.val(), value_type(c), expr2.val() ),
827  if_then_else( expr1.val() >= c, expr1.val(), value_type(c) ),
828  if_then_else( c >= expr2.val(), value_type(0.0), expr2.dx(i) ),
829  if_then_else( expr1.val() >= c, expr1.dx(i), value_type(0.0) ),
830  if_then_else( c >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
831  if_then_else( expr1.val() >= c, expr1.fastAccessDx(i), value_type(0.0) ) )
834  using Sacado::if_then_else;,
835  if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
836  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
837  if_then_else( expr1.val() <= expr2.val(), value_type(0.0), expr2.dx(i) ),
838  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), value_type(0.0) ),
839  if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
840  if_then_else( c <= expr2.val(), value_type(c), expr2.val() ),
841  if_then_else( expr1.val() <= c, expr1.val(), value_type(c) ),
842  if_then_else( c <= expr2.val(), value_type(0), expr2.dx(i) ),
843  if_then_else( expr1.val() <= c, expr1.dx(i), value_type(0) ),
844  if_then_else( c <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
845  if_then_else( expr1.val() <= c, expr1.fastAccessDx(i), value_type(0) ) )
846 
847 
848 // Special handling for std::pow() to provide specializations of PowerOp for
849 // "simd" value types that use if_then_else(). The only reason for not using
850 // if_then_else() always is to avoid evaluating the derivative if the value is
851 // zero to avoid throwing FPEs.
852 namespace Sacado {
853  namespace Fad {
854  namespace Exp {
855 
856  template <typename T1, typename T2,
857  bool is_const_T1, bool is_const_T2,
858  typename ExprSpec, typename Impl >
859  class PowerOp {};
860 
861  //
862  // Implementation for simd type using if_then_else()
863  //
864  template <typename T1, typename T2>
865  class PowerOp< T1, T2, false, false, ExprSpecDefault,
866  PowerImpl::Simd > :
867  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
868  PowerImpl::Simd > > {
869  public:
870 
871  typedef typename std::remove_cv<T1>::type ExprT1;
872  typedef typename std::remove_cv<T2>::type ExprT2;
873  typedef typename ExprT1::value_type value_type_1;
874  typedef typename ExprT2::value_type value_type_2;
875  typedef typename Sacado::Promote<value_type_1,
876  value_type_2>::type value_type;
877 
878  typedef typename ExprT1::scalar_type scalar_type_1;
879  typedef typename ExprT2::scalar_type scalar_type_2;
880  typedef typename Sacado::Promote<scalar_type_1,
881  scalar_type_2>::type scalar_type;
882 
883  typedef ExprSpecDefault expr_spec_type;
884 
886  PowerOp(const T1& expr1_, const T2& expr2_) :
887  expr1(expr1_), expr2(expr2_) {}
888 
890  int size() const {
891  const int sz1 = expr1.size(), sz2 = expr2.size();
892  return sz1 > sz2 ? sz1 : sz2;
893  }
894 
896  bool hasFastAccess() const {
897  return expr1.hasFastAccess() && expr2.hasFastAccess();
898  }
899 
901  value_type val() const {
902  using std::pow;
903  return pow(expr1.val(), expr2.val());
904  }
905 
907  value_type dx(int i) const {
908  using std::pow; using std::log; using Sacado::if_then_else;
909  const int sz1 = expr1.size(), sz2 = expr2.size();
910  if (sz1 > 0 && sz2 > 0)
911  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())) );
912  else if (sz1 > 0)
913  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
914  // It seems less accurate and caused convergence problems in some codes
915  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())) ));
916  else
917  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())) );
918  }
919 
921  value_type fastAccessDx(int i) const {
922  using std::pow; using std::log; using Sacado::if_then_else;
923  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())));
924  }
925 
926  protected:
927 
928  const T1& expr1;
929  const T2& expr2;
930 
931  };
932 
933  template <typename T1, typename T2>
934  class PowerOp< T1, T2, false, true, ExprSpecDefault,
935  PowerImpl::Simd > :
936  public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
937  PowerImpl::Simd > > {
938  public:
939 
940  typedef typename std::remove_cv<T1>::type ExprT1;
941  typedef T2 ConstT;
942  typedef typename ExprT1::value_type value_type;
943  typedef typename ExprT1::scalar_type scalar_type;
944 
945  typedef ExprSpecDefault expr_spec_type;
946 
948  PowerOp(const T1& expr1_, const ConstT& c_) :
949  expr1(expr1_), c(c_) {}
950 
952  int size() const {
953  return expr1.size();
954  }
955 
957  bool hasFastAccess() const {
958  return expr1.hasFastAccess();
959  }
960 
962  value_type val() const {
963  using std::pow;
964  return pow(expr1.val(), c);
965  }
966 
968  value_type dx(int i) const {
969  using std::pow; using Sacado::if_then_else;
970  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
971  // It seems less accurate and caused convergence problems in some codes
972  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)) ));
973  }
974 
976  value_type fastAccessDx(int i) const {
977  using std::pow; using Sacado::if_then_else;
978  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
979  // It seems less accurate and caused convergence problems in some codes
980  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)) ));
981  }
982 
983  protected:
984 
985  const T1& expr1;
986  const ConstT& c;
987  };
988 
989  template <typename T1, typename T2>
990  class PowerOp< T1, T2, true, false, ExprSpecDefault,
991  PowerImpl::Simd > :
992  public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
993  PowerImpl::Simd > > {
994  public:
995 
996  typedef typename std::remove_cv<T2>::type ExprT2;
997  typedef T1 ConstT;
998  typedef typename ExprT2::value_type value_type;
999  typedef typename ExprT2::scalar_type scalar_type;
1000 
1001  typedef ExprSpecDefault expr_spec_type;
1002 
1004  PowerOp(const ConstT& c_, const T2& expr2_) :
1005  c(c_), expr2(expr2_) {}
1006 
1008  int size() const {
1009  return expr2.size();
1010  }
1011 
1013  bool hasFastAccess() const {
1014  return expr2.hasFastAccess();
1015  }
1016 
1018  value_type val() const {
1019  using std::pow;
1020  return pow(c, expr2.val());
1021  }
1022 
1024  value_type dx(int i) const {
1025  using std::pow; using std::log; using Sacado::if_then_else;
1026  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) );
1027  }
1028 
1030  value_type fastAccessDx(int i) const {
1031  using std::pow; using std::log; using Sacado::if_then_else;
1032  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) );
1033  }
1034 
1035  protected:
1036 
1037  const ConstT& c;
1038  const T2& expr2;
1039  };
1040 
1041  //
1042  // Specialization for scalar types using ternary operator
1043  //
1044 
1045  template <typename T1, typename T2>
1046  class PowerOp< T1, T2, false, false, ExprSpecDefault,
1047  PowerImpl::Scalar > :
1048  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
1049  PowerImpl::Scalar > > {
1050  public:
1051 
1052  typedef typename std::remove_cv<T1>::type ExprT1;
1053  typedef typename std::remove_cv<T2>::type ExprT2;
1054  typedef typename ExprT1::value_type value_type_1;
1055  typedef typename ExprT2::value_type value_type_2;
1056  typedef typename Sacado::Promote<value_type_1,
1057  value_type_2>::type value_type;
1058 
1059  typedef typename ExprT1::scalar_type scalar_type_1;
1060  typedef typename ExprT2::scalar_type scalar_type_2;
1061  typedef typename Sacado::Promote<scalar_type_1,
1062  scalar_type_2>::type scalar_type;
1063 
1064  typedef ExprSpecDefault expr_spec_type;
1065 
1067  PowerOp(const T1& expr1_, const T2& expr2_) :
1068  expr1(expr1_), expr2(expr2_) {}
1069 
1071  int size() const {
1072  const int sz1 = expr1.size(), sz2 = expr2.size();
1073  return sz1 > sz2 ? sz1 : sz2;
1074  }
1075 
1077  bool hasFastAccess() const {
1078  return expr1.hasFastAccess() && expr2.hasFastAccess();
1079  }
1080 
1082  value_type val() const {
1083  using std::pow;
1084  return pow(expr1.val(), expr2.val());
1085  }
1086 
1088  value_type dx(int i) const {
1089  using std::pow; using std::log;
1090  const int sz1 = expr1.size(), sz2 = expr2.size();
1091  if (sz1 > 0 && sz2 > 0)
1092  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()));
1093  else if (sz1 > 0)
1094  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1095  // It seems less accurate and caused convergence problems in some codes
1096  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()));
1097  else
1098  return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
1099  }
1100 
1102  value_type fastAccessDx(int i) const {
1103  using std::pow; using std::log;
1104  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()));
1105  }
1106 
1107  protected:
1108 
1109  const T1& expr1;
1110  const T2& expr2;
1111 
1112  };
1113 
1114  template <typename T1, typename T2>
1115  class PowerOp< T1, T2, false, true, ExprSpecDefault,
1116  PowerImpl::Scalar > :
1117  public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
1118  PowerImpl::Scalar > > {
1119  public:
1120 
1121  typedef typename std::remove_cv<T1>::type ExprT1;
1122  typedef T2 ConstT;
1123  typedef typename ExprT1::value_type value_type;
1124  typedef typename ExprT1::scalar_type scalar_type;
1125 
1126  typedef ExprSpecDefault expr_spec_type;
1127 
1129  PowerOp(const T1& expr1_, const ConstT& c_) :
1130  expr1(expr1_), c(c_) {}
1131 
1133  int size() const {
1134  return expr1.size();
1135  }
1136 
1138  bool hasFastAccess() const {
1139  return expr1.hasFastAccess();
1140  }
1141 
1143  value_type val() const {
1144  using std::pow;
1145  return pow(expr1.val(), c);
1146  }
1147 
1149  value_type dx(int i) const {
1150  using std::pow;
1151  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1152  // It seems less accurate and caused convergence problems in some codes
1153  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));
1154  }
1155 
1157  value_type fastAccessDx(int i) const {
1158  using std::pow;
1159  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1160  // It seems less accurate and caused convergence problems in some codes
1161  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));
1162  }
1163 
1164  protected:
1165 
1166  const T1& expr1;
1167  const ConstT& c;
1168  };
1169 
1170  template <typename T1, typename T2>
1171  class PowerOp< T1, T2, true, false, ExprSpecDefault,
1172  PowerImpl::Scalar > :
1173  public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
1174  PowerImpl::Scalar > > {
1175  public:
1176 
1177  typedef typename std::remove_cv<T2>::type ExprT2;
1178  typedef T1 ConstT;
1179  typedef typename ExprT2::value_type value_type;
1180  typedef typename ExprT2::scalar_type scalar_type;
1181 
1182  typedef ExprSpecDefault expr_spec_type;
1183 
1185  PowerOp(const ConstT& c_, const T2& expr2_) :
1186  c(c_), expr2(expr2_) {}
1187 
1189  int size() const {
1190  return expr2.size();
1191  }
1192 
1194  bool hasFastAccess() const {
1195  return expr2.hasFastAccess();
1196  }
1197 
1199  value_type val() const {
1200  using std::pow;
1201  return pow(c, expr2.val());
1202  }
1203 
1205  value_type dx(int i) const {
1206  using std::pow; using std::log;
1207  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c)*pow(c,expr2.val()));
1208  }
1209 
1211  value_type fastAccessDx(int i) const {
1212  using std::pow; using std::log;
1213  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val()));
1214  }
1215 
1216  protected:
1217 
1218  const ConstT& c;
1219  const T2& expr2;
1220  };
1221 
1222  //
1223  // Specialization for nested derivatives. This version does not use
1224  // if_then_else/ternary-operator on the base so that nested derivatives
1225  // are correct.
1226  //
1227  template <typename T1, typename T2>
1228  class PowerOp< T1, T2, false, false, ExprSpecDefault,
1229  PowerImpl::Nested > :
1230  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
1231  PowerImpl::Nested > > {
1232  public:
1233 
1234  typedef typename std::remove_cv<T1>::type ExprT1;
1235  typedef typename std::remove_cv<T2>::type ExprT2;
1236  typedef typename ExprT1::value_type value_type_1;
1237  typedef typename ExprT2::value_type value_type_2;
1238  typedef typename Sacado::Promote<value_type_1,
1239  value_type_2>::type value_type;
1240 
1241  typedef typename ExprT1::scalar_type scalar_type_1;
1242  typedef typename ExprT2::scalar_type scalar_type_2;
1243  typedef typename Sacado::Promote<scalar_type_1,
1244  scalar_type_2>::type scalar_type;
1245 
1246  typedef ExprSpecDefault expr_spec_type;
1247 
1249  PowerOp(const T1& expr1_, const T2& expr2_) :
1250  expr1(expr1_), expr2(expr2_) {}
1251 
1253  int size() const {
1254  const int sz1 = expr1.size(), sz2 = expr2.size();
1255  return sz1 > sz2 ? sz1 : sz2;
1256  }
1257 
1259  bool hasFastAccess() const {
1260  return expr1.hasFastAccess() && expr2.hasFastAccess();
1261  }
1262 
1264  value_type val() const {
1265  using std::pow;
1266  return pow(expr1.val(), expr2.val());
1267  }
1268 
1270  value_type dx(int i) const {
1271  using std::pow; using std::log;
1272  const int sz1 = expr1.size(), sz2 = expr2.size();
1273  if (sz1 > 0 && sz2 > 0)
1274  return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1275  else if (sz1 > 0)
1276  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)));
1277  else
1278  return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1279  }
1280 
1282  value_type fastAccessDx(int i) const {
1283  using std::pow; using std::log;
1284  return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1285  }
1286 
1287  protected:
1288 
1289  const T1& expr1;
1290  const T2& expr2;
1291 
1292  };
1293 
1294  template <typename T1, typename T2>
1295  class PowerOp< T1, T2, false, true, ExprSpecDefault,
1296  PowerImpl::Nested > :
1297  public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
1298  PowerImpl::Nested > > {
1299  public:
1300 
1301  typedef typename std::remove_cv<T1>::type ExprT1;
1302  typedef T2 ConstT;
1303  typedef typename ExprT1::value_type value_type;
1304  typedef typename ExprT1::scalar_type scalar_type;
1305 
1306  typedef ExprSpecDefault expr_spec_type;
1307 
1309  PowerOp(const T1& expr1_, const ConstT& c_) :
1310  expr1(expr1_), c(c_) {}
1311 
1313  int size() const {
1314  return expr1.size();
1315  }
1316 
1318  bool hasFastAccess() const {
1319  return expr1.hasFastAccess();
1320  }
1321 
1323  value_type val() const {
1324  using std::pow;
1325  return pow(expr1.val(), c);
1326  }
1327 
1329  value_type dx(int i) const {
1330  using std::pow;
1331  return c == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)*pow(expr1.val(),c-scalar_type(1.0)));
1332  }
1333 
1335  value_type fastAccessDx(int i) const {
1336  using std::pow;
1337  return c == scalar_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)*pow(expr1.val(),c-scalar_type(1.0)));
1338  }
1339 
1340  protected:
1341 
1342  const T1& expr1;
1343  const ConstT& c;
1344  };
1345 
1346  template <typename T1, typename T2>
1347  class PowerOp<T1, T2, true, false, ExprSpecDefault,
1348  PowerImpl::Nested > :
1349  public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
1350  PowerImpl::Nested > > {
1351  public:
1352 
1353  typedef typename std::remove_cv<T2>::type ExprT2;
1354  typedef T1 ConstT;
1355  typedef typename ExprT2::value_type value_type;
1356  typedef typename ExprT2::scalar_type scalar_type;
1357 
1358  typedef ExprSpecDefault expr_spec_type;
1359 
1361  PowerOp(const ConstT& c_, const T2& expr2_) :
1362  c(c_), expr2(expr2_) {}
1363 
1365  int size() const {
1366  return expr2.size();
1367  }
1368 
1370  bool hasFastAccess() const {
1371  return expr2.hasFastAccess();
1372  }
1373 
1375  value_type val() const {
1376  using std::pow;
1377  return pow(c, expr2.val());
1378  }
1379 
1381  value_type dx(int i) const {
1382  using std::pow; using std::log;
1383  return expr2.dx(i)*log(c)*pow(c,expr2.val());
1384  }
1385 
1387  value_type fastAccessDx(int i) const {
1388  using std::pow; using std::log;
1389  return expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val());
1390  }
1391 
1392  protected:
1393 
1394  const ConstT& c;
1395  const T2& expr2;
1396  };
1397 
1398  //
1399  // Specialization for nested derivatives. This version does not use
1400  // if_then_else/ternary-operator on the base so that nested derivatives
1401  // are correct.
1402  //
1403  template <typename T1, typename T2>
1404  class PowerOp< T1, T2, false, false, ExprSpecDefault,
1405  PowerImpl::NestedSimd > :
1406  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault,
1407  PowerImpl::NestedSimd > > {
1408  public:
1409 
1410  typedef typename std::remove_cv<T1>::type ExprT1;
1411  typedef typename std::remove_cv<T2>::type ExprT2;
1412  typedef typename ExprT1::value_type value_type_1;
1413  typedef typename ExprT2::value_type value_type_2;
1414  typedef typename Sacado::Promote<value_type_1,
1415  value_type_2>::type value_type;
1416 
1417  typedef typename ExprT1::scalar_type scalar_type_1;
1418  typedef typename ExprT2::scalar_type scalar_type_2;
1419  typedef typename Sacado::Promote<scalar_type_1,
1420  scalar_type_2>::type scalar_type;
1421 
1422  typedef ExprSpecDefault expr_spec_type;
1423 
1425  PowerOp(const T1& expr1_, const T2& expr2_) :
1426  expr1(expr1_), expr2(expr2_) {}
1427 
1429  int size() const {
1430  const int sz1 = expr1.size(), sz2 = expr2.size();
1431  return sz1 > sz2 ? sz1 : sz2;
1432  }
1433 
1435  bool hasFastAccess() const {
1436  return expr1.hasFastAccess() && expr2.hasFastAccess();
1437  }
1438 
1440  value_type val() const {
1441  using std::pow;
1442  return pow(expr1.val(), expr2.val());
1443  }
1444 
1446  value_type dx(int i) const {
1447  using std::pow; using std::log; using Sacado::if_then_else;
1448  const int sz1 = expr1.size(), sz2 = expr2.size();
1449  if (sz1 > 0 && sz2 > 0)
1450  return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1451  else if (sz1 > 0)
1452  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))));
1453  else
1454  return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1455  }
1456 
1458  value_type fastAccessDx(int i) const {
1459  using std::pow; using std::log; using Sacado::if_then_else;
1460  return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1461  }
1462 
1463  protected:
1464 
1465  const T1& expr1;
1466  const T2& expr2;
1467 
1468  };
1469 
1470  template <typename T1, typename T2>
1471  class PowerOp< T1, T2, false, true, ExprSpecDefault,
1472  PowerImpl::NestedSimd > :
1473  public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault,
1474  PowerImpl::NestedSimd > > {
1475  public:
1476 
1477  typedef typename std::remove_cv<T1>::type ExprT1;
1478  typedef T2 ConstT;
1479  typedef typename ExprT1::value_type value_type;
1480  typedef typename ExprT1::scalar_type scalar_type;
1481 
1482  typedef ExprSpecDefault expr_spec_type;
1483 
1485  PowerOp(const T1& expr1_, const ConstT& c_) :
1486  expr1(expr1_), c(c_) {}
1487 
1489  int size() const {
1490  return expr1.size();
1491  }
1492 
1494  bool hasFastAccess() const {
1495  return expr1.hasFastAccess();
1496  }
1497 
1499  value_type val() const {
1500  using std::pow;
1501  return pow(expr1.val(), c);
1502  }
1503 
1505  value_type dx(int i) const {
1506  using std::pow; using Sacado::if_then_else;
1507  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))));
1508  }
1509 
1511  value_type fastAccessDx(int i) const {
1512  using std::pow; using Sacado::if_then_else;
1513  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))));
1514  }
1515 
1516  protected:
1517 
1518  const T1& expr1;
1519  const ConstT& c;
1520  };
1521 
1522  template <typename T1, typename T2>
1523  class PowerOp<T1, T2, true, false, ExprSpecDefault,
1524  PowerImpl::NestedSimd > :
1525  public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault,
1526  PowerImpl::NestedSimd > > {
1527  public:
1528 
1529  typedef typename std::remove_cv<T2>::type ExprT2;
1530  typedef T1 ConstT;
1531  typedef typename ExprT2::value_type value_type;
1532  typedef typename ExprT2::scalar_type scalar_type;
1533 
1534  typedef ExprSpecDefault expr_spec_type;
1535 
1537  PowerOp(const ConstT& c_, const T2& expr2_) :
1538  c(c_), expr2(expr2_) {}
1539 
1541  int size() const {
1542  return expr2.size();
1543  }
1544 
1546  bool hasFastAccess() const {
1547  return expr2.hasFastAccess();
1548  }
1549 
1551  value_type val() const {
1552  using std::pow;
1553  return pow(c, expr2.val());
1554  }
1555 
1557  value_type dx(int i) const {
1558  using std::pow; using std::log;
1559  return expr2.dx(i)*log(c)*pow(c,expr2.val());
1560  }
1561 
1563  value_type fastAccessDx(int i) const {
1564  using std::pow; using std::log;
1565  return expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val());
1566  }
1567 
1568  protected:
1569 
1570  const ConstT& c;
1571  const T2& expr2;
1572  };
1573 
1574  template <typename T1, typename T2>
1577  pow (const T1& expr1, const T2& expr2)
1578  {
1579  typedef PowerOp< typename Expr<T1>::derived_type,
1580  typename Expr<T2>::derived_type,
1581  false, false, typename T1::expr_spec_type > expr_t;
1582 
1583  return expr_t(expr1.derived(), expr2.derived());
1584  }
1585 
1586  template <typename T>
1588  PowerOp< typename T::value_type, typename Expr<T>::derived_type,
1589  true, false, typename T::expr_spec_type >
1590  pow (const typename T::value_type& c,
1591  const Expr<T>& expr)
1592  {
1593  typedef typename T::value_type ConstT;
1594  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1595  true, false, typename T::expr_spec_type > expr_t;
1596 
1597  return expr_t(c, expr.derived());
1598  }
1599 
1600  template <typename T>
1602  PowerOp< typename Expr<T>::derived_type, typename T::value_type,
1603  false, true, typename T::expr_spec_type >
1604  pow (const Expr<T>& expr,
1605  const typename T::value_type& c)
1606  {
1607  typedef typename T::value_type ConstT;
1608  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1609  false, true, typename T::expr_spec_type > expr_t;
1610 
1611  return expr_t(expr.derived(), c);
1612  }
1613 
1614  template <typename T>
1617  pow (const typename T::scalar_type& c,
1618  const Expr<T>& expr)
1619  {
1620  typedef typename T::scalar_type ConstT;
1621  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1622  true, false, typename T::expr_spec_type > expr_t;
1623 
1624  return expr_t(c, expr.derived());
1625  }
1626 
1627  template <typename T>
1630  pow (const Expr<T>& expr,
1631  const typename T::scalar_type& c)
1632  {
1633  typedef typename T::scalar_type ConstT;
1634  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1635  false, true, typename T::expr_spec_type > expr_t;
1636 
1637  return expr_t(expr.derived(), c);
1638  }
1639 
1640  template <typename T1, typename T2, bool c1, bool c2, typename E>
1641  struct ExprLevel< PowerOp< T1, T2, c1, c2, E > > {
1642  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1643  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1644  static constexpr unsigned value =
1645  value_1 >= value_2 ? value_1 : value_2;
1646  };
1647 
1648  template <typename T1, typename T2, bool c1, bool c2, typename E>
1649  struct IsFadExpr< PowerOp< T1, T2, c1, c2, E > > {
1650  static constexpr unsigned value = true;
1651  };
1652 
1653  }
1654  }
1655 
1656  template <typename T1, typename T2, bool c1, bool c2, typename E>
1657  struct IsExpr< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1658  static constexpr bool value = true;
1659  };
1660 
1661  template <typename T1, typename T2, bool c1, bool c2, typename E>
1662  struct BaseExprType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1663  typedef typename BaseExprType<T1>::type base_expr_1;
1664  typedef typename BaseExprType<T2>::type base_expr_2;
1665  typedef typename Sacado::Promote<base_expr_1,
1666  base_expr_2>::type type;
1667  };
1668 
1669  template <typename T1, typename T2, bool c1, bool c2, typename E>
1670  struct IsSimdType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1671  static const bool value =
1672  IsSimdType< typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::value_type >::value;
1673  };
1674 
1675  template <typename T1, typename T2, bool c1, bool c2, typename E>
1676  struct ValueType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1677  typedef typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::value_type type;
1678  };
1679 
1680  template <typename T1, typename T2, bool c1, bool c2, typename E>
1681  struct ScalarType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1682  typedef typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::scalar_type type;
1683  };
1684 
1685 }
1686 
1687 //--------------------------if_then_else operator -----------------------
1688 // Can't use the above macros because it is a ternary operator (sort of).
1689 // Also, relies on C++11
1690 
1691 namespace Sacado {
1692  namespace Fad {
1693  namespace Exp {
1694 
1695  template <typename CondT, typename T1, typename T2,
1696  bool is_const_T1, bool is_const_T2,
1697  typename ExprSpec>
1698  class IfThenElseOp {};
1699 
1700  template <typename CondT, typename T1, typename T2>
1702  public Expr< IfThenElseOp< CondT, T1, T2, false, false, ExprSpecDefault > > {
1703 
1704  public:
1705 
1706  typedef typename std::remove_cv<T1>::type ExprT1;
1707  typedef typename std::remove_cv<T2>::type ExprT2;
1708  typedef typename ExprT1::value_type value_type_1;
1709  typedef typename ExprT2::value_type value_type_2;
1710  typedef typename Sacado::Promote<value_type_1,
1712 
1713  typedef typename ExprT1::scalar_type scalar_type_1;
1714  typedef typename ExprT2::scalar_type scalar_type_2;
1715  typedef typename Sacado::Promote<scalar_type_1,
1717 
1719 
1721  IfThenElseOp(const CondT& cond_, const T1& expr1_, const T2& expr2_) :
1722  cond(cond_), expr1(expr1_), expr2(expr2_) {}
1723 
1725  int size() const {
1726  int sz1 = expr1.size(), sz2 = expr2.size();
1727  return sz1 > sz2 ? sz1 : sz2;
1728  }
1729 
1731  bool hasFastAccess() const {
1732  return expr1.hasFastAccess() && expr2.hasFastAccess();
1733  }
1734 
1736  value_type val() const {
1737  using Sacado::if_then_else;
1738  return if_then_else( cond, expr1.val(), expr2.val() );
1739  }
1740 
1742  value_type dx(int i) const {
1743  using Sacado::if_then_else;
1744  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
1745  }
1746 
1749  using Sacado::if_then_else;
1750  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
1751  }
1752 
1753  protected:
1754 
1755  const CondT& cond;
1756  const T1& expr1;
1757  const T2& expr2;
1758 
1759  };
1760 
1761  template <typename CondT, typename T1, typename T2>
1762  class IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > :
1763  public Expr< IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > > {
1764 
1765  public:
1766 
1767  typedef typename std::remove_cv<T1>::type ExprT1;
1768  typedef T2 ConstT;
1769  typedef typename ExprT1::value_type value_type;
1770  typedef typename ExprT1::scalar_type scalar_type;
1771 
1773 
1775  IfThenElseOp(const CondT& cond_, const T1& expr1_, const ConstT& c_) :
1776  cond(cond_), expr1(expr1_), c(c_) {}
1777 
1779  int size() const {
1780  return expr1.size();
1781  }
1782 
1784  bool hasFastAccess() const {
1785  return expr1.hasFastAccess();
1786  }
1787 
1789  value_type val() const {
1790  using Sacado::if_then_else;
1791  return if_then_else( cond, expr1.val(), c );
1792  }
1793 
1795  value_type dx(int i) const {
1796  using Sacado::if_then_else;
1797  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
1798  }
1799 
1802  using Sacado::if_then_else;
1803  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
1804  }
1805 
1806  protected:
1807 
1808  const CondT& cond;
1809  const T1& expr1;
1810  const ConstT& c;
1811  };
1812 
1813  template <typename CondT, typename T1, typename T2>
1814  class IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > :
1815  public Expr< IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > > {
1816 
1817  public:
1818 
1819  typedef typename std::remove_cv<T2>::type ExprT2;
1820  typedef T1 ConstT;
1821  typedef typename ExprT2::value_type value_type;
1822  typedef typename ExprT2::scalar_type scalar_type;
1823 
1825 
1827  IfThenElseOp(const CondT& cond_, const ConstT& c_, const T2& expr2_) :
1828  cond(cond_), c(c_), expr2(expr2_) {}
1829 
1831  int size() const {
1832  return expr2.size();
1833  }
1834 
1836  bool hasFastAccess() const {
1837  return expr2.hasFastAccess();
1838  }
1839 
1841  value_type val() const {
1842  using Sacado::if_then_else;
1843  return if_then_else( cond, c, expr2.val() );
1844  }
1845 
1847  value_type dx(int i) const {
1848  using Sacado::if_then_else;
1849  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
1850  }
1851 
1854  using Sacado::if_then_else;
1855  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
1856  }
1857 
1858  protected:
1859 
1860  const CondT& cond;
1861  const ConstT& c;
1862  const T2& expr2;
1863  };
1864 
1865  template <typename CondT, typename T1, typename T2>
1869  IfThenElseOp< CondT,
1870  typename Expr<T1>::derived_type,
1871  typename Expr<T2>::derived_type,
1872  false, false,
1873  typename T1::expr_spec_type >
1874  >::type
1875  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
1876  {
1878  typename Expr<T2>::derived_type,
1879  false, false, typename T1::expr_spec_type > expr_t;
1880 
1881  return expr_t(cond, expr1.derived(), expr2.derived());
1882  }
1883 
1884  template <typename CondT, typename T>
1887  true, false, typename T::expr_spec_type >
1888  if_then_else (const CondT& cond, const typename T::value_type& c,
1889  const Expr<T>& expr)
1890  {
1891  typedef typename T::value_type ConstT;
1893  true, false, typename T::expr_spec_type > expr_t;
1894 
1895  return expr_t(cond, c, expr.derived());
1896  }
1897 
1898  template <typename CondT, typename T>
1901  false, true, typename T::expr_spec_type >
1902  if_then_else (const CondT& cond, const Expr<T>& expr,
1903  const typename T::value_type& c)
1904  {
1905  typedef typename T::value_type ConstT;
1907  false, true, typename T::expr_spec_type > expr_t;
1908 
1909  return expr_t(cond, expr.derived(), c);
1910  }
1911 
1912  template <typename CondT, typename T>
1914  typename mpl::disable_if<
1915  std::is_same< typename T::value_type,
1916  typename T::scalar_type >,
1917  IfThenElseOp< CondT, typename T::scalar_type,
1918  typename Expr<T>::derived_type,
1919  true, false, typename T::expr_spec_type >
1920  >::type
1921  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
1922  const Expr<T>& expr)
1923  {
1924  typedef typename T::scalar_type ConstT;
1926  true, false, typename T::expr_spec_type > expr_t;
1927 
1928  return expr_t(cond, c, expr.derived());
1929  }
1930 
1931  template <typename CondT, typename T>
1933  typename mpl::disable_if<
1934  std::is_same< typename T::value_type,
1935  typename T::scalar_type >,
1937  typename T::scalar_type,
1938  false, true, typename T::expr_spec_type >
1939  >::type
1940  if_then_else (const CondT& cond, const Expr<T>& expr,
1941  const typename Expr<T>::scalar_type& c)
1942  {
1943  typedef typename T::scalar_type ConstT;
1945  false, true, typename T::expr_spec_type > expr_t;
1946 
1947  return expr_t(cond, expr.derived(), c);
1948  }
1949 
1950  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1951  typename E>
1952  struct ExprLevel< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1953  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1954  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1955  static constexpr unsigned value =
1956  value_1 >= value_2 ? value_1 : value_2;
1957  };
1958 
1959  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1960  typename E>
1961  struct IsFadExpr< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1962  static constexpr unsigned value = true;
1963  };
1964 
1965  }
1966  }
1967 
1968  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1969  typename E>
1970  struct IsExpr< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1971  static constexpr bool value = true;
1972  };
1973 
1974  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1975  typename E>
1976  struct BaseExprType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1979  typedef typename Sacado::Promote<base_expr_1,
1981  };
1982 
1983  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1984  typename E>
1985  struct IsSimdType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1986  static const bool value =
1988  };
1989 
1990  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1991  typename E>
1992  struct ValueType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1994  };
1995 
1996  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1997  typename E>
1998  struct ScalarType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
2000  };
2001 }
2002 
2003 #undef FAD_BINARYOP_MACRO
2004 
2005 //-------------------------- Relational Operators -----------------------
2006 
2007 namespace Sacado {
2008  namespace Fad {
2009  namespace Exp {
2010  namespace Impl {
2011  // Helper trait to determine return type of logical comparison operations
2012  // (==, !=, ...), usually bool but maybe something else for SIMD types.
2013  // Need to use SFINAE to restrict to types that define == operator in the
2014  // conditional overloads below, otherwise instantiating ConditionaReturnType
2015  // may fail during overload resolution.
2016  template <typename T1, typename T2 = T1,
2019 
2020  template <typename T1, typename T2>
2022  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
2023  };
2024  }
2025  }
2026  }
2027 }
2028 
2029 #define FAD_RELOP_MACRO(OP) \
2030 namespace Sacado { \
2031  namespace Fad { \
2032  namespace Exp { \
2033  template <typename T1, typename T2> \
2034  SACADO_INLINE_FUNCTION \
2035  typename mpl::enable_if_c< \
2036  IsFadExpr<T1>::value && IsFadExpr<T2>::value && \
2037  ExprLevel<T1>::value == ExprLevel<T2>::value, \
2038  typename Impl::ConditionalReturnType<typename T1::value_type, \
2039  typename T2::value_type>::type \
2040  >::type \
2041  operator OP (const T1& expr1, const T2& expr2) \
2042  { \
2043  return expr1.derived().val() OP expr2.derived().val(); \
2044  } \
2045  \
2046  template <typename T2> \
2047  SACADO_INLINE_FUNCTION \
2048  typename Impl::ConditionalReturnType<typename T2::value_type>::type \
2049  operator OP (const typename T2::value_type& a, \
2050  const Expr<T2>& expr2) \
2051  { \
2052  return a OP expr2.derived().val(); \
2053  } \
2054  \
2055  template <typename T1> \
2056  SACADO_INLINE_FUNCTION \
2057  typename Impl::ConditionalReturnType<typename T1::value_type>::type \
2058  operator OP (const Expr<T1>& expr1, \
2059  const typename T1::value_type& b) \
2060  { \
2061  return expr1.derived().val() OP b; \
2062  } \
2063  } \
2064  } \
2065 } \
2066 
2067 FAD_RELOP_MACRO(==)
2068 FAD_RELOP_MACRO(!=)
2069 FAD_RELOP_MACRO(<)
2070 FAD_RELOP_MACRO(>)
2071 FAD_RELOP_MACRO(<=)
2072 FAD_RELOP_MACRO(>=)
2073 FAD_RELOP_MACRO(<<=)
2074 FAD_RELOP_MACRO(>>=)
2075 FAD_RELOP_MACRO(&)
2076 FAD_RELOP_MACRO(|)
2077 
2078 #undef FAD_RELOP_MACRO
2079 
2080 namespace Sacado {
2081 
2082  namespace Fad {
2083  namespace Exp {
2084 
2085  template <typename ExprT>
2087  bool operator ! (const Expr<ExprT>& expr)
2088  {
2089  return ! expr.derived().val();
2090  }
2091 
2092  } // namespace Exp
2093  } // namespace Fad
2094 
2095 } // namespace Sacado
2096 
2097 //-------------------------- Boolean Operators -----------------------
2098 namespace Sacado {
2099 
2100  namespace Fad {
2101  namespace Exp {
2102 
2103  template <typename T>
2105  bool toBool(const Expr<T>& xx) {
2106  const typename Expr<T>::derived_type& x = xx.derived();
2107  bool is_zero = (x.val() == 0.0);
2108  for (int i=0; i<x.size(); i++)
2109  is_zero = is_zero && (x.dx(i) == 0.0);
2110  return !is_zero;
2111  }
2112 
2113  } // namespace Exp
2114  } // namespace Fad
2115 
2116 } // namespace Sacado
2117 
2118 #define FAD_BOOL_MACRO(OP) \
2119 namespace Sacado { \
2120  namespace Fad { \
2121  namespace Exp { \
2122  template <typename T1, typename T2> \
2123  SACADO_INLINE_FUNCTION \
2124  bool \
2125  operator OP (const Expr<T1>& expr1, \
2126  const Expr<T2>& expr2) \
2127  { \
2128  return toBool(expr1) OP toBool(expr2); \
2129  } \
2130  \
2131  template <typename T2> \
2132  SACADO_INLINE_FUNCTION \
2133  bool \
2134  operator OP (const typename Expr<T2>::value_type& a, \
2135  const Expr<T2>& expr2) \
2136  { \
2137  return a OP toBool(expr2); \
2138  } \
2139  \
2140  template <typename T1> \
2141  SACADO_INLINE_FUNCTION \
2142  bool \
2143  operator OP (const Expr<T1>& expr1, \
2144  const typename Expr<T1>::value_type& b) \
2145  { \
2146  return toBool(expr1) OP b; \
2147  } \
2148  } \
2149  } \
2150 }
2151 
2152 FAD_BOOL_MACRO(&&)
2153 FAD_BOOL_MACRO(||)
2154 
2155 #undef FAD_BOOL_MACRO
2156 
2157 //-------------------------- I/O Operators -----------------------
2158 
2159 namespace Sacado {
2160 
2161  namespace Fad {
2162  namespace Exp {
2163 
2164  template <typename T>
2165  std::ostream& operator << (std::ostream& os, const Expr<T>& xx) {
2166  const typename Expr<T>::derived_type& x = xx.derived();
2167  os << x.val() << " [";
2168 
2169  for (int i=0; i< x.size(); i++) {
2170  os << " " << x.dx(i);
2171  }
2172 
2173  os << " ]";
2174  return os;
2175  }
2176 
2177  } // namespace Exp
2178  } // namespace Fad
2179 
2180 } // namespace Sacado
2181 
2182 #if defined(HAVE_SACADO_KOKKOS)
2183 
2184 //-------------------------- Atomic Operators -----------------------
2185 
2186 namespace Sacado {
2187 
2188  namespace Fad {
2189  namespace Exp {
2190 
2191  // Overload of Kokkos::atomic_add for Fad types.
2192  template <typename S>
2194  void atomic_add(GeneralFad<S>* dst, const GeneralFad<S>& x) {
2195  using Kokkos::atomic_add;
2196 
2197  const int xsz = x.size();
2198  const int sz = dst->size();
2199 
2200  // We currently cannot handle resizing since that would need to be
2201  // done atomically.
2202  if (xsz > sz)
2203  Kokkos::abort(
2204  "Sacado error: Fad resize within atomic_add() not supported!");
2205 
2206  if (xsz != sz && sz > 0 && xsz > 0)
2207  Kokkos::abort(
2208  "Sacado error: Fad assignment of incompatiable sizes!");
2209 
2210 
2211  if (sz > 0 && xsz > 0) {
2213  atomic_add(&(dst->fastAccessDx(i)), x.fastAccessDx(i));
2214  }
2216  atomic_add(&(dst->val()), x.val());
2217  }
2218 
2219  } // namespace Exp
2220  } // namespace Fad
2221 
2222 } // namespace Sacado
2223 
2224 #endif // HAVE_SACADO_KOKKOS
2225 
2226 #endif // SACADO_FAD_OPS_HPP
Wrapper for a generic expression template.
cbrt(expr.val())
#define FAD_BOOL_MACRO(OP)
expr expr SinOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
asinh(expr.val())
char c_
SACADO_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E >::value_type type
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
asin(expr.val())
cosh(expr.val())
#define SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP)
expr expr dx(i)
abs(expr.val())
SACADO_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const ConstT &c_, const T2 &expr2_)
Base template specification for ScalarType.
#define SACADO_FAD_THREAD_SINGLE
atanh(expr.val())
expr expr CoshOp
expr expr ATanhOp
expr expr TanhOp
expr expr SqrtOp
expr expr ASinhOp
Determine whether a given type is an expression.
expr true
static const bool value
atan(expr.val())
Wrapper for a generic expression template.
Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E >::scalar_type type
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
Is a type an expression.
Base template specification for IsSimdType.
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 val()
tanh(expr.val())
expr expr CosOp
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
#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 SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP)
SACADO_INLINE_FUNCTION const derived_type & derived() const
Return derived object.
expr expr ATanOp
#define T2(r, f)
Definition: Sacado_rad.hpp:558
SACADO_INLINE_FUNCTION T if_then_else(const Cond cond, const T &a, const T &b)
SimpleFad< ValueT > min(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ACosOp
#define SACADO_FAD_DERIV_LOOP(I, SZ)
Get the base Fad type from a view/expression.
T derived_type
Typename of derived object, returned by derived()
sqrt(expr.val())
sinh(expr.val())
tan(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
SACADO_INLINE_FUNCTION mpl::enable_if_c< IsFadExpr< T1 >::value &&IsFadExpr< T2 >::value &&ExprLevel< T1 >::value==ExprLevel< T2 >::value, IfThenElseOp< CondT, typename Expr< T1 >::derived_type, typename Expr< T2 >::derived_type, false, false, typename T1::expr_spec_type > >::type if_then_else(const CondT &cond, const T1 &expr1, const T2 &expr2)
SACADO_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const T2 &expr2_)
#define T1(r, f)
Definition: Sacado_rad.hpp:583
static const bool value
Meta-function for determining nesting with an expression.
#define FAD_RELOP_MACRO(OP)
atan2(expr1.val(), expr2.val())
#define SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP)
sin(expr.val())
int value
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
SACADO_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const ConstT &c_)
SACADO_INLINE_FUNCTION T safe_sqrt(const T &x)
log(expr.val())
expr expr ACoshOp
expr expr Log10Op
expr expr SinhOp
SACADO_INLINE_FUNCTION bool toBool(const Expr< T > &xx)
acosh(expr.val())
if_then_else(expr.val() >=0, expr.dx(i), value_type(-expr.dx(i)))
acos(expr.val())
SimpleFad< ValueT > max(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ASinOp
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)
expr expr expr bar false
#define SACADO_INLINE_FUNCTION
exp(expr.val())
expr expr expr ExpOp
fabs(expr.val())
expr expr AbsOp
expr expr TanOp
log10(expr.val())
Base template specification for ValueType.
Base template specification for Promote.
cos(expr.val())
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