Sacado Package Browser (Single Doxygen Collection)  Version of the Day
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 
38 #include "Sacado_cmath.hpp"
39 
40 #if defined(HAVE_SACADO_KOKKOSCORE)
41 #include "Kokkos_Atomic.hpp"
42 #include "impl/Kokkos_Error.hpp"
43 #endif
44 
45 #define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
46 namespace Sacado { \
47  namespace Fad { \
48  namespace Exp { \
49  \
50  template <typename T, typename ExprSpec> \
51  class OP {}; \
52  \
53  template <typename T> \
54  class OP< T,ExprSpecDefault > : \
55  public Expr< OP<T,ExprSpecDefault> > { \
56  public: \
57  \
58  typedef typename std::remove_cv<T>::type ExprT; \
59  typedef typename ExprT::value_type value_type; \
60  typedef typename ExprT::scalar_type scalar_type; \
61  \
62  typedef ExprSpecDefault expr_spec_type; \
63  \
64  KOKKOS_INLINE_FUNCTION \
65  OP(const T& expr_) : expr(expr_) {} \
66  \
67  KOKKOS_INLINE_FUNCTION \
68  int size() const { return expr.size(); } \
69  \
70  KOKKOS_INLINE_FUNCTION \
71  bool hasFastAccess() const { \
72  return expr.hasFastAccess(); \
73  } \
74  \
75  KOKKOS_INLINE_FUNCTION \
76  value_type val() const { \
77  USING \
78  return VALUE; \
79  } \
80  \
81  KOKKOS_INLINE_FUNCTION \
82  value_type dx(int i) const { \
83  USING \
84  return DX; \
85  } \
86  \
87  KOKKOS_INLINE_FUNCTION \
88  value_type fastAccessDx(int i) const { \
89  USING \
90  return FASTACCESSDX; \
91  } \
92  \
93  protected: \
94  \
95  const T& expr; \
96  }; \
97  \
98  template <typename T> \
99  KOKKOS_INLINE_FUNCTION \
100  OP< typename Expr<T>::derived_type, \
101  typename T::expr_spec_type > \
102  OPNAME (const Expr<T>& expr) \
103  { \
104  typedef OP< typename Expr<T>::derived_type, \
105  typename T::expr_spec_type > expr_t; \
106  \
107  return expr_t(expr.derived()); \
108  } \
109  \
110  template <typename T, typename E> \
111  struct ExprLevel< OP< T,E > > { \
112  static const unsigned value = ExprLevel<T>::value; \
113  }; \
114  \
115  template <typename T, typename E> \
116  struct IsFadExpr< OP< T,E > > { \
117  static const unsigned value = true; \
118  }; \
119  \
120  } \
121  } \
122  \
123  template <typename T, typename E> \
124  struct IsExpr< Fad::Exp::OP< T,E > > { \
125  static const bool value = true; \
126  }; \
127  \
128  template <typename T, typename E> \
129  struct BaseExprType< Fad::Exp::OP< T,E > > { \
130  typedef typename BaseExprType<T>::type type; \
131  }; \
132  \
133 }
134 
135 FAD_UNARYOP_MACRO(operator+,
136  UnaryPlusOp,
137  ;,
138  expr.val(),
139  expr.dx(i),
140  expr.fastAccessDx(i))
141 FAD_UNARYOP_MACRO(operator-,
143  ;,
144  -expr.val(),
145  -expr.dx(i),
146  -expr.fastAccessDx(i))
149  using std::exp;,
150  exp(expr.val()),
151  exp(expr.val())*expr.dx(i),
152  exp(expr.val())*expr.fastAccessDx(i))
155  using std::log;,
156  log(expr.val()),
157  expr.dx(i)/expr.val(),
158  expr.fastAccessDx(i)/expr.val())
161  using std::log10; using std::log;,
162  log10(expr.val()),
163  expr.dx(i)/( log(value_type(10))*expr.val()),
164  expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
167  using std::sqrt;,
168  sqrt(expr.val()),
169  expr.dx(i)/(value_type(2)* sqrt(expr.val())),
170  expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
173  using std::cos; using std::sin;,
174  cos(expr.val()),
175  -expr.dx(i)* sin(expr.val()),
176  -expr.fastAccessDx(i)* sin(expr.val()))
179  using std::cos; using std::sin;,
180  sin(expr.val()),
181  expr.dx(i)* cos(expr.val()),
182  expr.fastAccessDx(i)* cos(expr.val()))
185  using std::tan;,
186  tan(expr.val()),
187  expr.dx(i)*
188  (value_type(1)+ tan(expr.val())* tan(expr.val())),
189  expr.fastAccessDx(i)*
190  (value_type(1)+ tan(expr.val())* tan(expr.val())))
193  using std::acos; using std::sqrt;,
194  acos(expr.val()),
195  -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
196  -expr.fastAccessDx(i) /
197  sqrt(value_type(1)-expr.val()*expr.val()))
200  using std::asin; using std::sqrt;,
201  asin(expr.val()),
202  expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
203  expr.fastAccessDx(i) /
204  sqrt(value_type(1)-expr.val()*expr.val()))
207  using std::atan;,
208  atan(expr.val()),
209  expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
210  expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
213  using std::cosh; using std::sinh;,
214  cosh(expr.val()),
215  expr.dx(i)* sinh(expr.val()),
216  expr.fastAccessDx(i)* sinh(expr.val()))
219  using std::cosh; using std::sinh;,
220  sinh(expr.val()),
221  expr.dx(i)* cosh(expr.val()),
222  expr.fastAccessDx(i)* cosh(expr.val()))
225  using std::tanh; using std::cosh;,
226  tanh(expr.val()),
227  expr.dx(i)/( cosh(expr.val())* cosh(expr.val())),
228  expr.fastAccessDx(i) /
229  ( cosh(expr.val())* cosh(expr.val())))
232  using std::acosh; using std::sqrt;,
233  acosh(expr.val()),
234  expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
235  (expr.val()+value_type(1))),
236  expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
237  (expr.val()+value_type(1))))
240  using std::asinh; using std::sqrt;,
241  asinh(expr.val()),
242  expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
243  expr.fastAccessDx(i)/ sqrt(value_type(1)+
244  expr.val()*expr.val()))
247  using std::atanh;,
248  atanh(expr.val()),
249  expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
250  expr.fastAccessDx(i)/(value_type(1)-
251  expr.val()*expr.val()))
254  using std::abs;,
255  abs(expr.val()),
256  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
257  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
260  using std::fabs;,
261  fabs(expr.val()),
262  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
263  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
266  using std::cbrt;,
267  cbrt(expr.val()),
268  expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
269  expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
270 
271 #undef FAD_UNARYOP_MACRO
272 
273 #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) \
274 namespace Sacado { \
275  namespace Fad { \
276  namespace Exp { \
277  \
278  template <typename T1, typename T2, \
279  bool is_const_T1, bool is_const_T2, \
280  typename ExprSpec > \
281  class OP {}; \
282  \
283  template <typename T1, typename T2> \
284  class OP< T1, T2, false, false, ExprSpecDefault > : \
285  public Expr< OP< T1, T2, false, false, ExprSpecDefault > > { \
286  public: \
287  \
288  typedef typename std::remove_cv<T1>::type ExprT1; \
289  typedef typename std::remove_cv<T2>::type ExprT2; \
290  typedef typename ExprT1::value_type value_type_1; \
291  typedef typename ExprT2::value_type value_type_2; \
292  typedef typename Sacado::Promote<value_type_1, \
293  value_type_2>::type value_type; \
294  \
295  typedef typename ExprT1::scalar_type scalar_type_1; \
296  typedef typename ExprT2::scalar_type scalar_type_2; \
297  typedef typename Sacado::Promote<scalar_type_1, \
298  scalar_type_2>::type scalar_type; \
299  \
300  typedef ExprSpecDefault expr_spec_type; \
301  \
302  KOKKOS_INLINE_FUNCTION \
303  OP(const T1& expr1_, const T2& expr2_) : \
304  expr1(expr1_), expr2(expr2_) {} \
305  \
306  KOKKOS_INLINE_FUNCTION \
307  int size() const { \
308  const int sz1 = expr1.size(), sz2 = expr2.size(); \
309  return sz1 > sz2 ? sz1 : sz2; \
310  } \
311  \
312  KOKKOS_INLINE_FUNCTION \
313  bool hasFastAccess() const { \
314  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
315  } \
316  \
317  KOKKOS_INLINE_FUNCTION \
318  value_type val() const { \
319  USING \
320  return VALUE; \
321  } \
322  \
323  KOKKOS_INLINE_FUNCTION \
324  value_type dx(int i) const { \
325  USING \
326  const int sz1 = expr1.size(), sz2 = expr2.size(); \
327  if (sz1 > 0 && sz2 > 0) \
328  return DX; \
329  else if (sz1 > 0) \
330  return CDX2; \
331  else \
332  return CDX1; \
333  } \
334  \
335  KOKKOS_INLINE_FUNCTION \
336  value_type fastAccessDx(int i) const { \
337  USING \
338  return FASTACCESSDX; \
339  } \
340  \
341  protected: \
342  \
343  const T1& expr1; \
344  const T2& expr2; \
345  \
346  }; \
347  \
348  template <typename T1, typename T2> \
349  class OP< T1, T2, false, true, ExprSpecDefault > \
350  : public Expr< OP< T1, T2, false, true, ExprSpecDefault > > { \
351  public: \
352  \
353  typedef typename std::remove_cv<T1>::type ExprT1; \
354  typedef T2 ConstT; \
355  typedef typename ExprT1::value_type value_type; \
356  typedef typename ExprT1::scalar_type scalar_type; \
357  \
358  typedef ExprSpecDefault expr_spec_type; \
359  \
360  KOKKOS_INLINE_FUNCTION \
361  OP(const T1& expr1_, const ConstT& c_) : \
362  expr1(expr1_), c(c_) {} \
363  \
364  KOKKOS_INLINE_FUNCTION \
365  int size() const { \
366  return expr1.size(); \
367  } \
368  \
369  KOKKOS_INLINE_FUNCTION \
370  bool hasFastAccess() const { \
371  return expr1.hasFastAccess(); \
372  } \
373  \
374  KOKKOS_INLINE_FUNCTION \
375  value_type val() const { \
376  USING \
377  return VAL_CONST_DX_2; \
378  } \
379  \
380  KOKKOS_INLINE_FUNCTION \
381  value_type dx(int i) const { \
382  USING \
383  return CONST_DX_2; \
384  } \
385  \
386  KOKKOS_INLINE_FUNCTION \
387  value_type fastAccessDx(int i) const { \
388  USING \
389  return CONST_FASTACCESSDX_2; \
390  } \
391  \
392  protected: \
393  \
394  const T1& expr1; \
395  const ConstT& c; \
396  }; \
397  \
398  template <typename T1, typename T2> \
399  class OP< T1, T2, true, false, ExprSpecDefault > \
400  : public Expr< OP< T1, T2, true, false, ExprSpecDefault > > { \
401  public: \
402  \
403  typedef typename std::remove_cv<T2>::type ExprT2; \
404  typedef T1 ConstT; \
405  typedef typename ExprT2::value_type value_type; \
406  typedef typename ExprT2::scalar_type scalar_type; \
407  \
408  typedef ExprSpecDefault expr_spec_type; \
409  \
410  KOKKOS_INLINE_FUNCTION \
411  OP(const ConstT& c_, const T2& expr2_) : \
412  c(c_), expr2(expr2_) {} \
413  \
414  KOKKOS_INLINE_FUNCTION \
415  int size() const { \
416  return expr2.size(); \
417  } \
418  \
419  KOKKOS_INLINE_FUNCTION \
420  bool hasFastAccess() const { \
421  return expr2.hasFastAccess(); \
422  } \
423  \
424  KOKKOS_INLINE_FUNCTION \
425  value_type val() const { \
426  USING \
427  return VAL_CONST_DX_1; \
428  } \
429  \
430  KOKKOS_INLINE_FUNCTION \
431  value_type dx(int i) const { \
432  USING \
433  return CONST_DX_1; \
434  } \
435  \
436  KOKKOS_INLINE_FUNCTION \
437  value_type fastAccessDx(int i) const { \
438  USING \
439  return CONST_FASTACCESSDX_1; \
440  } \
441  \
442  protected: \
443  \
444  const ConstT& c; \
445  const T2& expr2; \
446  }; \
447  \
448  template <typename T1, typename T2> \
449  KOKKOS_INLINE_FUNCTION \
450  SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP) \
451  OPNAME (const T1& expr1, const T2& expr2) \
452  { \
453  typedef OP< typename Expr<T1>::derived_type, \
454  typename Expr<T2>::derived_type, \
455  false, false, typename T1::expr_spec_type > expr_t; \
456  \
457  return expr_t(expr1.derived(), expr2.derived()); \
458  } \
459  \
460  template <typename T> \
461  KOKKOS_INLINE_FUNCTION \
462  OP< typename T::value_type, typename Expr<T>::derived_type, \
463  true, false, typename T::expr_spec_type > \
464  OPNAME (const typename T::value_type& c, \
465  const Expr<T>& expr) \
466  { \
467  typedef typename T::value_type ConstT; \
468  typedef OP< ConstT, typename Expr<T>::derived_type, \
469  true, false, typename T::expr_spec_type > expr_t; \
470  \
471  return expr_t(c, expr.derived()); \
472  } \
473  \
474  template <typename T> \
475  KOKKOS_INLINE_FUNCTION \
476  OP< typename Expr<T>::derived_type, typename T::value_type, \
477  false, true, typename T::expr_spec_type > \
478  OPNAME (const Expr<T>& expr, \
479  const typename T::value_type& c) \
480  { \
481  typedef typename T::value_type ConstT; \
482  typedef OP< typename Expr<T>::derived_type, ConstT, \
483  false, true, typename T::expr_spec_type > expr_t; \
484  \
485  return expr_t(expr.derived(), c); \
486  } \
487  \
488  template <typename T> \
489  KOKKOS_INLINE_FUNCTION \
490  SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP) \
491  OPNAME (const typename T::scalar_type& c, \
492  const Expr<T>& expr) \
493  { \
494  typedef typename T::scalar_type ConstT; \
495  typedef OP< ConstT, typename Expr<T>::derived_type, \
496  true, false, typename T::expr_spec_type > expr_t; \
497  \
498  return expr_t(c, expr.derived()); \
499  } \
500  \
501  template <typename T> \
502  KOKKOS_INLINE_FUNCTION \
503  SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP) \
504  OPNAME (const Expr<T>& expr, \
505  const typename T::scalar_type& c) \
506  { \
507  typedef typename T::scalar_type ConstT; \
508  typedef OP< typename Expr<T>::derived_type, ConstT, \
509  false, true, typename T::expr_spec_type > expr_t; \
510  \
511  return expr_t(expr.derived(), c); \
512  } \
513  \
514  template <typename T1, typename T2, bool c1, bool c2, typename E> \
515  struct ExprLevel< OP< T1, T2, c1, c2, E > > { \
516  static constexpr unsigned value_1 = ExprLevel<T1>::value; \
517  static constexpr unsigned value_2 = ExprLevel<T2>::value; \
518  static constexpr unsigned value = \
519  value_1 >= value_2 ? value_1 : value_2; \
520  }; \
521  \
522  template <typename T1, typename T2, bool c1, bool c2, typename E> \
523  struct IsFadExpr< OP< T1, T2, c1, c2, E > > { \
524  static constexpr unsigned value = true; \
525  }; \
526  \
527  } \
528  } \
529  \
530  template <typename T1, typename T2, bool c1, bool c2, typename E> \
531  struct IsExpr< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
532  static constexpr bool value = true; \
533  }; \
534  \
535  template <typename T1, typename T2, bool c1, bool c2, typename E> \
536  struct BaseExprType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
537  typedef typename BaseExprType<T1>::type base_expr_1; \
538  typedef typename BaseExprType<T2>::type base_expr_2; \
539  typedef typename Sacado::Promote<base_expr_1, \
540  base_expr_2>::type type; \
541  }; \
542  \
543 }
544 
545 
546 FAD_BINARYOP_MACRO(operator+,
548  ;,
549  expr1.val() + expr2.val(),
550  expr1.dx(i) + expr2.dx(i),
551  expr2.dx(i),
552  expr1.dx(i),
553  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
554  c + expr2.val(),
555  expr1.val() + c,
556  expr2.dx(i),
557  expr1.dx(i),
558  expr2.fastAccessDx(i),
559  expr1.fastAccessDx(i))
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.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
579  expr1.val()*expr2.dx(i),
580  expr1.dx(i)*expr2.val(),
581  expr1.val()*expr2.fastAccessDx(i) +
582  expr1.fastAccessDx(i)*expr2.val(),
583  c * expr2.val(),
584  expr1.val() * c,
585  c*expr2.dx(i),
586  expr1.dx(i)*c,
587  c*expr2.fastAccessDx(i),
588  expr1.fastAccessDx(i)*c)
589 FAD_BINARYOP_MACRO(operator/,
591  ;,
592  expr1.val() / expr2.val(),
593  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
594  (expr2.val()*expr2.val()),
595  -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
596  expr1.dx(i)/expr2.val(),
597  (expr1.fastAccessDx(i)*expr2.val() -
598  expr2.fastAccessDx(i)*expr1.val()) /
599  (expr2.val()*expr2.val()),
600  c / expr2.val(),
601  expr1.val() / c,
602  -expr2.dx(i)*c / (expr2.val()*expr2.val()),
603  expr1.dx(i)/c,
604  -expr2.fastAccessDx(i)*c / (expr2.val()*expr2.val()),
605  expr1.fastAccessDx(i)/c)
608  using std::atan2;,
609  atan2(expr1.val(), expr2.val()),
610  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
611  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
612  -expr1.val()*expr2.dx(i)/
613  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
614  expr2.val()*expr1.dx(i)/
615  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
616  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
617  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
618  atan2(c, expr2.val()),
619  atan2(expr1.val(), c),
620  (-c*expr2.dx(i)) / (c*c + expr2.val()*expr2.val()),
621  (c*expr1.dx(i))/ (expr1.val()*expr1.val() + c*c),
622  (-c*expr2.fastAccessDx(i))/ (c*c + expr2.val()*expr2.val()),
623  (c*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c*c))
626  using std::pow; using std::log;,
627  pow(expr1.val(), expr2.val()),
628  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())) ),
629  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())) ),
630  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())) ),
631  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())) ),
632  pow(c, expr2.val()),
633  pow(expr1.val(), c),
634  if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) ),
635  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)) ),
636  if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) ),
637  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))) )
640  ;,
641  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
642  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
643  if_then_else( expr1.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
644  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), value_type(0.0) ),
645  if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
646  if_then_else( c >= expr2.val(), value_type(c), expr2.val() ),
647  if_then_else( expr1.val() >= c, expr1.val(), value_type(c) ),
648  if_then_else( c >= expr2.val(), value_type(0.0), expr2.dx(i) ),
649  if_then_else( expr1.val() >= c, expr1.dx(i), value_type(0.0) ),
650  if_then_else( c >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
651  if_then_else( expr1.val() >= c, expr1.fastAccessDx(i), value_type(0.0) ) )
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), expr2.dx(i) ),
663  if_then_else( expr1.val() <= c, expr1.dx(i), value_type(0) ),
664  if_then_else( c <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
665  if_then_else( expr1.val() <= c, expr1.fastAccessDx(i), value_type(0) ) )
666 
667 //--------------------------if_then_else operator -----------------------
668 // Can't use the above macros because it is a ternary operator (sort of).
669 // Also, relies on C++11
670 
671 namespace Sacado {
672  namespace Fad {
673  namespace Exp {
674 
675  template <typename CondT, typename T1, typename T2,
676  bool is_const_T1, bool is_const_T2,
677  typename ExprSpec>
678  class IfThenElseOp {};
679 
680  template <typename CondT, typename T1, typename T2>
681  class IfThenElseOp< CondT, T1, T2, false, false, ExprSpecDefault > :
682  public Expr< IfThenElseOp< CondT, T1, T2, false, false, ExprSpecDefault > > {
683 
684  public:
685 
686  typedef typename std::remove_cv<T1>::type ExprT1;
687  typedef typename std::remove_cv<T2>::type ExprT2;
688  typedef typename ExprT1::value_type value_type_1;
689  typedef typename ExprT2::value_type value_type_2;
690  typedef typename Sacado::Promote<value_type_1,
691  value_type_2>::type value_type;
692 
693  typedef typename ExprT1::scalar_type scalar_type_1;
694  typedef typename ExprT2::scalar_type scalar_type_2;
695  typedef typename Sacado::Promote<scalar_type_1,
696  scalar_type_2>::type scalar_type;
697 
698  typedef ExprSpecDefault expr_spec_type;
699 
701  IfThenElseOp(const CondT& cond_, const T1& expr1_, const T2& expr2_) :
702  cond(cond_), expr1(expr1_), expr2(expr2_) {}
703 
705  int size() const {
706  int sz1 = expr1.size(), sz2 = expr2.size();
707  return sz1 > sz2 ? sz1 : sz2;
708  }
709 
711  bool hasFastAccess() const {
712  return expr1.hasFastAccess() && expr2.hasFastAccess();
713  }
714 
716  value_type val() const {
717  return if_then_else( cond, expr1.val(), expr2.val() );
718  }
719 
721  value_type dx(int i) const {
722  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
723  }
724 
726  value_type fastAccessDx(int i) const {
727  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
728  }
729 
730  protected:
731 
732  const CondT& cond;
733  const T1& expr1;
734  const T2& expr2;
735 
736  };
737 
738  template <typename CondT, typename T1, typename T2>
739  class IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > :
740  public Expr< IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > > {
741 
742  public:
743 
744  typedef typename std::remove_cv<T1>::type ExprT1;
745  typedef T2 ConstT;
746  typedef typename ExprT1::value_type value_type;
747  typedef typename ExprT1::scalar_type scalar_type;
748 
749  typedef ExprSpecDefault expr_spec_type;
750 
752  IfThenElseOp(const CondT& cond_, const T1& expr1_, const ConstT& c_) :
753  cond(cond_), expr1(expr1_), c(c_) {}
754 
756  int size() const {
757  return expr1.size();
758  }
759 
761  bool hasFastAccess() const {
762  return expr1.hasFastAccess();
763  }
764 
766  value_type val() const {
767  return if_then_else( cond, expr1.val(), c );
768  }
769 
771  value_type dx(int i) const {
772  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
773  }
774 
776  value_type fastAccessDx(int i) const {
777  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
778  }
779 
780  protected:
781 
782  const CondT& cond;
783  const T1& expr1;
784  const ConstT& c;
785  };
786 
787  template <typename CondT, typename T1, typename T2>
788  class IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > :
789  public Expr< IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > > {
790 
791  public:
792 
793  typedef typename std::remove_cv<T2>::type ExprT2;
794  typedef T1 ConstT;
795  typedef typename ExprT2::value_type value_type;
796  typedef typename ExprT2::scalar_type scalar_type;
797 
798  typedef ExprSpecDefault expr_spec_type;
799 
801  IfThenElseOp(const CondT& cond_, const ConstT& c_, const T2& expr2_) :
802  cond(cond_), c(c_), expr2(expr2_) {}
803 
805  int size() const {
806  return expr2.size();
807  }
808 
810  bool hasFastAccess() const {
811  return expr2.hasFastAccess();
812  }
813 
815  value_type val() const {
816  return if_then_else( cond, c, expr2.val() );
817  }
818 
820  value_type dx(int i) const {
821  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
822  }
823 
825  value_type fastAccessDx(int i) const {
826  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
827  }
828 
829  protected:
830 
831  const CondT& cond;
832  const ConstT& c;
833  const T2& expr2;
834  };
835 
836  template <typename CondT, typename T1, typename T2>
838  typename mpl::enable_if_c< IsFadExpr<T1>::value && IsFadExpr<T2>::value &&
839  ExprLevel<T1>::value == ExprLevel<T2>::value,
840  IfThenElseOp< CondT,
841  typename Expr<T1>::derived_type,
842  typename Expr<T2>::derived_type,
843  false, false,
844  typename T1::expr_spec_type >
845  >::type
846  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
847  {
848  typedef IfThenElseOp< CondT, typename Expr<T1>::derived_type,
849  typename Expr<T2>::derived_type,
850  false, false, typename T1::expr_spec_type > expr_t;
851 
852  return expr_t(cond, expr1.derived(), expr2.derived());
853  }
854 
855  template <typename CondT, typename T>
857  IfThenElseOp< CondT, typename T::value_type, typename Expr<T>::derived_type,
858  true, false, typename T::expr_spec_type >
859  if_then_else (const CondT& cond, const typename T::value_type& c,
860  const Expr<T>& expr)
861  {
862  typedef typename T::value_type ConstT;
863  typedef IfThenElseOp< CondT, ConstT, typename Expr<T>::derived_type,
864  true, false, typename T::expr_spec_type > expr_t;
865 
866  return expr_t(cond, c, expr.derived());
867  }
868 
869  template <typename CondT, typename T>
871  IfThenElseOp< CondT, typename Expr<T>::derived_type, typename T::value_type,
872  false, true, typename T::expr_spec_type >
873  if_then_else (const CondT& cond, const Expr<T>& expr,
874  const typename T::value_type& c)
875  {
876  typedef typename T::value_type ConstT;
877  typedef IfThenElseOp< CondT, typename Expr<T>::derived_type, ConstT,
878  false, true, typename T::expr_spec_type > expr_t;
879 
880  return expr_t(cond, expr.derived(), c);
881  }
882 
883  template <typename CondT, typename T>
885  typename mpl::disable_if<
886  mpl::is_same< typename T::value_type,
887  typename T::scalar_type >,
888  IfThenElseOp< CondT, typename T::scalar_type,
889  typename Expr<T>::derived_type,
890  true, false, typename T::expr_spec_type >
891  >::type
892  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
893  const Expr<T>& expr)
894  {
895  typedef typename T::scalar_type ConstT;
896  typedef IfThenElseOp< CondT, ConstT, typename Expr<T>::derived_type,
897  true, false, typename T::expr_spec_type > expr_t;
898 
899  return expr_t(cond, c, expr.derived());
900  }
901 
902  template <typename CondT, typename T>
904  typename mpl::disable_if<
905  mpl::is_same< typename T::value_type,
906  typename T::scalar_type >,
907  IfThenElseOp< CondT, typename Expr<T>::derived_type,
908  typename T::scalar_type,
909  false, true, typename T::expr_spec_type >
910  >::type
911  if_then_else (const CondT& cond, const Expr<T>& expr,
912  const typename Expr<T>::scalar_type& c)
913  {
914  typedef typename T::scalar_type ConstT;
915  typedef IfThenElseOp< CondT, typename Expr<T>::derived_type, ConstT,
916  false, true, typename T::expr_spec_type > expr_t;
917 
918  return expr_t(cond, expr.derived(), c);
919  }
920 
921  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
922  typename E>
923  struct ExprLevel< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
924  static constexpr unsigned value_1 = ExprLevel<T1>::value;
925  static constexpr unsigned value_2 = ExprLevel<T2>::value;
926  static constexpr unsigned value =
927  value_1 >= value_2 ? value_1 : value_2;
928  };
929 
930  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
931  typename E>
932  struct IsFadExpr< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
933  static constexpr unsigned value = true;
934  };
935 
936  }
937  }
938 
939  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
940  typename E>
941  struct IsExpr< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
942  static constexpr bool value = true;
943  };
944 
945  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
946  typename E>
947  struct BaseExprType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
948  typedef typename BaseExprType<T1>::type base_expr_1;
949  typedef typename BaseExprType<T2>::type base_expr_2;
950  typedef typename Sacado::Promote<base_expr_1,
951  base_expr_2>::type type;
952  };
953 }
954 
955 #undef FAD_BINARYOP_MACRO
956 
957 //-------------------------- Relational Operators -----------------------
958 
959 namespace Sacado {
960  namespace Fad {
961  namespace Exp {
962  template <typename T1, typename T2 = T1>
964  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
965  };
966  }
967  }
968 }
969 
970 #define FAD_RELOP_MACRO(OP) \
971 namespace Sacado { \
972  namespace Fad { \
973  namespace Exp { \
974  template <typename T1, typename T2> \
975  KOKKOS_INLINE_FUNCTION \
976  typename mpl::enable_if_c< \
977  IsFadExpr<T1>::value && IsFadExpr<T2>::value && \
978  ExprLevel<T1>::value == ExprLevel<T2>::value, \
979  typename ConditionalReturnType<typename T1::value_type, \
980  typename T2::value_type>::type \
981  >::type \
982  operator OP (const T1& expr1, const T2& expr2) \
983  { \
984  return expr1.derived().val() OP expr2.derived().val(); \
985  } \
986  \
987  template <typename T2> \
988  KOKKOS_INLINE_FUNCTION \
989  typename ConditionalReturnType<typename T2::value_type>::type \
990  operator OP (const typename T2::value_type& a, \
991  const Expr<T2>& expr2) \
992  { \
993  return a OP expr2.derived().val(); \
994  } \
995  \
996  template <typename T1> \
997  KOKKOS_INLINE_FUNCTION \
998  typename ConditionalReturnType<typename T1::value_type>::type \
999  operator OP (const Expr<T1>& expr1, \
1000  const typename T1::value_type& b) \
1001  { \
1002  return expr1.derived().val() OP b; \
1003  } \
1004  } \
1005  } \
1006 } \
1007 
1008 FAD_RELOP_MACRO(==)
1009 FAD_RELOP_MACRO(!=)
1010 FAD_RELOP_MACRO(<)
1011 FAD_RELOP_MACRO(>)
1012 FAD_RELOP_MACRO(<=)
1013 FAD_RELOP_MACRO(>=)
1014 FAD_RELOP_MACRO(<<=)
1015 FAD_RELOP_MACRO(>>=)
1016 FAD_RELOP_MACRO(&)
1017 FAD_RELOP_MACRO(|)
1018 
1019 #undef FAD_RELOP_MACRO
1020 
1021 namespace Sacado {
1022 
1023  namespace Fad {
1024  namespace Exp {
1025 
1026  template <typename ExprT>
1028  bool operator ! (const Expr<ExprT>& expr)
1029  {
1030  return ! expr.derived().val();
1031  }
1032 
1033  } // namespace Exp
1034  } // namespace Fad
1035 
1036 } // namespace Sacado
1037 
1038 //-------------------------- Boolean Operators -----------------------
1039 namespace Sacado {
1040 
1041  namespace Fad {
1042  namespace Exp {
1043 
1044  template <typename T>
1046  bool toBool(const Expr<T>& xx) {
1047  const typename Expr<T>::derived_type& x = xx.derived();
1048  bool is_zero = (x.val() == 0.0);
1049  for (int i=0; i<x.size(); i++)
1050  is_zero = is_zero && (x.dx(i) == 0.0);
1051  return !is_zero;
1052  }
1053 
1054  } // namespace Exp
1055  } // namespace Fad
1056 
1057 } // namespace Sacado
1058 
1059 #define FAD_BOOL_MACRO(OP) \
1060 namespace Sacado { \
1061  namespace Fad { \
1062  namespace Exp { \
1063  template <typename T1, typename T2> \
1064  KOKKOS_INLINE_FUNCTION \
1065  bool \
1066  operator OP (const Expr<T1>& expr1, \
1067  const Expr<T2>& expr2) \
1068  { \
1069  return toBool(expr1) OP toBool(expr2); \
1070  } \
1071  \
1072  template <typename T2> \
1073  KOKKOS_INLINE_FUNCTION \
1074  bool \
1075  operator OP (const typename Expr<T2>::value_type& a, \
1076  const Expr<T2>& expr2) \
1077  { \
1078  return a OP toBool(expr2); \
1079  } \
1080  \
1081  template <typename T1> \
1082  KOKKOS_INLINE_FUNCTION \
1083  bool \
1084  operator OP (const Expr<T1>& expr1, \
1085  const typename Expr<T1>::value_type& b) \
1086  { \
1087  return toBool(expr1) OP b; \
1088  } \
1089  } \
1090  } \
1091 }
1092 
1093 FAD_BOOL_MACRO(&&)
1094 FAD_BOOL_MACRO(||)
1095 
1096 #undef FAD_BOOL_MACRO
1097 
1098 //-------------------------- I/O Operators -----------------------
1099 
1100 namespace Sacado {
1101 
1102  namespace Fad {
1103  namespace Exp {
1104 
1105  template <typename T>
1106  std::ostream& operator << (std::ostream& os, const Expr<T>& xx) {
1107  const typename Expr<T>::derived_type& x = xx.derived();
1108  os << x.val() << " [";
1109 
1110  for (int i=0; i< x.size(); i++) {
1111  os << " " << x.dx(i);
1112  }
1113 
1114  os << " ]";
1115  return os;
1116  }
1117 
1118  } // namespace Exp
1119  } // namespace Fad
1120 
1121 } // namespace Sacado
1122 
1123 #if defined(HAVE_SACADO_KOKKOSCORE)
1124 
1125 //-------------------------- Atomic Operators -----------------------
1126 
1127 namespace Sacado {
1128 
1129  namespace Fad {
1130  namespace Exp {
1131 
1132  // Overload of Kokkos::atomic_add for Fad types.
1133  template <typename S>
1135  void atomic_add(GeneralFad<S>* dst, const GeneralFad<S>& x) {
1136  using Kokkos::atomic_add;
1137 
1138  const int xsz = x.size();
1139  const int sz = dst->size();
1140 
1141  // We currently cannot handle resizing since that would need to be
1142  // done atomically.
1143  if (xsz > sz)
1144  Kokkos::abort(
1145  "Sacado error: Fad resize within atomic_add() not supported!");
1146 
1147  if (xsz != sz && sz > 0 && xsz > 0)
1148  Kokkos::abort(
1149  "Sacado error: Fad assignment of incompatiable sizes!");
1150 
1151 
1152  if (sz > 0 && xsz > 0) {
1153  SACADO_FAD_DERIV_LOOP(i,sz)
1154  atomic_add(&(dst->fastAccessDx(i)), x.fastAccessDx(i));
1155  }
1157  atomic_add(&(dst->val()), x.val());
1158  }
1159 
1160  } // namespace Exp
1161  } // namespace Fad
1162 
1163 } // namespace Sacado
1164 
1165 #endif // HAVE_SACADO_KOKKOSCORE
1166 
1167 #endif // SACADO_FAD_OPS_HPP
cbrt(expr.val())
#define FAD_BOOL_MACRO(OP)
expr expr AdditionOp
expr expr SinOp
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())
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
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
expr expr TanhOp
pow(expr1.val(), expr2.val())
expr expr SqrtOp
expr expr ASinhOp
expr true
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
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)
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)
T derived_type
Typename of derived object, returned by derived()
sqrt(expr.val())
sinh(expr.val())
tan(expr.val())
KOKKOS_INLINE_FUNCTION const derived_type & derived() const
Return derived object.
#define T1(r, f)
Definition: Sacado_rad.hpp:603
#define FAD_RELOP_MACRO(OP)
atan2(expr1.val(), expr2.val())
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
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
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 PowerOp
fabs(expr.val())
expr expr AbsOp
expr expr TanOp
log10(expr.val())
Base template specification for Promote.
cos(expr.val())