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 
42 #if defined(HAVE_SACADO_KOKKOSCORE)
43 #include "Kokkos_Atomic.hpp"
44 #include "impl/Kokkos_Error.hpp"
45 #endif
46 
47 #define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
48 namespace Sacado { \
49  namespace Fad { \
50  namespace Exp { \
51  \
52  template <typename T, typename ExprSpec> \
53  class OP {}; \
54  \
55  template <typename T> \
56  class OP< T,ExprSpecDefault > : \
57  public Expr< OP<T,ExprSpecDefault> > { \
58  public: \
59  \
60  typedef typename std::remove_cv<T>::type ExprT; \
61  typedef typename ExprT::value_type value_type; \
62  typedef typename ExprT::scalar_type scalar_type; \
63  \
64  typedef ExprSpecDefault expr_spec_type; \
65  \
66  KOKKOS_INLINE_FUNCTION \
67  OP(const T& expr_) : expr(expr_) {} \
68  \
69  KOKKOS_INLINE_FUNCTION \
70  int size() const { return expr.size(); } \
71  \
72  KOKKOS_INLINE_FUNCTION \
73  bool hasFastAccess() const { \
74  return expr.hasFastAccess(); \
75  } \
76  \
77  KOKKOS_INLINE_FUNCTION \
78  value_type val() const { \
79  USING \
80  return VALUE; \
81  } \
82  \
83  KOKKOS_INLINE_FUNCTION \
84  value_type dx(int i) const { \
85  USING \
86  return DX; \
87  } \
88  \
89  KOKKOS_INLINE_FUNCTION \
90  value_type fastAccessDx(int i) const { \
91  USING \
92  return FASTACCESSDX; \
93  } \
94  \
95  protected: \
96  \
97  const T& expr; \
98  }; \
99  \
100  template <typename T> \
101  KOKKOS_INLINE_FUNCTION \
102  OP< typename Expr<T>::derived_type, \
103  typename T::expr_spec_type > \
104  OPNAME (const Expr<T>& expr) \
105  { \
106  typedef OP< typename Expr<T>::derived_type, \
107  typename T::expr_spec_type > expr_t; \
108  \
109  return expr_t(expr.derived()); \
110  } \
111  \
112  template <typename T, typename E> \
113  struct ExprLevel< OP< T,E > > { \
114  static const unsigned value = ExprLevel<T>::value; \
115  }; \
116  \
117  template <typename T, typename E> \
118  struct IsFadExpr< OP< T,E > > { \
119  static const unsigned value = true; \
120  }; \
121  \
122  } \
123  } \
124  \
125  template <typename T, typename E> \
126  struct IsExpr< Fad::Exp::OP< T,E > > { \
127  static const bool value = true; \
128  }; \
129  \
130  template <typename T, typename E> \
131  struct BaseExprType< Fad::Exp::OP< T,E > > { \
132  typedef typename BaseExprType<T>::type type; \
133  }; \
134  \
135  template <typename T, typename E> \
136  struct IsSimdType< Fad::Exp::OP< T,E > > { \
137  static const bool value = \
138  IsSimdType< typename Fad::Exp::OP< T,E >::scalar_type >::value; \
139  }; \
140  \
141 }
142 
143 FAD_UNARYOP_MACRO(operator+,
144  UnaryPlusOp,
145  ;,
146  expr.val(),
147  expr.dx(i),
148  expr.fastAccessDx(i))
149 FAD_UNARYOP_MACRO(operator-,
151  ;,
152  -expr.val(),
153  -expr.dx(i),
154  -expr.fastAccessDx(i))
157  using std::exp;,
158  exp(expr.val()),
159  exp(expr.val())*expr.dx(i),
160  exp(expr.val())*expr.fastAccessDx(i))
163  using std::log;,
164  log(expr.val()),
165  expr.dx(i)/expr.val(),
166  expr.fastAccessDx(i)/expr.val())
169  using std::log10; using std::log;,
170  log10(expr.val()),
171  expr.dx(i)/( log(value_type(10))*expr.val()),
172  expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
175  using std::sqrt;,
176  sqrt(expr.val()),
177  expr.dx(i)/(value_type(2)* sqrt(expr.val())),
178  expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
181  using std::cos; using std::sin;,
182  cos(expr.val()),
183  -expr.dx(i)* sin(expr.val()),
184  -expr.fastAccessDx(i)* sin(expr.val()))
187  using std::cos; using std::sin;,
188  sin(expr.val()),
189  expr.dx(i)* cos(expr.val()),
190  expr.fastAccessDx(i)* cos(expr.val()))
193  using std::tan;,
194  tan(expr.val()),
195  expr.dx(i)*
196  (value_type(1)+ tan(expr.val())* tan(expr.val())),
197  expr.fastAccessDx(i)*
198  (value_type(1)+ tan(expr.val())* tan(expr.val())))
201  using std::acos; using std::sqrt;,
202  acos(expr.val()),
203  -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
204  -expr.fastAccessDx(i) /
205  sqrt(value_type(1)-expr.val()*expr.val()))
208  using std::asin; using std::sqrt;,
209  asin(expr.val()),
210  expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
211  expr.fastAccessDx(i) /
212  sqrt(value_type(1)-expr.val()*expr.val()))
215  using std::atan;,
216  atan(expr.val()),
217  expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
218  expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
221  using std::cosh; using std::sinh;,
222  cosh(expr.val()),
223  expr.dx(i)* sinh(expr.val()),
224  expr.fastAccessDx(i)* sinh(expr.val()))
225 FAD_UNARYOP_MACRO(sinh,
227  using std::cosh; using std::sinh;,
228  sinh(expr.val()),
229  expr.dx(i)* cosh(expr.val()),
230  expr.fastAccessDx(i)* cosh(expr.val()))
233  using std::tanh; using std::cosh;,
234  tanh(expr.val()),
235  expr.dx(i)/( cosh(expr.val())* cosh(expr.val())),
236  expr.fastAccessDx(i) /
237  ( cosh(expr.val())* cosh(expr.val())))
240  using std::acosh; using std::sqrt;,
241  acosh(expr.val()),
242  expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
243  (expr.val()+value_type(1))),
244  expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
245  (expr.val()+value_type(1))))
248  using std::asinh; using std::sqrt;,
249  asinh(expr.val()),
250  expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
251  expr.fastAccessDx(i)/ sqrt(value_type(1)+
252  expr.val()*expr.val()))
255  using std::atanh;,
256  atanh(expr.val()),
257  expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
258  expr.fastAccessDx(i)/(value_type(1)-
259  expr.val()*expr.val()))
262  using std::abs;,
263  abs(expr.val()),
264  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
265  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
268  using std::fabs;,
269  fabs(expr.val()),
270  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
271  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
274  using std::cbrt;,
275  cbrt(expr.val()),
276  expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
277  expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
278 
279 #undef FAD_UNARYOP_MACRO
280 
281 #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) \
282 namespace Sacado { \
283  namespace Fad { \
284  namespace Exp { \
285  \
286  template <typename T1, typename T2, \
287  bool is_const_T1, bool is_const_T2, \
288  typename ExprSpec > \
289  class OP {}; \
290  \
291  template <typename T1, typename T2> \
292  class OP< T1, T2, false, false, ExprSpecDefault > : \
293  public Expr< OP< T1, T2, false, false, ExprSpecDefault > > { \
294  public: \
295  \
296  typedef typename std::remove_cv<T1>::type ExprT1; \
297  typedef typename std::remove_cv<T2>::type ExprT2; \
298  typedef typename ExprT1::value_type value_type_1; \
299  typedef typename ExprT2::value_type value_type_2; \
300  typedef typename Sacado::Promote<value_type_1, \
301  value_type_2>::type value_type; \
302  \
303  typedef typename ExprT1::scalar_type scalar_type_1; \
304  typedef typename ExprT2::scalar_type scalar_type_2; \
305  typedef typename Sacado::Promote<scalar_type_1, \
306  scalar_type_2>::type scalar_type; \
307  \
308  typedef ExprSpecDefault expr_spec_type; \
309  \
310  KOKKOS_INLINE_FUNCTION \
311  OP(const T1& expr1_, const T2& expr2_) : \
312  expr1(expr1_), expr2(expr2_) {} \
313  \
314  KOKKOS_INLINE_FUNCTION \
315  int size() const { \
316  const int sz1 = expr1.size(), sz2 = expr2.size(); \
317  return sz1 > sz2 ? sz1 : sz2; \
318  } \
319  \
320  KOKKOS_INLINE_FUNCTION \
321  bool hasFastAccess() const { \
322  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
323  } \
324  \
325  KOKKOS_INLINE_FUNCTION \
326  value_type val() const { \
327  USING \
328  return VALUE; \
329  } \
330  \
331  KOKKOS_INLINE_FUNCTION \
332  value_type dx(int i) const { \
333  USING \
334  const int sz1 = expr1.size(), sz2 = expr2.size(); \
335  if (sz1 > 0 && sz2 > 0) \
336  return DX; \
337  else if (sz1 > 0) \
338  return CDX2; \
339  else \
340  return CDX1; \
341  } \
342  \
343  KOKKOS_INLINE_FUNCTION \
344  value_type fastAccessDx(int i) const { \
345  USING \
346  return FASTACCESSDX; \
347  } \
348  \
349  protected: \
350  \
351  const T1& expr1; \
352  const T2& expr2; \
353  \
354  }; \
355  \
356  template <typename T1, typename T2> \
357  class OP< T1, T2, false, true, ExprSpecDefault > \
358  : public Expr< OP< T1, T2, false, true, ExprSpecDefault > > { \
359  public: \
360  \
361  typedef typename std::remove_cv<T1>::type ExprT1; \
362  typedef T2 ConstT; \
363  typedef typename ExprT1::value_type value_type; \
364  typedef typename ExprT1::scalar_type scalar_type; \
365  \
366  typedef ExprSpecDefault expr_spec_type; \
367  \
368  KOKKOS_INLINE_FUNCTION \
369  OP(const T1& expr1_, const ConstT& c_) : \
370  expr1(expr1_), c(c_) {} \
371  \
372  KOKKOS_INLINE_FUNCTION \
373  int size() const { \
374  return expr1.size(); \
375  } \
376  \
377  KOKKOS_INLINE_FUNCTION \
378  bool hasFastAccess() const { \
379  return expr1.hasFastAccess(); \
380  } \
381  \
382  KOKKOS_INLINE_FUNCTION \
383  value_type val() const { \
384  USING \
385  return VAL_CONST_DX_2; \
386  } \
387  \
388  KOKKOS_INLINE_FUNCTION \
389  value_type dx(int i) const { \
390  USING \
391  return CONST_DX_2; \
392  } \
393  \
394  KOKKOS_INLINE_FUNCTION \
395  value_type fastAccessDx(int i) const { \
396  USING \
397  return CONST_FASTACCESSDX_2; \
398  } \
399  \
400  protected: \
401  \
402  const T1& expr1; \
403  const ConstT& c; \
404  }; \
405  \
406  template <typename T1, typename T2> \
407  class OP< T1, T2, true, false, ExprSpecDefault > \
408  : public Expr< OP< T1, T2, true, false, ExprSpecDefault > > { \
409  public: \
410  \
411  typedef typename std::remove_cv<T2>::type ExprT2; \
412  typedef T1 ConstT; \
413  typedef typename ExprT2::value_type value_type; \
414  typedef typename ExprT2::scalar_type scalar_type; \
415  \
416  typedef ExprSpecDefault expr_spec_type; \
417  \
418  KOKKOS_INLINE_FUNCTION \
419  OP(const ConstT& c_, const T2& expr2_) : \
420  c(c_), expr2(expr2_) {} \
421  \
422  KOKKOS_INLINE_FUNCTION \
423  int size() const { \
424  return expr2.size(); \
425  } \
426  \
427  KOKKOS_INLINE_FUNCTION \
428  bool hasFastAccess() const { \
429  return expr2.hasFastAccess(); \
430  } \
431  \
432  KOKKOS_INLINE_FUNCTION \
433  value_type val() const { \
434  USING \
435  return VAL_CONST_DX_1; \
436  } \
437  \
438  KOKKOS_INLINE_FUNCTION \
439  value_type dx(int i) const { \
440  USING \
441  return CONST_DX_1; \
442  } \
443  \
444  KOKKOS_INLINE_FUNCTION \
445  value_type fastAccessDx(int i) const { \
446  USING \
447  return CONST_FASTACCESSDX_1; \
448  } \
449  \
450  protected: \
451  \
452  const ConstT& c; \
453  const T2& expr2; \
454  }; \
455  \
456  template <typename T1, typename T2> \
457  KOKKOS_INLINE_FUNCTION \
458  SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP) \
459  OPNAME (const T1& expr1, const T2& expr2) \
460  { \
461  typedef OP< typename Expr<T1>::derived_type, \
462  typename Expr<T2>::derived_type, \
463  false, false, typename T1::expr_spec_type > expr_t; \
464  \
465  return expr_t(expr1.derived(), expr2.derived()); \
466  } \
467  \
468  template <typename T> \
469  KOKKOS_INLINE_FUNCTION \
470  OP< typename T::value_type, typename Expr<T>::derived_type, \
471  true, false, typename T::expr_spec_type > \
472  OPNAME (const typename T::value_type& c, \
473  const Expr<T>& expr) \
474  { \
475  typedef typename T::value_type ConstT; \
476  typedef OP< ConstT, typename Expr<T>::derived_type, \
477  true, false, typename T::expr_spec_type > expr_t; \
478  \
479  return expr_t(c, expr.derived()); \
480  } \
481  \
482  template <typename T> \
483  KOKKOS_INLINE_FUNCTION \
484  OP< typename Expr<T>::derived_type, typename T::value_type, \
485  false, true, typename T::expr_spec_type > \
486  OPNAME (const Expr<T>& expr, \
487  const typename T::value_type& c) \
488  { \
489  typedef typename T::value_type ConstT; \
490  typedef OP< typename Expr<T>::derived_type, ConstT, \
491  false, true, typename T::expr_spec_type > expr_t; \
492  \
493  return expr_t(expr.derived(), c); \
494  } \
495  \
496  template <typename T> \
497  KOKKOS_INLINE_FUNCTION \
498  SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP) \
499  OPNAME (const typename T::scalar_type& c, \
500  const Expr<T>& expr) \
501  { \
502  typedef typename T::scalar_type ConstT; \
503  typedef OP< ConstT, typename Expr<T>::derived_type, \
504  true, false, typename T::expr_spec_type > expr_t; \
505  \
506  return expr_t(c, expr.derived()); \
507  } \
508  \
509  template <typename T> \
510  KOKKOS_INLINE_FUNCTION \
511  SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP) \
512  OPNAME (const Expr<T>& expr, \
513  const typename T::scalar_type& c) \
514  { \
515  typedef typename T::scalar_type ConstT; \
516  typedef OP< typename Expr<T>::derived_type, ConstT, \
517  false, true, typename T::expr_spec_type > expr_t; \
518  \
519  return expr_t(expr.derived(), c); \
520  } \
521  \
522  template <typename T1, typename T2, bool c1, bool c2, typename E> \
523  struct ExprLevel< OP< T1, T2, c1, c2, E > > { \
524  static constexpr unsigned value_1 = ExprLevel<T1>::value; \
525  static constexpr unsigned value_2 = ExprLevel<T2>::value; \
526  static constexpr unsigned value = \
527  value_1 >= value_2 ? value_1 : value_2; \
528  }; \
529  \
530  template <typename T1, typename T2, bool c1, bool c2, typename E> \
531  struct IsFadExpr< OP< T1, T2, c1, c2, E > > { \
532  static constexpr unsigned value = true; \
533  }; \
534  \
535  } \
536  } \
537  \
538  template <typename T1, typename T2, bool c1, bool c2, typename E> \
539  struct IsExpr< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
540  static constexpr bool value = true; \
541  }; \
542  \
543  template <typename T1, typename T2, bool c1, bool c2, typename E> \
544  struct BaseExprType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
545  typedef typename BaseExprType<T1>::type base_expr_1; \
546  typedef typename BaseExprType<T2>::type base_expr_2; \
547  typedef typename Sacado::Promote<base_expr_1, \
548  base_expr_2>::type type; \
549  }; \
550  \
551  template <typename T1, typename T2, bool c1, bool c2, typename E> \
552  struct IsSimdType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
553  static const bool value = \
554  IsSimdType< typename Fad::Exp::OP< T1, T2, c1, c2, E >::value_type >::value; \
555  }; \
556  \
557 }
558 
559 
560 FAD_BINARYOP_MACRO(operator+,
562  ;,
563  expr1.val() + expr2.val(),
564  expr1.dx(i) + expr2.dx(i),
565  expr2.dx(i),
566  expr1.dx(i),
567  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
568  c + expr2.val(),
569  expr1.val() + c,
570  expr2.dx(i),
571  expr1.dx(i),
572  expr2.fastAccessDx(i),
573  expr1.fastAccessDx(i))
574 FAD_BINARYOP_MACRO(operator-,
576  ;,
577  expr1.val() - expr2.val(),
578  expr1.dx(i) - expr2.dx(i),
579  -expr2.dx(i),
580  expr1.dx(i),
581  expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
582  c - expr2.val(),
583  expr1.val() - c,
584  -expr2.dx(i),
585  expr1.dx(i),
586  -expr2.fastAccessDx(i),
587  expr1.fastAccessDx(i))
588 FAD_BINARYOP_MACRO(operator*,
590  ;,
591  expr1.val() * expr2.val(),
592  expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
593  expr1.val()*expr2.dx(i),
594  expr1.dx(i)*expr2.val(),
595  expr1.val()*expr2.fastAccessDx(i) +
596  expr1.fastAccessDx(i)*expr2.val(),
597  c * expr2.val(),
598  expr1.val() * c,
599  c*expr2.dx(i),
600  expr1.dx(i)*c,
601  c*expr2.fastAccessDx(i),
602  expr1.fastAccessDx(i)*c)
603 FAD_BINARYOP_MACRO(operator/,
605  ;,
606  expr1.val() / expr2.val(),
607  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
608  (expr2.val()*expr2.val()),
609  -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
610  expr1.dx(i)/expr2.val(),
611  (expr1.fastAccessDx(i)*expr2.val() -
612  expr2.fastAccessDx(i)*expr1.val()) /
613  (expr2.val()*expr2.val()),
614  c / expr2.val(),
615  expr1.val() / c,
616  -expr2.dx(i)*c / (expr2.val()*expr2.val()),
617  expr1.dx(i)/c,
618  -expr2.fastAccessDx(i)*c / (expr2.val()*expr2.val()),
619  expr1.fastAccessDx(i)/c)
622  using std::atan2;,
623  atan2(expr1.val(), expr2.val()),
624  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
625  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
626  -expr1.val()*expr2.dx(i)/
627  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
628  expr2.val()*expr1.dx(i)/
629  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
630  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
631  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
632  atan2(c, expr2.val()),
633  atan2(expr1.val(), c),
634  (-c*expr2.dx(i)) / (c*c + expr2.val()*expr2.val()),
635  (c*expr1.dx(i))/ (expr1.val()*expr1.val() + c*c),
636  (-c*expr2.fastAccessDx(i))/ (c*c + expr2.val()*expr2.val()),
637  (c*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c*c))
638 // FAD_BINARYOP_MACRO(pow,
639 // PowerOp,
640 // using std::pow; using std::log;,
641 // pow(expr1.val(), expr2.val()),
642 // 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())) ),
643 // 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())) ),
644 // 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())) ),
645 // 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())) ),
646 // pow(c, expr2.val()),
647 // pow(expr1.val(), c),
648 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) ),
649 // 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)) ),
650 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) ),
651 // 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))) )
654  ;,
655  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
656  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
657  if_then_else( expr1.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
658  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), value_type(0.0) ),
659  if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
660  if_then_else( c >= expr2.val(), value_type(c), expr2.val() ),
661  if_then_else( expr1.val() >= c, expr1.val(), value_type(c) ),
662  if_then_else( c >= expr2.val(), value_type(0.0), expr2.dx(i) ),
663  if_then_else( expr1.val() >= c, expr1.dx(i), value_type(0.0) ),
664  if_then_else( c >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
665  if_then_else( expr1.val() >= c, expr1.fastAccessDx(i), value_type(0.0) ) )
668  ;,
669  if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
670  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
671  if_then_else( expr1.val() <= expr2.val(), value_type(0.0), expr2.dx(i) ),
672  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), value_type(0.0) ),
673  if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
674  if_then_else( c <= expr2.val(), value_type(c), expr2.val() ),
675  if_then_else( expr1.val() <= c, expr1.val(), value_type(c) ),
676  if_then_else( c <= expr2.val(), value_type(0), expr2.dx(i) ),
677  if_then_else( expr1.val() <= c, expr1.dx(i), value_type(0) ),
678  if_then_else( c <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
679  if_then_else( expr1.val() <= c, expr1.fastAccessDx(i), value_type(0) ) )
680 
681 
682 // Special handling for std::pow() to provide specializations of PowerOp for
683 // "simd" value types that use if_then_else(). The only reason for not using
684 // if_then_else() always is to avoid evaluating the derivative if the value is
685 // zero to avoid throwing FPEs.
686 namespace Sacado {
687  namespace Fad {
688  namespace Exp {
689 
690  template <typename T1, typename T2,
691  bool is_const_T1, bool is_const_T2,
692  typename ExprSpec, bool is_simd >
693  class PowerOp {};
694 
695  //
696  // Implementation for simd type using if_then_else()
697  //
698  template <typename T1, typename T2>
699  class PowerOp< T1, T2, false, false, ExprSpecDefault, true > :
700  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault, true > > {
701  public:
702 
703  typedef typename std::remove_cv<T1>::type ExprT1;
704  typedef typename std::remove_cv<T2>::type ExprT2;
705  typedef typename ExprT1::value_type value_type_1;
706  typedef typename ExprT2::value_type value_type_2;
707  typedef typename Sacado::Promote<value_type_1,
708  value_type_2>::type value_type;
709 
710  typedef typename ExprT1::scalar_type scalar_type_1;
711  typedef typename ExprT2::scalar_type scalar_type_2;
712  typedef typename Sacado::Promote<scalar_type_1,
713  scalar_type_2>::type scalar_type;
714 
715  typedef ExprSpecDefault expr_spec_type;
716 
718  PowerOp(const T1& expr1_, const T2& expr2_) :
719  expr1(expr1_), expr2(expr2_) {}
720 
722  int size() const {
723  const int sz1 = expr1.size(), sz2 = expr2.size();
724  return sz1 > sz2 ? sz1 : sz2;
725  }
726 
728  bool hasFastAccess() const {
729  return expr1.hasFastAccess() && expr2.hasFastAccess();
730  }
731 
733  value_type val() const {
734  using std::pow;
735  return pow(expr1.val(), expr2.val());
736  }
737 
739  value_type dx(int i) const {
740  using std::pow; using std::log;
741  const int sz1 = expr1.size(), sz2 = expr2.size();
742  if (sz1 > 0 && sz2 > 0)
743  return 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())) );
744  else if (sz1 > 0)
745  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
746  // It seems less accurate and caused convergence problems in some codes
747  return 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())) );
748  else
749  return 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())) );
750  }
751 
753  value_type fastAccessDx(int i) const {
754  using std::pow; using std::log;
755  return 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())) );
756  }
757 
758  protected:
759 
760  const T1& expr1;
761  const T2& expr2;
762 
763  };
764 
765  template <typename T1, typename T2>
766  class PowerOp< T1, T2, false, true, ExprSpecDefault, true >
767  : public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault, true > > {
768  public:
769 
770  typedef typename std::remove_cv<T1>::type ExprT1;
771  typedef T2 ConstT;
772  typedef typename ExprT1::value_type value_type;
773  typedef typename ExprT1::scalar_type scalar_type;
774 
775  typedef ExprSpecDefault expr_spec_type;
776 
778  PowerOp(const T1& expr1_, const ConstT& c_) :
779  expr1(expr1_), c(c_) {}
780 
782  int size() const {
783  return expr1.size();
784  }
785 
787  bool hasFastAccess() const {
788  return expr1.hasFastAccess();
789  }
790 
792  value_type val() const {
793  using std::pow;
794  return pow(expr1.val(), c);
795  }
796 
798  value_type dx(int i) const {
799  using std::pow;
800  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
801  // It seems less accurate and caused convergence problems in some codes
802  return 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)) );
803  }
804 
806  value_type fastAccessDx(int i) const {
807  using std::pow;
808  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
809  // It seems less accurate and caused convergence problems in some codes
810  return 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)) );
811  }
812 
813  protected:
814 
815  const T1& expr1;
816  const ConstT& c;
817  };
818 
819  template <typename T1, typename T2>
820  class PowerOp< T1, T2, true, false, ExprSpecDefault, true >
821  : public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault, true > > {
822  public:
823 
824  typedef typename std::remove_cv<T2>::type ExprT2;
825  typedef T1 ConstT;
826  typedef typename ExprT2::value_type value_type;
827  typedef typename ExprT2::scalar_type scalar_type;
828 
829  typedef ExprSpecDefault expr_spec_type;
830 
832  PowerOp(const ConstT& c_, const T2& expr2_) :
833  c(c_), expr2(expr2_) {}
834 
836  int size() const {
837  return expr2.size();
838  }
839 
841  bool hasFastAccess() const {
842  return expr2.hasFastAccess();
843  }
844 
846  value_type val() const {
847  using std::pow;
848  return pow(c, expr2.val());
849  }
850 
852  value_type dx(int i) const {
853  using std::pow; using std::log;
854  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) );
855  }
856 
858  value_type fastAccessDx(int i) const {
859  using std::pow; using std::log;
860  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) );
861  }
862 
863  protected:
864 
865  const ConstT& c;
866  const T2& expr2;
867  };
868 
869  //
870  // Specialization for scalar types using ternary operator
871  //
872 
873  template <typename T1, typename T2>
874  class PowerOp< T1, T2, false, false, ExprSpecDefault, false > :
875  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault, false > > {
876  public:
877 
878  typedef typename std::remove_cv<T1>::type ExprT1;
879  typedef typename std::remove_cv<T2>::type ExprT2;
880  typedef typename ExprT1::value_type value_type_1;
881  typedef typename ExprT2::value_type value_type_2;
882  typedef typename Sacado::Promote<value_type_1,
883  value_type_2>::type value_type;
884 
885  typedef typename ExprT1::scalar_type scalar_type_1;
886  typedef typename ExprT2::scalar_type scalar_type_2;
887  typedef typename Sacado::Promote<scalar_type_1,
888  scalar_type_2>::type scalar_type;
889 
890  typedef ExprSpecDefault expr_spec_type;
891 
893  PowerOp(const T1& expr1_, const T2& expr2_) :
894  expr1(expr1_), expr2(expr2_) {}
895 
897  int size() const {
898  const int sz1 = expr1.size(), sz2 = expr2.size();
899  return sz1 > sz2 ? sz1 : sz2;
900  }
901 
903  bool hasFastAccess() const {
904  return expr1.hasFastAccess() && expr2.hasFastAccess();
905  }
906 
908  value_type val() const {
909  using std::pow;
910  return pow(expr1.val(), expr2.val());
911  }
912 
914  value_type dx(int i) const {
915  using std::pow; using std::log;
916  const int sz1 = expr1.size(), sz2 = expr2.size();
917  if (sz1 > 0 && sz2 > 0)
918  return 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()));
919  else if (sz1 > 0)
920  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
921  // It seems less accurate and caused convergence problems in some codes
922  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val()));
923  else
924  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
925  }
926 
928  value_type fastAccessDx(int i) const {
929  using std::pow; using std::log;
930  return 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()));
931  }
932 
933  protected:
934 
935  const T1& expr1;
936  const T2& expr2;
937 
938  };
939 
940  template <typename T1, typename T2>
941  class PowerOp< T1, T2, false, true, ExprSpecDefault, false >
942  : public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault, false > > {
943  public:
944 
945  typedef typename std::remove_cv<T1>::type ExprT1;
946  typedef T2 ConstT;
947  typedef typename ExprT1::value_type value_type;
948  typedef typename ExprT1::scalar_type scalar_type;
949 
950  typedef ExprSpecDefault expr_spec_type;
951 
953  PowerOp(const T1& expr1_, const ConstT& c_) :
954  expr1(expr1_), c(c_) {}
955 
957  int size() const {
958  return expr1.size();
959  }
960 
962  bool hasFastAccess() const {
963  return expr1.hasFastAccess();
964  }
965 
967  value_type val() const {
968  using std::pow;
969  return pow(expr1.val(), c);
970  }
971 
973  value_type dx(int i) const {
974  using std::pow;
975  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
976  // It seems less accurate and caused convergence problems in some codes
977  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)/expr1.val()*pow(expr1.val(),c));
978  }
979 
981  value_type fastAccessDx(int i) const {
982  using std::pow;
983  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
984  // It seems less accurate and caused convergence problems in some codes
985  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c));
986  }
987 
988  protected:
989 
990  const T1& expr1;
991  const ConstT& c;
992  };
993 
994  template <typename T1, typename T2>
995  class PowerOp< T1, T2, true, false, ExprSpecDefault, false >
996  : public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault, false > > {
997  public:
998 
999  typedef typename std::remove_cv<T2>::type ExprT2;
1000  typedef T1 ConstT;
1001  typedef typename ExprT2::value_type value_type;
1002  typedef typename ExprT2::scalar_type scalar_type;
1003 
1004  typedef ExprSpecDefault expr_spec_type;
1005 
1007  PowerOp(const ConstT& c_, const T2& expr2_) :
1008  c(c_), expr2(expr2_) {}
1009 
1011  int size() const {
1012  return expr2.size();
1013  }
1014 
1016  bool hasFastAccess() const {
1017  return expr2.hasFastAccess();
1018  }
1019 
1021  value_type val() const {
1022  using std::pow;
1023  return pow(c, expr2.val());
1024  }
1025 
1027  value_type dx(int i) const {
1028  using std::pow; using std::log;
1029  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c)*pow(c,expr2.val()));
1030  }
1031 
1033  value_type fastAccessDx(int i) const {
1034  using std::pow; using std::log;
1035  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val()));
1036  }
1037 
1038  protected:
1039 
1040  const ConstT& c;
1041  const T2& expr2;
1042  };
1043 
1044  template <typename T1, typename T2>
1047  pow (const T1& expr1, const T2& expr2)
1048  {
1049  typedef PowerOp< typename Expr<T1>::derived_type,
1050  typename Expr<T2>::derived_type,
1051  false, false, typename T1::expr_spec_type > expr_t;
1052 
1053  return expr_t(expr1.derived(), expr2.derived());
1054  }
1055 
1056  template <typename T>
1058  PowerOp< typename T::value_type, typename Expr<T>::derived_type,
1059  true, false, typename T::expr_spec_type >
1060  pow (const typename T::value_type& c,
1061  const Expr<T>& expr)
1062  {
1063  typedef typename T::value_type ConstT;
1064  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1065  true, false, typename T::expr_spec_type > expr_t;
1066 
1067  return expr_t(c, expr.derived());
1068  }
1069 
1070  template <typename T>
1072  PowerOp< typename Expr<T>::derived_type, typename T::value_type,
1073  false, true, typename T::expr_spec_type >
1074  pow (const Expr<T>& expr,
1075  const typename T::value_type& c)
1076  {
1077  typedef typename T::value_type ConstT;
1078  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1079  false, true, typename T::expr_spec_type > expr_t;
1080 
1081  return expr_t(expr.derived(), c);
1082  }
1083 
1084  template <typename T>
1087  pow (const typename T::scalar_type& c,
1088  const Expr<T>& expr)
1089  {
1090  typedef typename T::scalar_type ConstT;
1091  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1092  true, false, typename T::expr_spec_type > expr_t;
1093 
1094  return expr_t(c, expr.derived());
1095  }
1096 
1097  template <typename T>
1100  pow (const Expr<T>& expr,
1101  const typename T::scalar_type& c)
1102  {
1103  typedef typename T::scalar_type ConstT;
1104  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1105  false, true, typename T::expr_spec_type > expr_t;
1106 
1107  return expr_t(expr.derived(), c);
1108  }
1109 
1110  template <typename T1, typename T2, bool c1, bool c2, typename E>
1111  struct ExprLevel< PowerOp< T1, T2, c1, c2, E > > {
1112  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1113  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1114  static constexpr unsigned value =
1115  value_1 >= value_2 ? value_1 : value_2;
1116  };
1117 
1118  template <typename T1, typename T2, bool c1, bool c2, typename E>
1119  struct IsFadExpr< PowerOp< T1, T2, c1, c2, E > > {
1120  static constexpr unsigned value = true;
1121  };
1122 
1123  }
1124  }
1125 
1126  template <typename T1, typename T2, bool c1, bool c2, typename E>
1127  struct IsExpr< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1128  static constexpr bool value = true;
1129  };
1130 
1131  template <typename T1, typename T2, bool c1, bool c2, typename E>
1132  struct BaseExprType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1133  typedef typename BaseExprType<T1>::type base_expr_1;
1134  typedef typename BaseExprType<T2>::type base_expr_2;
1135  typedef typename Sacado::Promote<base_expr_1,
1136  base_expr_2>::type type;
1137  };
1138 
1139  template <typename T1, typename T2, bool c1, bool c2, typename E>
1140  struct IsSimdType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1141  static const bool value =
1142  IsSimdType< typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::value_type >::value;
1143  };
1144 
1145 }
1146 
1147 //--------------------------if_then_else operator -----------------------
1148 // Can't use the above macros because it is a ternary operator (sort of).
1149 // Also, relies on C++11
1150 
1151 namespace Sacado {
1152  namespace Fad {
1153  namespace Exp {
1154 
1155  template <typename CondT, typename T1, typename T2,
1156  bool is_const_T1, bool is_const_T2,
1157  typename ExprSpec>
1158  class IfThenElseOp {};
1159 
1160  template <typename CondT, typename T1, typename T2>
1162  public Expr< IfThenElseOp< CondT, T1, T2, false, false, ExprSpecDefault > > {
1163 
1164  public:
1165 
1166  typedef typename std::remove_cv<T1>::type ExprT1;
1167  typedef typename std::remove_cv<T2>::type ExprT2;
1168  typedef typename ExprT1::value_type value_type_1;
1169  typedef typename ExprT2::value_type value_type_2;
1170  typedef typename Sacado::Promote<value_type_1,
1172 
1173  typedef typename ExprT1::scalar_type scalar_type_1;
1174  typedef typename ExprT2::scalar_type scalar_type_2;
1175  typedef typename Sacado::Promote<scalar_type_1,
1177 
1179 
1181  IfThenElseOp(const CondT& cond_, const T1& expr1_, const T2& expr2_) :
1182  cond(cond_), expr1(expr1_), expr2(expr2_) {}
1183 
1185  int size() const {
1186  int sz1 = expr1.size(), sz2 = expr2.size();
1187  return sz1 > sz2 ? sz1 : sz2;
1188  }
1189 
1191  bool hasFastAccess() const {
1192  return expr1.hasFastAccess() && expr2.hasFastAccess();
1193  }
1194 
1196  value_type val() const {
1197  return if_then_else( cond, expr1.val(), expr2.val() );
1198  }
1199 
1201  value_type dx(int i) const {
1202  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
1203  }
1204 
1206  value_type fastAccessDx(int i) const {
1207  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
1208  }
1209 
1210  protected:
1211 
1212  const CondT& cond;
1213  const T1& expr1;
1214  const T2& expr2;
1215 
1216  };
1217 
1218  template <typename CondT, typename T1, typename T2>
1219  class IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > :
1220  public Expr< IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > > {
1221 
1222  public:
1223 
1224  typedef typename std::remove_cv<T1>::type ExprT1;
1225  typedef T2 ConstT;
1226  typedef typename ExprT1::value_type value_type;
1227  typedef typename ExprT1::scalar_type scalar_type;
1228 
1230 
1232  IfThenElseOp(const CondT& cond_, const T1& expr1_, const ConstT& c_) :
1233  cond(cond_), expr1(expr1_), c(c_) {}
1234 
1236  int size() const {
1237  return expr1.size();
1238  }
1239 
1241  bool hasFastAccess() const {
1242  return expr1.hasFastAccess();
1243  }
1244 
1246  value_type val() const {
1247  return if_then_else( cond, expr1.val(), c );
1248  }
1249 
1251  value_type dx(int i) const {
1252  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
1253  }
1254 
1256  value_type fastAccessDx(int i) const {
1257  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
1258  }
1259 
1260  protected:
1261 
1262  const CondT& cond;
1263  const T1& expr1;
1264  const ConstT& c;
1265  };
1266 
1267  template <typename CondT, typename T1, typename T2>
1268  class IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > :
1269  public Expr< IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > > {
1270 
1271  public:
1272 
1273  typedef typename std::remove_cv<T2>::type ExprT2;
1274  typedef T1 ConstT;
1275  typedef typename ExprT2::value_type value_type;
1276  typedef typename ExprT2::scalar_type scalar_type;
1277 
1279 
1281  IfThenElseOp(const CondT& cond_, const ConstT& c_, const T2& expr2_) :
1282  cond(cond_), c(c_), expr2(expr2_) {}
1283 
1285  int size() const {
1286  return expr2.size();
1287  }
1288 
1290  bool hasFastAccess() const {
1291  return expr2.hasFastAccess();
1292  }
1293 
1295  value_type val() const {
1296  return if_then_else( cond, c, expr2.val() );
1297  }
1298 
1300  value_type dx(int i) const {
1301  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
1302  }
1303 
1305  value_type fastAccessDx(int i) const {
1306  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
1307  }
1308 
1309  protected:
1310 
1311  const CondT& cond;
1312  const ConstT& c;
1313  const T2& expr2;
1314  };
1315 
1316  template <typename CondT, typename T1, typename T2>
1320  IfThenElseOp< CondT,
1321  typename Expr<T1>::derived_type,
1322  typename Expr<T2>::derived_type,
1323  false, false,
1324  typename T1::expr_spec_type >
1325  >::type
1326  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
1327  {
1329  typename Expr<T2>::derived_type,
1330  false, false, typename T1::expr_spec_type > expr_t;
1331 
1332  return expr_t(cond, expr1.derived(), expr2.derived());
1333  }
1334 
1335  template <typename CondT, typename T>
1337  IfThenElseOp< CondT, typename T::value_type, typename Expr<T>::derived_type,
1338  true, false, typename T::expr_spec_type >
1339  if_then_else (const CondT& cond, const typename T::value_type& c,
1340  const Expr<T>& expr)
1341  {
1342  typedef typename T::value_type ConstT;
1344  true, false, typename T::expr_spec_type > expr_t;
1345 
1346  return expr_t(cond, c, expr.derived());
1347  }
1348 
1349  template <typename CondT, typename T>
1351  IfThenElseOp< CondT, typename Expr<T>::derived_type, typename T::value_type,
1352  false, true, typename T::expr_spec_type >
1353  if_then_else (const CondT& cond, const Expr<T>& expr,
1354  const typename T::value_type& c)
1355  {
1356  typedef typename T::value_type ConstT;
1358  false, true, typename T::expr_spec_type > expr_t;
1359 
1360  return expr_t(cond, expr.derived(), c);
1361  }
1362 
1363  template <typename CondT, typename T>
1365  typename mpl::disable_if<
1366  mpl::is_same< typename T::value_type,
1367  typename T::scalar_type >,
1368  IfThenElseOp< CondT, typename T::scalar_type,
1369  typename Expr<T>::derived_type,
1370  true, false, typename T::expr_spec_type >
1371  >::type
1372  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
1373  const Expr<T>& expr)
1374  {
1375  typedef typename T::scalar_type ConstT;
1377  true, false, typename T::expr_spec_type > expr_t;
1378 
1379  return expr_t(cond, c, expr.derived());
1380  }
1381 
1382  template <typename CondT, typename T>
1384  typename mpl::disable_if<
1385  mpl::is_same< typename T::value_type,
1386  typename T::scalar_type >,
1387  IfThenElseOp< CondT, typename Expr<T>::derived_type,
1388  typename T::scalar_type,
1389  false, true, typename T::expr_spec_type >
1390  >::type
1391  if_then_else (const CondT& cond, const Expr<T>& expr,
1392  const typename Expr<T>::scalar_type& c)
1393  {
1394  typedef typename T::scalar_type ConstT;
1396  false, true, typename T::expr_spec_type > expr_t;
1397 
1398  return expr_t(cond, expr.derived(), c);
1399  }
1400 
1401  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1402  typename E>
1403  struct ExprLevel< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1404  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1405  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1406  static constexpr unsigned value =
1407  value_1 >= value_2 ? value_1 : value_2;
1408  };
1409 
1410  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1411  typename E>
1412  struct IsFadExpr< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1413  static constexpr unsigned value = true;
1414  };
1415 
1416  }
1417  }
1418 
1419  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1420  typename E>
1421  struct IsExpr< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1422  static constexpr bool value = true;
1423  };
1424 
1425  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1426  typename E>
1427  struct BaseExprType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1430  typedef typename Sacado::Promote<base_expr_1,
1432  };
1433 }
1434 
1435 #undef FAD_BINARYOP_MACRO
1436 
1437 //-------------------------- Relational Operators -----------------------
1438 
1439 namespace Sacado {
1440  namespace Fad {
1441  namespace Exp {
1442  template <typename T1, typename T2 = T1>
1444  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
1445  };
1446  }
1447  }
1448 }
1449 
1450 #define FAD_RELOP_MACRO(OP) \
1451 namespace Sacado { \
1452  namespace Fad { \
1453  namespace Exp { \
1454  template <typename T1, typename T2> \
1455  KOKKOS_INLINE_FUNCTION \
1456  typename mpl::enable_if_c< \
1457  IsFadExpr<T1>::value && IsFadExpr<T2>::value && \
1458  ExprLevel<T1>::value == ExprLevel<T2>::value, \
1459  typename ConditionalReturnType<typename T1::value_type, \
1460  typename T2::value_type>::type \
1461  >::type \
1462  operator OP (const T1& expr1, const T2& expr2) \
1463  { \
1464  return expr1.derived().val() OP expr2.derived().val(); \
1465  } \
1466  \
1467  template <typename T2> \
1468  KOKKOS_INLINE_FUNCTION \
1469  typename ConditionalReturnType<typename T2::value_type>::type \
1470  operator OP (const typename T2::value_type& a, \
1471  const Expr<T2>& expr2) \
1472  { \
1473  return a OP expr2.derived().val(); \
1474  } \
1475  \
1476  template <typename T1> \
1477  KOKKOS_INLINE_FUNCTION \
1478  typename ConditionalReturnType<typename T1::value_type>::type \
1479  operator OP (const Expr<T1>& expr1, \
1480  const typename T1::value_type& b) \
1481  { \
1482  return expr1.derived().val() OP b; \
1483  } \
1484  } \
1485  } \
1486 } \
1487 
1488 FAD_RELOP_MACRO(==)
1489 FAD_RELOP_MACRO(!=)
1490 FAD_RELOP_MACRO(<)
1491 FAD_RELOP_MACRO(>)
1492 FAD_RELOP_MACRO(<=)
1493 FAD_RELOP_MACRO(>=)
1494 FAD_RELOP_MACRO(<<=)
1495 FAD_RELOP_MACRO(>>=)
1496 FAD_RELOP_MACRO(&)
1497 FAD_RELOP_MACRO(|)
1498 
1499 #undef FAD_RELOP_MACRO
1500 
1501 namespace Sacado {
1502 
1503  namespace Fad {
1504  namespace Exp {
1505 
1506  template <typename ExprT>
1508  bool operator ! (const Expr<ExprT>& expr)
1509  {
1510  return ! expr.derived().val();
1511  }
1512 
1513  } // namespace Exp
1514  } // namespace Fad
1515 
1516 } // namespace Sacado
1517 
1518 //-------------------------- Boolean Operators -----------------------
1519 namespace Sacado {
1520 
1521  namespace Fad {
1522  namespace Exp {
1523 
1524  template <typename T>
1526  bool toBool(const Expr<T>& xx) {
1527  const typename Expr<T>::derived_type& x = xx.derived();
1528  bool is_zero = (x.val() == 0.0);
1529  for (int i=0; i<x.size(); i++)
1530  is_zero = is_zero && (x.dx(i) == 0.0);
1531  return !is_zero;
1532  }
1533 
1534  } // namespace Exp
1535  } // namespace Fad
1536 
1537 } // namespace Sacado
1538 
1539 #define FAD_BOOL_MACRO(OP) \
1540 namespace Sacado { \
1541  namespace Fad { \
1542  namespace Exp { \
1543  template <typename T1, typename T2> \
1544  KOKKOS_INLINE_FUNCTION \
1545  bool \
1546  operator OP (const Expr<T1>& expr1, \
1547  const Expr<T2>& expr2) \
1548  { \
1549  return toBool(expr1) OP toBool(expr2); \
1550  } \
1551  \
1552  template <typename T2> \
1553  KOKKOS_INLINE_FUNCTION \
1554  bool \
1555  operator OP (const typename Expr<T2>::value_type& a, \
1556  const Expr<T2>& expr2) \
1557  { \
1558  return a OP toBool(expr2); \
1559  } \
1560  \
1561  template <typename T1> \
1562  KOKKOS_INLINE_FUNCTION \
1563  bool \
1564  operator OP (const Expr<T1>& expr1, \
1565  const typename Expr<T1>::value_type& b) \
1566  { \
1567  return toBool(expr1) OP b; \
1568  } \
1569  } \
1570  } \
1571 }
1572 
1573 FAD_BOOL_MACRO(&&)
1574 FAD_BOOL_MACRO(||)
1575 
1576 #undef FAD_BOOL_MACRO
1577 
1578 //-------------------------- I/O Operators -----------------------
1579 
1580 namespace Sacado {
1581 
1582  namespace Fad {
1583  namespace Exp {
1584 
1585  template <typename T>
1586  std::ostream& operator << (std::ostream& os, const Expr<T>& xx) {
1587  const typename Expr<T>::derived_type& x = xx.derived();
1588  os << x.val() << " [";
1589 
1590  for (int i=0; i< x.size(); i++) {
1591  os << " " << x.dx(i);
1592  }
1593 
1594  os << " ]";
1595  return os;
1596  }
1597 
1598  } // namespace Exp
1599  } // namespace Fad
1600 
1601 } // namespace Sacado
1602 
1603 #if defined(HAVE_SACADO_KOKKOSCORE)
1604 
1605 //-------------------------- Atomic Operators -----------------------
1606 
1607 namespace Sacado {
1608 
1609  namespace Fad {
1610  namespace Exp {
1611 
1612  // Overload of Kokkos::atomic_add for Fad types.
1613  template <typename S>
1615  void atomic_add(GeneralFad<S>* dst, const GeneralFad<S>& x) {
1616  using Kokkos::atomic_add;
1617 
1618  const int xsz = x.size();
1619  const int sz = dst->size();
1620 
1621  // We currently cannot handle resizing since that would need to be
1622  // done atomically.
1623  if (xsz > sz)
1624  Kokkos::abort(
1625  "Sacado error: Fad resize within atomic_add() not supported!");
1626 
1627  if (xsz != sz && sz > 0 && xsz > 0)
1628  Kokkos::abort(
1629  "Sacado error: Fad assignment of incompatiable sizes!");
1630 
1631 
1632  if (sz > 0 && xsz > 0) {
1633  SACADO_FAD_DERIV_LOOP(i,sz)
1634  atomic_add(&(dst->fastAccessDx(i)), x.fastAccessDx(i));
1635  }
1637  atomic_add(&(dst->val()), x.val());
1638  }
1639 
1640  } // namespace Exp
1641  } // namespace Fad
1642 
1643 } // namespace Sacado
1644 
1645 #endif // HAVE_SACADO_KOKKOSCORE
1646 
1647 #endif // SACADO_FAD_OPS_HPP
Wrapper for a generic expression template.
cbrt(expr.val())
#define FAD_BOOL_MACRO(OP)
expr expr AdditionOp
expr expr SinOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
asinh(expr.val())
expr expr expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
#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())
decltype(std::declval< T1 >()==std::declval< T2 >()) typedef type
KOKKOS_INLINE_FUNCTION bool toBool(const Expr< T > &xx)
#define SACADO_FAD_THREAD_SINGLE
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MinOp
atanh(expr.val())
expr expr expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 Atan2Op
expr expr CoshOp
expr expr ATanhOp
KOKKOS_INLINE_FUNCTION const derived_type & derived() const
Return derived object.
expr expr TanhOp
expr expr SqrtOp
expr expr ASinhOp
Determine whether a given type is an expression.
expr true
KOKKOS_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)
KOKKOS_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const ConstT &c_)
static const bool value
atan(expr.val())
Wrapper for a generic expression template.
KOKKOS_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
expr expr expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 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
KOKKOS_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)
Is a type an expression.
expr val()
#define KOKKOS_INLINE_FUNCTION
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)
expr expr ATanOp
#define T2(r, f)
Definition: Sacado_rad.hpp:578
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())
#define T1(r, f)
Definition: Sacado_rad.hpp:603
Meta-function for determining nesting with an expression.
#define FAD_RELOP_MACRO(OP)
atan2(expr1.val(), expr2.val())
KOKKOS_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const T2 &expr2_)
#define SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP)
sin(expr.val())
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
log(expr.val())
expr expr ACoshOp
expr expr Log10Op
expr expr SinhOp
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
KOKKOS_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const ConstT &c_, const T2 &expr2_)
expr expr expr bar false
exp(expr.val())
expr expr expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 SubtractionOp
expr expr expr ExpOp
fabs(expr.val())
expr expr AbsOp
expr expr TanOp
log10(expr.val())
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