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 //
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_KOKKOS)
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) \
50 namespace 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 
155 FAD_UNARYOP_MACRO(operator+,
156  UnaryPlusOp,
157  ;,
158  expr.val(),
159  expr.dx(i),
160  expr.fastAccessDx(i))
161 FAD_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()))
237 FAD_UNARYOP_MACRO(sinh,
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.
294 namespace 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) \
459 namespace 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 
747 FAD_BINARYOP_MACRO(operator+,
748  AdditionOp,
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))
761 FAD_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))
775 FAD_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)
790 FAD_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.
873 namespace 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 
1712 namespace 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>
1719  class IfThenElseOp {};
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 {
1758  using Sacado::if_then_else;
1759  return if_then_else( cond, expr1.val(), expr2.val() );
1760  }
1761 
1763  value_type dx(int i) const {
1764  using Sacado::if_then_else;
1765  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
1766  }
1767 
1770  using Sacado::if_then_else;
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>
1783  class IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > :
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 {
1811  using Sacado::if_then_else;
1812  return if_then_else( cond, expr1.val(), c );
1813  }
1814 
1816  value_type dx(int i) const {
1817  using Sacado::if_then_else;
1818  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
1819  }
1820 
1823  using Sacado::if_then_else;
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>
1835  class IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > :
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 {
1863  using Sacado::if_then_else;
1864  return if_then_else( cond, c, expr2.val() );
1865  }
1866 
1868  value_type dx(int i) const {
1869  using Sacado::if_then_else;
1870  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
1871  }
1872 
1875  using Sacado::if_then_else;
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>
1890  IfThenElseOp< CondT,
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  std::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  std::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 
2028 namespace 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) \
2051 namespace 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 
2088 FAD_RELOP_MACRO(==)
2089 FAD_RELOP_MACRO(!=)
2090 FAD_RELOP_MACRO(<)
2091 FAD_RELOP_MACRO(>)
2092 FAD_RELOP_MACRO(<=)
2093 FAD_RELOP_MACRO(>=)
2094 FAD_RELOP_MACRO(<<=)
2095 FAD_RELOP_MACRO(>>=)
2096 FAD_RELOP_MACRO(&)
2097 FAD_RELOP_MACRO(|)
2098 
2099 #undef FAD_RELOP_MACRO
2100 
2101 namespace 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 -----------------------
2119 namespace 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) \
2140 namespace 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 
2173 FAD_BOOL_MACRO(&&)
2174 FAD_BOOL_MACRO(||)
2175 
2176 #undef FAD_BOOL_MACRO
2177 
2178 //-------------------------- I/O Operators -----------------------
2179 
2180 namespace 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_KOKKOS)
2204 
2205 //-------------------------- Atomic Operators -----------------------
2206 
2207 namespace 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_KOKKOS
2246 
2247 #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:578
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:603
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