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_KOKKOSCORE)
45 #include "Kokkos_Atomic.hpp"
46 #include "impl/Kokkos_Error.hpp"
47 #endif
48 
49 #define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
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  KOKKOS_INLINE_FUNCTION \
69  OP(const T& expr_) : expr(expr_) {} \
70  \
71  KOKKOS_INLINE_FUNCTION \
72  int size() const { return expr.size(); } \
73  \
74  KOKKOS_INLINE_FUNCTION \
75  bool hasFastAccess() const { \
76  return expr.hasFastAccess(); \
77  } \
78  \
79  KOKKOS_INLINE_FUNCTION \
80  value_type val() const { \
81  USING \
82  return VALUE; \
83  } \
84  \
85  KOKKOS_INLINE_FUNCTION \
86  value_type dx(int i) const { \
87  USING \
88  return DX; \
89  } \
90  \
91  KOKKOS_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  KOKKOS_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 }
144 
145 FAD_UNARYOP_MACRO(operator+,
146  UnaryPlusOp,
147  ;,
148  expr.val(),
149  expr.dx(i),
150  expr.fastAccessDx(i))
151 FAD_UNARYOP_MACRO(operator-,
153  ;,
154  -expr.val(),
155  -expr.dx(i),
156  -expr.fastAccessDx(i))
159  using std::exp;,
160  exp(expr.val()),
161  exp(expr.val())*expr.dx(i),
162  exp(expr.val())*expr.fastAccessDx(i))
165  using std::log;,
166  log(expr.val()),
167  expr.dx(i)/expr.val(),
168  expr.fastAccessDx(i)/expr.val())
171  using std::log10; using std::log;,
172  log10(expr.val()),
173  expr.dx(i)/( log(value_type(10))*expr.val()),
174  expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
177  using std::sqrt;,
178  sqrt(expr.val()),
179  expr.dx(i)/(value_type(2)* sqrt(expr.val())),
180  expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
183  using std::cos; using std::sin;,
184  cos(expr.val()),
185  -expr.dx(i)* sin(expr.val()),
186  -expr.fastAccessDx(i)* sin(expr.val()))
189  using std::cos; using std::sin;,
190  sin(expr.val()),
191  expr.dx(i)* cos(expr.val()),
192  expr.fastAccessDx(i)* cos(expr.val()))
195  using std::tan;,
196  tan(expr.val()),
197  expr.dx(i)*
198  (value_type(1)+ tan(expr.val())* tan(expr.val())),
199  expr.fastAccessDx(i)*
200  (value_type(1)+ tan(expr.val())* tan(expr.val())))
203  using std::acos; using std::sqrt;,
204  acos(expr.val()),
205  -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
206  -expr.fastAccessDx(i) /
207  sqrt(value_type(1)-expr.val()*expr.val()))
210  using std::asin; using std::sqrt;,
211  asin(expr.val()),
212  expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
213  expr.fastAccessDx(i) /
214  sqrt(value_type(1)-expr.val()*expr.val()))
217  using std::atan;,
218  atan(expr.val()),
219  expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
220  expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
223  using std::cosh; using std::sinh;,
224  cosh(expr.val()),
225  expr.dx(i)* sinh(expr.val()),
226  expr.fastAccessDx(i)* sinh(expr.val()))
227 FAD_UNARYOP_MACRO(sinh,
229  using std::cosh; using std::sinh;,
230  sinh(expr.val()),
231  expr.dx(i)* cosh(expr.val()),
232  expr.fastAccessDx(i)* cosh(expr.val()))
235  using std::tanh; using std::cosh;,
236  tanh(expr.val()),
237  expr.dx(i)/( cosh(expr.val())* cosh(expr.val())),
238  expr.fastAccessDx(i) /
239  ( cosh(expr.val())* cosh(expr.val())))
242  using std::acosh; using std::sqrt;,
243  acosh(expr.val()),
244  expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
245  (expr.val()+value_type(1))),
246  expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
247  (expr.val()+value_type(1))))
250  using std::asinh; using std::sqrt;,
251  asinh(expr.val()),
252  expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
253  expr.fastAccessDx(i)/ sqrt(value_type(1)+
254  expr.val()*expr.val()))
257  using std::atanh;,
258  atanh(expr.val()),
259  expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
260  expr.fastAccessDx(i)/(value_type(1)-
261  expr.val()*expr.val()))
264  using std::abs;,
265  abs(expr.val()),
266  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
267  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
270  using std::fabs;,
271  fabs(expr.val()),
272  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
273  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
276  using std::cbrt;,
277  cbrt(expr.val()),
278  expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
279  expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
280 
281 #undef FAD_UNARYOP_MACRO
282 
283 #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) \
284 namespace Sacado { \
285  namespace Fad { \
286  namespace Exp { \
287  \
288  template <typename T1, typename T2, \
289  bool is_const_T1, bool is_const_T2, \
290  typename ExprSpec > \
291  class OP {}; \
292  \
293  template <typename T1, typename T2> \
294  class OP< T1, T2, false, false, ExprSpecDefault > : \
295  public Expr< OP< T1, T2, false, false, ExprSpecDefault > > { \
296  public: \
297  \
298  typedef typename std::remove_cv<T1>::type ExprT1; \
299  typedef typename std::remove_cv<T2>::type ExprT2; \
300  typedef typename ExprT1::value_type value_type_1; \
301  typedef typename ExprT2::value_type value_type_2; \
302  typedef typename Sacado::Promote<value_type_1, \
303  value_type_2>::type value_type; \
304  \
305  typedef typename ExprT1::scalar_type scalar_type_1; \
306  typedef typename ExprT2::scalar_type scalar_type_2; \
307  typedef typename Sacado::Promote<scalar_type_1, \
308  scalar_type_2>::type scalar_type; \
309  \
310  typedef ExprSpecDefault expr_spec_type; \
311  \
312  KOKKOS_INLINE_FUNCTION \
313  OP(const T1& expr1_, const T2& expr2_) : \
314  expr1(expr1_), expr2(expr2_) {} \
315  \
316  KOKKOS_INLINE_FUNCTION \
317  int size() const { \
318  const int sz1 = expr1.size(), sz2 = expr2.size(); \
319  return sz1 > sz2 ? sz1 : sz2; \
320  } \
321  \
322  KOKKOS_INLINE_FUNCTION \
323  bool hasFastAccess() const { \
324  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
325  } \
326  \
327  KOKKOS_INLINE_FUNCTION \
328  value_type val() const { \
329  USING \
330  return VALUE; \
331  } \
332  \
333  KOKKOS_INLINE_FUNCTION \
334  value_type dx(int i) const { \
335  USING \
336  const int sz1 = expr1.size(), sz2 = expr2.size(); \
337  if (sz1 > 0 && sz2 > 0) \
338  return DX; \
339  else if (sz1 > 0) \
340  return CDX2; \
341  else \
342  return CDX1; \
343  } \
344  \
345  KOKKOS_INLINE_FUNCTION \
346  value_type fastAccessDx(int i) const { \
347  USING \
348  return FASTACCESSDX; \
349  } \
350  \
351  protected: \
352  \
353  const T1& expr1; \
354  const T2& expr2; \
355  \
356  }; \
357  \
358  template <typename T1, typename T2> \
359  class OP< T1, T2, false, true, ExprSpecDefault > \
360  : public Expr< OP< T1, T2, false, true, ExprSpecDefault > > { \
361  public: \
362  \
363  typedef typename std::remove_cv<T1>::type ExprT1; \
364  typedef T2 ConstT; \
365  typedef typename ExprT1::value_type value_type; \
366  typedef typename ExprT1::scalar_type scalar_type; \
367  \
368  typedef ExprSpecDefault expr_spec_type; \
369  \
370  KOKKOS_INLINE_FUNCTION \
371  OP(const T1& expr1_, const ConstT& c_) : \
372  expr1(expr1_), c(c_) {} \
373  \
374  KOKKOS_INLINE_FUNCTION \
375  int size() const { \
376  return expr1.size(); \
377  } \
378  \
379  KOKKOS_INLINE_FUNCTION \
380  bool hasFastAccess() const { \
381  return expr1.hasFastAccess(); \
382  } \
383  \
384  KOKKOS_INLINE_FUNCTION \
385  value_type val() const { \
386  USING \
387  return VAL_CONST_DX_2; \
388  } \
389  \
390  KOKKOS_INLINE_FUNCTION \
391  value_type dx(int i) const { \
392  USING \
393  return CONST_DX_2; \
394  } \
395  \
396  KOKKOS_INLINE_FUNCTION \
397  value_type fastAccessDx(int i) const { \
398  USING \
399  return CONST_FASTACCESSDX_2; \
400  } \
401  \
402  protected: \
403  \
404  const T1& expr1; \
405  const ConstT& c; \
406  }; \
407  \
408  template <typename T1, typename T2> \
409  class OP< T1, T2, true, false, ExprSpecDefault > \
410  : public Expr< OP< T1, T2, true, false, ExprSpecDefault > > { \
411  public: \
412  \
413  typedef typename std::remove_cv<T2>::type ExprT2; \
414  typedef T1 ConstT; \
415  typedef typename ExprT2::value_type value_type; \
416  typedef typename ExprT2::scalar_type scalar_type; \
417  \
418  typedef ExprSpecDefault expr_spec_type; \
419  \
420  KOKKOS_INLINE_FUNCTION \
421  OP(const ConstT& c_, const T2& expr2_) : \
422  c(c_), expr2(expr2_) {} \
423  \
424  KOKKOS_INLINE_FUNCTION \
425  int size() const { \
426  return expr2.size(); \
427  } \
428  \
429  KOKKOS_INLINE_FUNCTION \
430  bool hasFastAccess() const { \
431  return expr2.hasFastAccess(); \
432  } \
433  \
434  KOKKOS_INLINE_FUNCTION \
435  value_type val() const { \
436  USING \
437  return VAL_CONST_DX_1; \
438  } \
439  \
440  KOKKOS_INLINE_FUNCTION \
441  value_type dx(int i) const { \
442  USING \
443  return CONST_DX_1; \
444  } \
445  \
446  KOKKOS_INLINE_FUNCTION \
447  value_type fastAccessDx(int i) const { \
448  USING \
449  return CONST_FASTACCESSDX_1; \
450  } \
451  \
452  protected: \
453  \
454  const ConstT& c; \
455  const T2& expr2; \
456  }; \
457  \
458  template <typename T1, typename T2> \
459  KOKKOS_INLINE_FUNCTION \
460  SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP) \
461  OPNAME (const T1& expr1, const T2& expr2) \
462  { \
463  typedef OP< typename Expr<T1>::derived_type, \
464  typename Expr<T2>::derived_type, \
465  false, false, typename T1::expr_spec_type > expr_t; \
466  \
467  return expr_t(expr1.derived(), expr2.derived()); \
468  } \
469  \
470  template <typename T> \
471  KOKKOS_INLINE_FUNCTION \
472  OP< typename T::value_type, typename Expr<T>::derived_type, \
473  true, false, typename T::expr_spec_type > \
474  OPNAME (const typename T::value_type& c, \
475  const Expr<T>& expr) \
476  { \
477  typedef typename T::value_type ConstT; \
478  typedef OP< ConstT, typename Expr<T>::derived_type, \
479  true, false, typename T::expr_spec_type > expr_t; \
480  \
481  return expr_t(c, expr.derived()); \
482  } \
483  \
484  template <typename T> \
485  KOKKOS_INLINE_FUNCTION \
486  OP< typename Expr<T>::derived_type, typename T::value_type, \
487  false, true, typename T::expr_spec_type > \
488  OPNAME (const Expr<T>& expr, \
489  const typename T::value_type& c) \
490  { \
491  typedef typename T::value_type ConstT; \
492  typedef OP< typename Expr<T>::derived_type, ConstT, \
493  false, true, typename T::expr_spec_type > expr_t; \
494  \
495  return expr_t(expr.derived(), c); \
496  } \
497  \
498  template <typename T> \
499  KOKKOS_INLINE_FUNCTION \
500  SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP) \
501  OPNAME (const typename T::scalar_type& c, \
502  const Expr<T>& expr) \
503  { \
504  typedef typename T::scalar_type ConstT; \
505  typedef OP< ConstT, typename Expr<T>::derived_type, \
506  true, false, typename T::expr_spec_type > expr_t; \
507  \
508  return expr_t(c, expr.derived()); \
509  } \
510  \
511  template <typename T> \
512  KOKKOS_INLINE_FUNCTION \
513  SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP) \
514  OPNAME (const Expr<T>& expr, \
515  const typename T::scalar_type& c) \
516  { \
517  typedef typename T::scalar_type ConstT; \
518  typedef OP< typename Expr<T>::derived_type, ConstT, \
519  false, true, typename T::expr_spec_type > expr_t; \
520  \
521  return expr_t(expr.derived(), c); \
522  } \
523  \
524  template <typename T1, typename T2, bool c1, bool c2, typename E> \
525  struct ExprLevel< OP< T1, T2, c1, c2, E > > { \
526  static constexpr unsigned value_1 = ExprLevel<T1>::value; \
527  static constexpr unsigned value_2 = ExprLevel<T2>::value; \
528  static constexpr unsigned value = \
529  value_1 >= value_2 ? value_1 : value_2; \
530  }; \
531  \
532  template <typename T1, typename T2, bool c1, bool c2, typename E> \
533  struct IsFadExpr< OP< T1, T2, c1, c2, E > > { \
534  static constexpr unsigned value = true; \
535  }; \
536  \
537  } \
538  } \
539  \
540  template <typename T1, typename T2, bool c1, bool c2, typename E> \
541  struct IsExpr< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
542  static constexpr bool value = true; \
543  }; \
544  \
545  template <typename T1, typename T2, bool c1, bool c2, typename E> \
546  struct BaseExprType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
547  typedef typename BaseExprType<T1>::type base_expr_1; \
548  typedef typename BaseExprType<T2>::type base_expr_2; \
549  typedef typename Sacado::Promote<base_expr_1, \
550  base_expr_2>::type type; \
551  }; \
552  \
553  template <typename T1, typename T2, bool c1, bool c2, typename E> \
554  struct IsSimdType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
555  static const bool value = \
556  IsSimdType< typename Fad::Exp::OP< T1, T2, c1, c2, E >::value_type >::value; \
557  }; \
558  \
559 }
560 
561 
562 FAD_BINARYOP_MACRO(operator+,
564  ;,
565  expr1.val() + expr2.val(),
566  expr1.dx(i) + expr2.dx(i),
567  expr2.dx(i),
568  expr1.dx(i),
569  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
570  c + expr2.val(),
571  expr1.val() + c,
572  expr2.dx(i),
573  expr1.dx(i),
574  expr2.fastAccessDx(i),
575  expr1.fastAccessDx(i))
576 FAD_BINARYOP_MACRO(operator-,
578  ;,
579  expr1.val() - expr2.val(),
580  expr1.dx(i) - expr2.dx(i),
581  -expr2.dx(i),
582  expr1.dx(i),
583  expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
584  c - expr2.val(),
585  expr1.val() - c,
586  -expr2.dx(i),
587  expr1.dx(i),
588  -expr2.fastAccessDx(i),
589  expr1.fastAccessDx(i))
590 FAD_BINARYOP_MACRO(operator*,
592  ;,
593  expr1.val() * expr2.val(),
594  expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
595  expr1.val()*expr2.dx(i),
596  expr1.dx(i)*expr2.val(),
597  expr1.val()*expr2.fastAccessDx(i) +
598  expr1.fastAccessDx(i)*expr2.val(),
599  c * expr2.val(),
600  expr1.val() * c,
601  c*expr2.dx(i),
602  expr1.dx(i)*c,
603  c*expr2.fastAccessDx(i),
604  expr1.fastAccessDx(i)*c)
605 FAD_BINARYOP_MACRO(operator/,
607  ;,
608  expr1.val() / expr2.val(),
609  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
610  (expr2.val()*expr2.val()),
611  -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
612  expr1.dx(i)/expr2.val(),
613  (expr1.fastAccessDx(i)*expr2.val() -
614  expr2.fastAccessDx(i)*expr1.val()) /
615  (expr2.val()*expr2.val()),
616  c / expr2.val(),
617  expr1.val() / c,
618  -expr2.dx(i)*c / (expr2.val()*expr2.val()),
619  expr1.dx(i)/c,
620  -expr2.fastAccessDx(i)*c / (expr2.val()*expr2.val()),
621  expr1.fastAccessDx(i)/c)
624  using std::atan2;,
625  atan2(expr1.val(), expr2.val()),
626  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
627  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
628  -expr1.val()*expr2.dx(i)/
629  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
630  expr2.val()*expr1.dx(i)/
631  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
632  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
633  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
634  atan2(c, expr2.val()),
635  atan2(expr1.val(), c),
636  (-c*expr2.dx(i)) / (c*c + expr2.val()*expr2.val()),
637  (c*expr1.dx(i))/ (expr1.val()*expr1.val() + c*c),
638  (-c*expr2.fastAccessDx(i))/ (c*c + expr2.val()*expr2.val()),
639  (c*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c*c))
640 // FAD_BINARYOP_MACRO(pow,
641 // PowerOp,
642 // using std::pow; using std::log;,
643 // pow(expr1.val(), expr2.val()),
644 // 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())) ),
645 // 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())) ),
646 // 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())) ),
647 // 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())) ),
648 // pow(c, expr2.val()),
649 // pow(expr1.val(), c),
650 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.dx(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.dx(i)/expr1.val()*pow(expr1.val(),c)) ),
652 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) ),
653 // 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))) )
656  ;,
657  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
658  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
659  if_then_else( expr1.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
660  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), value_type(0.0) ),
661  if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
662  if_then_else( c >= expr2.val(), value_type(c), expr2.val() ),
663  if_then_else( expr1.val() >= c, expr1.val(), value_type(c) ),
664  if_then_else( c >= expr2.val(), value_type(0.0), expr2.dx(i) ),
665  if_then_else( expr1.val() >= c, expr1.dx(i), value_type(0.0) ),
666  if_then_else( c >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
667  if_then_else( expr1.val() >= c, expr1.fastAccessDx(i), value_type(0.0) ) )
670  ;,
671  if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
672  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
673  if_then_else( expr1.val() <= expr2.val(), value_type(0.0), expr2.dx(i) ),
674  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), value_type(0.0) ),
675  if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
676  if_then_else( c <= expr2.val(), value_type(c), expr2.val() ),
677  if_then_else( expr1.val() <= c, expr1.val(), value_type(c) ),
678  if_then_else( c <= expr2.val(), value_type(0), expr2.dx(i) ),
679  if_then_else( expr1.val() <= c, expr1.dx(i), value_type(0) ),
680  if_then_else( c <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
681  if_then_else( expr1.val() <= c, expr1.fastAccessDx(i), value_type(0) ) )
682 
683 
684 // Special handling for std::pow() to provide specializations of PowerOp for
685 // "simd" value types that use if_then_else(). The only reason for not using
686 // if_then_else() always is to avoid evaluating the derivative if the value is
687 // zero to avoid throwing FPEs.
688 namespace Sacado {
689  namespace Fad {
690  namespace Exp {
691 
692  template <typename T1, typename T2,
693  bool is_const_T1, bool is_const_T2,
694  typename ExprSpec, bool is_simd >
695  class PowerOp {};
696 
697  //
698  // Implementation for simd type using if_then_else()
699  //
700  template <typename T1, typename T2>
701  class PowerOp< T1, T2, false, false, ExprSpecDefault, true > :
702  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault, true > > {
703  public:
704 
705  typedef typename std::remove_cv<T1>::type ExprT1;
706  typedef typename std::remove_cv<T2>::type ExprT2;
707  typedef typename ExprT1::value_type value_type_1;
708  typedef typename ExprT2::value_type value_type_2;
709  typedef typename Sacado::Promote<value_type_1,
710  value_type_2>::type value_type;
711 
712  typedef typename ExprT1::scalar_type scalar_type_1;
713  typedef typename ExprT2::scalar_type scalar_type_2;
714  typedef typename Sacado::Promote<scalar_type_1,
715  scalar_type_2>::type scalar_type;
716 
717  typedef ExprSpecDefault expr_spec_type;
718 
720  PowerOp(const T1& expr1_, const T2& expr2_) :
721  expr1(expr1_), expr2(expr2_) {}
722 
724  int size() const {
725  const int sz1 = expr1.size(), sz2 = expr2.size();
726  return sz1 > sz2 ? sz1 : sz2;
727  }
728 
730  bool hasFastAccess() const {
731  return expr1.hasFastAccess() && expr2.hasFastAccess();
732  }
733 
735  value_type val() const {
736  using std::pow;
737  return pow(expr1.val(), expr2.val());
738  }
739 
741  value_type dx(int i) const {
742  using std::pow; using std::log;
743  const int sz1 = expr1.size(), sz2 = expr2.size();
744  if (sz1 > 0 && sz2 > 0)
745  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())) );
746  else if (sz1 > 0)
747  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
748  // It seems less accurate and caused convergence problems in some codes
749  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())) );
750  else
751  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())) );
752  }
753 
755  value_type fastAccessDx(int i) const {
756  using std::pow; using std::log;
757  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())) );
758  }
759 
760  protected:
761 
762  const T1& expr1;
763  const T2& expr2;
764 
765  };
766 
767  template <typename T1, typename T2>
768  class PowerOp< T1, T2, false, true, ExprSpecDefault, true >
769  : public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault, true > > {
770  public:
771 
772  typedef typename std::remove_cv<T1>::type ExprT1;
773  typedef T2 ConstT;
774  typedef typename ExprT1::value_type value_type;
775  typedef typename ExprT1::scalar_type scalar_type;
776 
777  typedef ExprSpecDefault expr_spec_type;
778 
780  PowerOp(const T1& expr1_, const ConstT& c_) :
781  expr1(expr1_), c(c_) {}
782 
784  int size() const {
785  return expr1.size();
786  }
787 
789  bool hasFastAccess() const {
790  return expr1.hasFastAccess();
791  }
792 
794  value_type val() const {
795  using std::pow;
796  return pow(expr1.val(), c);
797  }
798 
800  value_type dx(int i) const {
801  using std::pow;
802  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
803  // It seems less accurate and caused convergence problems in some codes
804  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)) );
805  }
806 
808  value_type fastAccessDx(int i) const {
809  using std::pow;
810  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
811  // It seems less accurate and caused convergence problems in some codes
812  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)) );
813  }
814 
815  protected:
816 
817  const T1& expr1;
818  const ConstT& c;
819  };
820 
821  template <typename T1, typename T2>
822  class PowerOp< T1, T2, true, false, ExprSpecDefault, true >
823  : public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault, true > > {
824  public:
825 
826  typedef typename std::remove_cv<T2>::type ExprT2;
827  typedef T1 ConstT;
828  typedef typename ExprT2::value_type value_type;
829  typedef typename ExprT2::scalar_type scalar_type;
830 
831  typedef ExprSpecDefault expr_spec_type;
832 
834  PowerOp(const ConstT& c_, const T2& expr2_) :
835  c(c_), expr2(expr2_) {}
836 
838  int size() const {
839  return expr2.size();
840  }
841 
843  bool hasFastAccess() const {
844  return expr2.hasFastAccess();
845  }
846 
848  value_type val() const {
849  using std::pow;
850  return pow(c, expr2.val());
851  }
852 
854  value_type dx(int i) const {
855  using std::pow; using std::log;
856  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) );
857  }
858 
860  value_type fastAccessDx(int i) const {
861  using std::pow; using std::log;
862  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) );
863  }
864 
865  protected:
866 
867  const ConstT& c;
868  const T2& expr2;
869  };
870 
871  //
872  // Specialization for scalar types using ternary operator
873  //
874 
875  template <typename T1, typename T2>
876  class PowerOp< T1, T2, false, false, ExprSpecDefault, false > :
877  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault, false > > {
878  public:
879 
880  typedef typename std::remove_cv<T1>::type ExprT1;
881  typedef typename std::remove_cv<T2>::type ExprT2;
882  typedef typename ExprT1::value_type value_type_1;
883  typedef typename ExprT2::value_type value_type_2;
884  typedef typename Sacado::Promote<value_type_1,
885  value_type_2>::type value_type;
886 
887  typedef typename ExprT1::scalar_type scalar_type_1;
888  typedef typename ExprT2::scalar_type scalar_type_2;
889  typedef typename Sacado::Promote<scalar_type_1,
890  scalar_type_2>::type scalar_type;
891 
892  typedef ExprSpecDefault expr_spec_type;
893 
895  PowerOp(const T1& expr1_, const T2& expr2_) :
896  expr1(expr1_), expr2(expr2_) {}
897 
899  int size() const {
900  const int sz1 = expr1.size(), sz2 = expr2.size();
901  return sz1 > sz2 ? sz1 : sz2;
902  }
903 
905  bool hasFastAccess() const {
906  return expr1.hasFastAccess() && expr2.hasFastAccess();
907  }
908 
910  value_type val() const {
911  using std::pow;
912  return pow(expr1.val(), expr2.val());
913  }
914 
916  value_type dx(int i) const {
917  using std::pow; using std::log;
918  const int sz1 = expr1.size(), sz2 = expr2.size();
919  if (sz1 > 0 && sz2 > 0)
920  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()));
921  else if (sz1 > 0)
922  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
923  // It seems less accurate and caused convergence problems in some codes
924  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()));
925  else
926  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
927  }
928 
930  value_type fastAccessDx(int i) const {
931  using std::pow; using std::log;
932  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()));
933  }
934 
935  protected:
936 
937  const T1& expr1;
938  const T2& expr2;
939 
940  };
941 
942  template <typename T1, typename T2>
943  class PowerOp< T1, T2, false, true, ExprSpecDefault, false >
944  : public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault, false > > {
945  public:
946 
947  typedef typename std::remove_cv<T1>::type ExprT1;
948  typedef T2 ConstT;
949  typedef typename ExprT1::value_type value_type;
950  typedef typename ExprT1::scalar_type scalar_type;
951 
952  typedef ExprSpecDefault expr_spec_type;
953 
955  PowerOp(const T1& expr1_, const ConstT& c_) :
956  expr1(expr1_), c(c_) {}
957 
959  int size() const {
960  return expr1.size();
961  }
962 
964  bool hasFastAccess() const {
965  return expr1.hasFastAccess();
966  }
967 
969  value_type val() const {
970  using std::pow;
971  return pow(expr1.val(), c);
972  }
973 
975  value_type dx(int i) const {
976  using std::pow;
977  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
978  // It seems less accurate and caused convergence problems in some codes
979  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)/expr1.val()*pow(expr1.val(),c));
980  }
981 
983  value_type fastAccessDx(int i) const {
984  using std::pow;
985  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
986  // It seems less accurate and caused convergence problems in some codes
987  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c));
988  }
989 
990  protected:
991 
992  const T1& expr1;
993  const ConstT& c;
994  };
995 
996  template <typename T1, typename T2>
997  class PowerOp< T1, T2, true, false, ExprSpecDefault, false >
998  : public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault, false > > {
999  public:
1000 
1001  typedef typename std::remove_cv<T2>::type ExprT2;
1002  typedef T1 ConstT;
1003  typedef typename ExprT2::value_type value_type;
1004  typedef typename ExprT2::scalar_type scalar_type;
1005 
1006  typedef ExprSpecDefault expr_spec_type;
1007 
1009  PowerOp(const ConstT& c_, const T2& expr2_) :
1010  c(c_), expr2(expr2_) {}
1011 
1013  int size() const {
1014  return expr2.size();
1015  }
1016 
1018  bool hasFastAccess() const {
1019  return expr2.hasFastAccess();
1020  }
1021 
1023  value_type val() const {
1024  using std::pow;
1025  return pow(c, expr2.val());
1026  }
1027 
1029  value_type dx(int i) const {
1030  using std::pow; using std::log;
1031  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c)*pow(c,expr2.val()));
1032  }
1033 
1035  value_type fastAccessDx(int i) const {
1036  using std::pow; using std::log;
1037  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val()));
1038  }
1039 
1040  protected:
1041 
1042  const ConstT& c;
1043  const T2& expr2;
1044  };
1045 
1046  template <typename T1, typename T2>
1049  pow (const T1& expr1, const T2& expr2)
1050  {
1051  typedef PowerOp< typename Expr<T1>::derived_type,
1052  typename Expr<T2>::derived_type,
1053  false, false, typename T1::expr_spec_type > expr_t;
1054 
1055  return expr_t(expr1.derived(), expr2.derived());
1056  }
1057 
1058  template <typename T>
1060  PowerOp< typename T::value_type, typename Expr<T>::derived_type,
1061  true, false, typename T::expr_spec_type >
1062  pow (const typename T::value_type& c,
1063  const Expr<T>& expr)
1064  {
1065  typedef typename T::value_type ConstT;
1066  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1067  true, false, typename T::expr_spec_type > expr_t;
1068 
1069  return expr_t(c, expr.derived());
1070  }
1071 
1072  template <typename T>
1074  PowerOp< typename Expr<T>::derived_type, typename T::value_type,
1075  false, true, typename T::expr_spec_type >
1076  pow (const Expr<T>& expr,
1077  const typename T::value_type& c)
1078  {
1079  typedef typename T::value_type ConstT;
1080  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1081  false, true, typename T::expr_spec_type > expr_t;
1082 
1083  return expr_t(expr.derived(), c);
1084  }
1085 
1086  template <typename T>
1089  pow (const typename T::scalar_type& c,
1090  const Expr<T>& expr)
1091  {
1092  typedef typename T::scalar_type ConstT;
1093  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1094  true, false, typename T::expr_spec_type > expr_t;
1095 
1096  return expr_t(c, expr.derived());
1097  }
1098 
1099  template <typename T>
1102  pow (const Expr<T>& expr,
1103  const typename T::scalar_type& c)
1104  {
1105  typedef typename T::scalar_type ConstT;
1106  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1107  false, true, typename T::expr_spec_type > expr_t;
1108 
1109  return expr_t(expr.derived(), c);
1110  }
1111 
1112  template <typename T1, typename T2, bool c1, bool c2, typename E>
1113  struct ExprLevel< PowerOp< T1, T2, c1, c2, E > > {
1114  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1115  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1116  static constexpr unsigned value =
1117  value_1 >= value_2 ? value_1 : value_2;
1118  };
1119 
1120  template <typename T1, typename T2, bool c1, bool c2, typename E>
1121  struct IsFadExpr< PowerOp< T1, T2, c1, c2, E > > {
1122  static constexpr unsigned value = true;
1123  };
1124 
1125  }
1126  }
1127 
1128  template <typename T1, typename T2, bool c1, bool c2, typename E>
1129  struct IsExpr< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1130  static constexpr bool value = true;
1131  };
1132 
1133  template <typename T1, typename T2, bool c1, bool c2, typename E>
1134  struct BaseExprType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1135  typedef typename BaseExprType<T1>::type base_expr_1;
1136  typedef typename BaseExprType<T2>::type base_expr_2;
1137  typedef typename Sacado::Promote<base_expr_1,
1138  base_expr_2>::type type;
1139  };
1140 
1141  template <typename T1, typename T2, bool c1, bool c2, typename E>
1142  struct IsSimdType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1143  static const bool value =
1144  IsSimdType< typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::value_type >::value;
1145  };
1146 
1147 }
1148 
1149 //--------------------------if_then_else operator -----------------------
1150 // Can't use the above macros because it is a ternary operator (sort of).
1151 // Also, relies on C++11
1152 
1153 namespace Sacado {
1154  namespace Fad {
1155  namespace Exp {
1156 
1157  template <typename CondT, typename T1, typename T2,
1158  bool is_const_T1, bool is_const_T2,
1159  typename ExprSpec>
1160  class IfThenElseOp {};
1161 
1162  template <typename CondT, typename T1, typename T2>
1164  public Expr< IfThenElseOp< CondT, T1, T2, false, false, ExprSpecDefault > > {
1165 
1166  public:
1167 
1168  typedef typename std::remove_cv<T1>::type ExprT1;
1169  typedef typename std::remove_cv<T2>::type ExprT2;
1170  typedef typename ExprT1::value_type value_type_1;
1171  typedef typename ExprT2::value_type value_type_2;
1172  typedef typename Sacado::Promote<value_type_1,
1174 
1175  typedef typename ExprT1::scalar_type scalar_type_1;
1176  typedef typename ExprT2::scalar_type scalar_type_2;
1177  typedef typename Sacado::Promote<scalar_type_1,
1179 
1181 
1183  IfThenElseOp(const CondT& cond_, const T1& expr1_, const T2& expr2_) :
1184  cond(cond_), expr1(expr1_), expr2(expr2_) {}
1185 
1187  int size() const {
1188  int sz1 = expr1.size(), sz2 = expr2.size();
1189  return sz1 > sz2 ? sz1 : sz2;
1190  }
1191 
1193  bool hasFastAccess() const {
1194  return expr1.hasFastAccess() && expr2.hasFastAccess();
1195  }
1196 
1198  value_type val() const {
1199  return if_then_else( cond, expr1.val(), expr2.val() );
1200  }
1201 
1203  value_type dx(int i) const {
1204  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
1205  }
1206 
1208  value_type fastAccessDx(int i) const {
1209  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
1210  }
1211 
1212  protected:
1213 
1214  const CondT& cond;
1215  const T1& expr1;
1216  const T2& expr2;
1217 
1218  };
1219 
1220  template <typename CondT, typename T1, typename T2>
1221  class IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > :
1222  public Expr< IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > > {
1223 
1224  public:
1225 
1226  typedef typename std::remove_cv<T1>::type ExprT1;
1227  typedef T2 ConstT;
1228  typedef typename ExprT1::value_type value_type;
1229  typedef typename ExprT1::scalar_type scalar_type;
1230 
1232 
1234  IfThenElseOp(const CondT& cond_, const T1& expr1_, const ConstT& c_) :
1235  cond(cond_), expr1(expr1_), c(c_) {}
1236 
1238  int size() const {
1239  return expr1.size();
1240  }
1241 
1243  bool hasFastAccess() const {
1244  return expr1.hasFastAccess();
1245  }
1246 
1248  value_type val() const {
1249  return if_then_else( cond, expr1.val(), c );
1250  }
1251 
1253  value_type dx(int i) const {
1254  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
1255  }
1256 
1258  value_type fastAccessDx(int i) const {
1259  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
1260  }
1261 
1262  protected:
1263 
1264  const CondT& cond;
1265  const T1& expr1;
1266  const ConstT& c;
1267  };
1268 
1269  template <typename CondT, typename T1, typename T2>
1270  class IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > :
1271  public Expr< IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > > {
1272 
1273  public:
1274 
1275  typedef typename std::remove_cv<T2>::type ExprT2;
1276  typedef T1 ConstT;
1277  typedef typename ExprT2::value_type value_type;
1278  typedef typename ExprT2::scalar_type scalar_type;
1279 
1281 
1283  IfThenElseOp(const CondT& cond_, const ConstT& c_, const T2& expr2_) :
1284  cond(cond_), c(c_), expr2(expr2_) {}
1285 
1287  int size() const {
1288  return expr2.size();
1289  }
1290 
1292  bool hasFastAccess() const {
1293  return expr2.hasFastAccess();
1294  }
1295 
1297  value_type val() const {
1298  return if_then_else( cond, c, expr2.val() );
1299  }
1300 
1302  value_type dx(int i) const {
1303  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
1304  }
1305 
1307  value_type fastAccessDx(int i) const {
1308  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
1309  }
1310 
1311  protected:
1312 
1313  const CondT& cond;
1314  const ConstT& c;
1315  const T2& expr2;
1316  };
1317 
1318  template <typename CondT, typename T1, typename T2>
1322  IfThenElseOp< CondT,
1323  typename Expr<T1>::derived_type,
1324  typename Expr<T2>::derived_type,
1325  false, false,
1326  typename T1::expr_spec_type >
1327  >::type
1328  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
1329  {
1331  typename Expr<T2>::derived_type,
1332  false, false, typename T1::expr_spec_type > expr_t;
1333 
1334  return expr_t(cond, expr1.derived(), expr2.derived());
1335  }
1336 
1337  template <typename CondT, typename T>
1339  IfThenElseOp< CondT, typename T::value_type, typename Expr<T>::derived_type,
1340  true, false, typename T::expr_spec_type >
1341  if_then_else (const CondT& cond, const typename T::value_type& c,
1342  const Expr<T>& expr)
1343  {
1344  typedef typename T::value_type ConstT;
1346  true, false, typename T::expr_spec_type > expr_t;
1347 
1348  return expr_t(cond, c, expr.derived());
1349  }
1350 
1351  template <typename CondT, typename T>
1353  IfThenElseOp< CondT, typename Expr<T>::derived_type, typename T::value_type,
1354  false, true, typename T::expr_spec_type >
1355  if_then_else (const CondT& cond, const Expr<T>& expr,
1356  const typename T::value_type& c)
1357  {
1358  typedef typename T::value_type ConstT;
1360  false, true, typename T::expr_spec_type > expr_t;
1361 
1362  return expr_t(cond, expr.derived(), c);
1363  }
1364 
1365  template <typename CondT, typename T>
1367  typename mpl::disable_if<
1368  mpl::is_same< typename T::value_type,
1369  typename T::scalar_type >,
1370  IfThenElseOp< CondT, typename T::scalar_type,
1371  typename Expr<T>::derived_type,
1372  true, false, typename T::expr_spec_type >
1373  >::type
1374  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
1375  const Expr<T>& expr)
1376  {
1377  typedef typename T::scalar_type ConstT;
1379  true, false, typename T::expr_spec_type > expr_t;
1380 
1381  return expr_t(cond, c, expr.derived());
1382  }
1383 
1384  template <typename CondT, typename T>
1386  typename mpl::disable_if<
1387  mpl::is_same< typename T::value_type,
1388  typename T::scalar_type >,
1389  IfThenElseOp< CondT, typename Expr<T>::derived_type,
1390  typename T::scalar_type,
1391  false, true, typename T::expr_spec_type >
1392  >::type
1393  if_then_else (const CondT& cond, const Expr<T>& expr,
1394  const typename Expr<T>::scalar_type& c)
1395  {
1396  typedef typename T::scalar_type ConstT;
1398  false, true, typename T::expr_spec_type > expr_t;
1399 
1400  return expr_t(cond, expr.derived(), c);
1401  }
1402 
1403  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1404  typename E>
1405  struct ExprLevel< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1406  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1407  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1408  static constexpr unsigned value =
1409  value_1 >= value_2 ? value_1 : value_2;
1410  };
1411 
1412  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1413  typename E>
1414  struct IsFadExpr< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1415  static constexpr unsigned value = true;
1416  };
1417 
1418  }
1419  }
1420 
1421  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1422  typename E>
1423  struct IsExpr< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1424  static constexpr bool value = true;
1425  };
1426 
1427  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1428  typename E>
1429  struct BaseExprType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1432  typedef typename Sacado::Promote<base_expr_1,
1434  };
1435 }
1436 
1437 #undef FAD_BINARYOP_MACRO
1438 
1439 //-------------------------- Relational Operators -----------------------
1440 
1441 namespace Sacado {
1442  namespace Fad {
1443  namespace Exp {
1444  namespace Impl {
1445  // Helper trait to determine return type of logical comparison operations
1446  // (==, !=, ...), usually bool but maybe something else for SIMD types.
1447  // Need to use SFINAE to restrict to types that define == operator in the
1448  // conditional overloads below, otherwise instantiating ConditionaReturnType
1449  // may fail during overload resolution.
1450  template <typename T1, typename T2 = T1,
1453 
1454  template <typename T1, typename T2>
1456  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
1457  };
1458  }
1459  }
1460  }
1461 }
1462 
1463 #define FAD_RELOP_MACRO(OP) \
1464 namespace Sacado { \
1465  namespace Fad { \
1466  namespace Exp { \
1467  template <typename T1, typename T2> \
1468  KOKKOS_INLINE_FUNCTION \
1469  typename mpl::enable_if_c< \
1470  IsFadExpr<T1>::value && IsFadExpr<T2>::value && \
1471  ExprLevel<T1>::value == ExprLevel<T2>::value, \
1472  typename Impl::ConditionalReturnType<typename T1::value_type, \
1473  typename T2::value_type>::type \
1474  >::type \
1475  operator OP (const T1& expr1, const T2& expr2) \
1476  { \
1477  return expr1.derived().val() OP expr2.derived().val(); \
1478  } \
1479  \
1480  template <typename T2> \
1481  KOKKOS_INLINE_FUNCTION \
1482  typename Impl::ConditionalReturnType<typename T2::value_type>::type \
1483  operator OP (const typename T2::value_type& a, \
1484  const Expr<T2>& expr2) \
1485  { \
1486  return a OP expr2.derived().val(); \
1487  } \
1488  \
1489  template <typename T1> \
1490  KOKKOS_INLINE_FUNCTION \
1491  typename Impl::ConditionalReturnType<typename T1::value_type>::type \
1492  operator OP (const Expr<T1>& expr1, \
1493  const typename T1::value_type& b) \
1494  { \
1495  return expr1.derived().val() OP b; \
1496  } \
1497  } \
1498  } \
1499 } \
1500 
1501 FAD_RELOP_MACRO(==)
1502 FAD_RELOP_MACRO(!=)
1503 FAD_RELOP_MACRO(<)
1504 FAD_RELOP_MACRO(>)
1505 FAD_RELOP_MACRO(<=)
1506 FAD_RELOP_MACRO(>=)
1507 FAD_RELOP_MACRO(<<=)
1508 FAD_RELOP_MACRO(>>=)
1509 FAD_RELOP_MACRO(&)
1510 FAD_RELOP_MACRO(|)
1511 
1512 #undef FAD_RELOP_MACRO
1513 
1514 namespace Sacado {
1515 
1516  namespace Fad {
1517  namespace Exp {
1518 
1519  template <typename ExprT>
1521  bool operator ! (const Expr<ExprT>& expr)
1522  {
1523  return ! expr.derived().val();
1524  }
1525 
1526  } // namespace Exp
1527  } // namespace Fad
1528 
1529 } // namespace Sacado
1530 
1531 //-------------------------- Boolean Operators -----------------------
1532 namespace Sacado {
1533 
1534  namespace Fad {
1535  namespace Exp {
1536 
1537  template <typename T>
1539  bool toBool(const Expr<T>& xx) {
1540  const typename Expr<T>::derived_type& x = xx.derived();
1541  bool is_zero = (x.val() == 0.0);
1542  for (int i=0; i<x.size(); i++)
1543  is_zero = is_zero && (x.dx(i) == 0.0);
1544  return !is_zero;
1545  }
1546 
1547  } // namespace Exp
1548  } // namespace Fad
1549 
1550 } // namespace Sacado
1551 
1552 #define FAD_BOOL_MACRO(OP) \
1553 namespace Sacado { \
1554  namespace Fad { \
1555  namespace Exp { \
1556  template <typename T1, typename T2> \
1557  KOKKOS_INLINE_FUNCTION \
1558  bool \
1559  operator OP (const Expr<T1>& expr1, \
1560  const Expr<T2>& expr2) \
1561  { \
1562  return toBool(expr1) OP toBool(expr2); \
1563  } \
1564  \
1565  template <typename T2> \
1566  KOKKOS_INLINE_FUNCTION \
1567  bool \
1568  operator OP (const typename Expr<T2>::value_type& a, \
1569  const Expr<T2>& expr2) \
1570  { \
1571  return a OP toBool(expr2); \
1572  } \
1573  \
1574  template <typename T1> \
1575  KOKKOS_INLINE_FUNCTION \
1576  bool \
1577  operator OP (const Expr<T1>& expr1, \
1578  const typename Expr<T1>::value_type& b) \
1579  { \
1580  return toBool(expr1) OP b; \
1581  } \
1582  } \
1583  } \
1584 }
1585 
1586 FAD_BOOL_MACRO(&&)
1587 FAD_BOOL_MACRO(||)
1588 
1589 #undef FAD_BOOL_MACRO
1590 
1591 //-------------------------- I/O Operators -----------------------
1592 
1593 namespace Sacado {
1594 
1595  namespace Fad {
1596  namespace Exp {
1597 
1598  template <typename T>
1599  std::ostream& operator << (std::ostream& os, const Expr<T>& xx) {
1600  const typename Expr<T>::derived_type& x = xx.derived();
1601  os << x.val() << " [";
1602 
1603  for (int i=0; i< x.size(); i++) {
1604  os << " " << x.dx(i);
1605  }
1606 
1607  os << " ]";
1608  return os;
1609  }
1610 
1611  } // namespace Exp
1612  } // namespace Fad
1613 
1614 } // namespace Sacado
1615 
1616 #if defined(HAVE_SACADO_KOKKOSCORE)
1617 
1618 //-------------------------- Atomic Operators -----------------------
1619 
1620 namespace Sacado {
1621 
1622  namespace Fad {
1623  namespace Exp {
1624 
1625  // Overload of Kokkos::atomic_add for Fad types.
1626  template <typename S>
1628  void atomic_add(GeneralFad<S>* dst, const GeneralFad<S>& x) {
1629  using Kokkos::atomic_add;
1630 
1631  const int xsz = x.size();
1632  const int sz = dst->size();
1633 
1634  // We currently cannot handle resizing since that would need to be
1635  // done atomically.
1636  if (xsz > sz)
1637  Kokkos::abort(
1638  "Sacado error: Fad resize within atomic_add() not supported!");
1639 
1640  if (xsz != sz && sz > 0 && xsz > 0)
1641  Kokkos::abort(
1642  "Sacado error: Fad assignment of incompatiable sizes!");
1643 
1644 
1645  if (sz > 0 && xsz > 0) {
1646  SACADO_FAD_DERIV_LOOP(i,sz)
1647  atomic_add(&(dst->fastAccessDx(i)), x.fastAccessDx(i));
1648  }
1650  atomic_add(&(dst->val()), x.val());
1651  }
1652 
1653  } // namespace Exp
1654  } // namespace Fad
1655 
1656 } // namespace Sacado
1657 
1658 #endif // HAVE_SACADO_KOKKOSCORE
1659 
1660 #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())
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