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  explicit 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;,
236  tanh(expr.val()),
237  expr.dx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())),
238  expr.fastAccessDx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())))
241  using std::acosh; using std::sqrt;,
242  acosh(expr.val()),
243  expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
244  (expr.val()+value_type(1))),
245  expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
246  (expr.val()+value_type(1))))
249  using std::asinh; using std::sqrt;,
250  asinh(expr.val()),
251  expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
252  expr.fastAccessDx(i)/ sqrt(value_type(1)+
253  expr.val()*expr.val()))
256  using std::atanh;,
257  atanh(expr.val()),
258  expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
259  expr.fastAccessDx(i)/(value_type(1)-
260  expr.val()*expr.val()))
263  using std::abs;,
264  abs(expr.val()),
265  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
266  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
269  using std::fabs;,
270  fabs(expr.val()),
271  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
272  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
275  using std::cbrt;,
276  cbrt(expr.val()),
277  expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
278  expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
279 
280 // Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
281 // "simd" value types that use if_then_else(). The only reason for not using
282 // if_then_else() always is to avoid evaluating the derivative if the value is
283 // zero to avoid throwing FPEs.
284 namespace Sacado {
285  namespace Fad {
286  namespace Exp {
287 
288  template <typename T, typename ExprSpec, bool is_simd>
289  class SafeSqrtOp {};
290 
291  //
292  // Implementation for simd type using if_then_else()
293  //
294  template <typename T>
295  class SafeSqrtOp< T,ExprSpecDefault,true > :
296  public Expr< SafeSqrtOp<T,ExprSpecDefault> > {
297  public:
298 
299  typedef typename std::remove_cv<T>::type ExprT;
300  typedef typename ExprT::value_type value_type;
301  typedef typename ExprT::scalar_type scalar_type;
302 
303  typedef ExprSpecDefault expr_spec_type;
304 
306  explicit SafeSqrtOp(const T& expr_) : expr(expr_) {}
307 
309  int size() const { return expr.size(); }
310 
312  bool hasFastAccess() const {
313  return expr.hasFastAccess();
314  }
315 
317  value_type val() const {
318  using std::sqrt;
319  return sqrt(expr.val());
320  }
321 
323  value_type dx(int i) const {
324  using std::sqrt;
325  return if_then_else(
326  expr.val() == value_type(0.0), value_type(0.0),
327  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
328  }
329 
331  value_type fastAccessDx(int i) const {
332  using std::sqrt;
333  return if_then_else(
334  expr.val() == value_type(0.0), value_type(0.0),
335  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
336  }
337 
338  protected:
339 
340  const T& expr;
341  };
342 
343  //
344  // Specialization for scalar types using ternary operator
345  //
346  template <typename T>
347  class SafeSqrtOp< T,ExprSpecDefault,false > :
348  public Expr< SafeSqrtOp<T,ExprSpecDefault> > {
349  public:
350 
351  typedef typename std::remove_cv<T>::type ExprT;
352  typedef typename ExprT::value_type value_type;
353  typedef typename ExprT::scalar_type scalar_type;
354 
355  typedef ExprSpecDefault expr_spec_type;
356 
358  explicit SafeSqrtOp(const T& expr_) : expr(expr_) {}
359 
361  int size() const { return expr.size(); }
362 
364  bool hasFastAccess() const {
365  return expr.hasFastAccess();
366  }
367 
369  value_type val() const {
370  using std::sqrt;
371  return sqrt(expr.val());
372  }
373 
375  value_type dx(int i) const {
376  using std::sqrt;
377  return expr.val() == value_type(0.0) ? value_type(0.0) :
378  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
379  }
380 
382  value_type fastAccessDx(int i) const {
383  using std::sqrt;
384  return expr.val() == value_type(0.0) ? value_type(0.0) :
385  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
386  }
387 
388  protected:
389 
390  const T& expr;
391  };
392 
393  template <typename T>
395  SafeSqrtOp< typename Expr<T>::derived_type,
396  typename T::expr_spec_type >
397  safe_sqrt (const Expr<T>& expr)
398  {
399  typedef SafeSqrtOp< typename Expr<T>::derived_type,
400  typename T::expr_spec_type > expr_t;
401 
402  return expr_t(expr.derived());
403  }
404 
405  template <typename T, typename E>
406  struct ExprLevel< SafeSqrtOp< T,E > > {
407  static const unsigned value = ExprLevel<T>::value;
408  };
409 
410  template <typename T, typename E>
411  struct IsFadExpr< SafeSqrtOp< T,E > > {
412  static const unsigned value = true;
413  };
414 
415  }
416  }
417 
418  template <typename T, typename E>
419  struct IsExpr< Fad::Exp::SafeSqrtOp< T,E > > {
420  static const bool value = true;
421  };
422 
423  template <typename T, typename E>
424  struct BaseExprType< Fad::Exp::SafeSqrtOp< T,E > > {
425  typedef typename BaseExprType<T>::type type;
426  };
427 
428  template <typename T, typename E>
429  struct IsSimdType< Fad::Exp::SafeSqrtOp< T,E > > {
430  static const bool value =
431  IsSimdType< typename Fad::Exp::SafeSqrtOp< T,E >::scalar_type >::value;
432  };
433 
434 }
435 
436 #undef FAD_UNARYOP_MACRO
437 
438 #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) \
439 namespace Sacado { \
440  namespace Fad { \
441  namespace Exp { \
442  \
443  template <typename T1, typename T2, \
444  bool is_const_T1, bool is_const_T2, \
445  typename ExprSpec > \
446  class OP {}; \
447  \
448  template <typename T1, typename T2> \
449  class OP< T1, T2, false, false, ExprSpecDefault > : \
450  public Expr< OP< T1, T2, false, false, ExprSpecDefault > > { \
451  public: \
452  \
453  typedef typename std::remove_cv<T1>::type ExprT1; \
454  typedef typename std::remove_cv<T2>::type ExprT2; \
455  typedef typename ExprT1::value_type value_type_1; \
456  typedef typename ExprT2::value_type value_type_2; \
457  typedef typename Sacado::Promote<value_type_1, \
458  value_type_2>::type value_type; \
459  \
460  typedef typename ExprT1::scalar_type scalar_type_1; \
461  typedef typename ExprT2::scalar_type scalar_type_2; \
462  typedef typename Sacado::Promote<scalar_type_1, \
463  scalar_type_2>::type scalar_type; \
464  \
465  typedef ExprSpecDefault expr_spec_type; \
466  \
467  KOKKOS_INLINE_FUNCTION \
468  OP(const T1& expr1_, const T2& expr2_) : \
469  expr1(expr1_), expr2(expr2_) {} \
470  \
471  KOKKOS_INLINE_FUNCTION \
472  int size() const { \
473  const int sz1 = expr1.size(), sz2 = expr2.size(); \
474  return sz1 > sz2 ? sz1 : sz2; \
475  } \
476  \
477  KOKKOS_INLINE_FUNCTION \
478  bool hasFastAccess() const { \
479  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
480  } \
481  \
482  KOKKOS_INLINE_FUNCTION \
483  value_type val() const { \
484  USING \
485  return VALUE; \
486  } \
487  \
488  KOKKOS_INLINE_FUNCTION \
489  value_type dx(int i) const { \
490  USING \
491  const int sz1 = expr1.size(), sz2 = expr2.size(); \
492  if (sz1 > 0 && sz2 > 0) \
493  return DX; \
494  else if (sz1 > 0) \
495  return CDX2; \
496  else \
497  return CDX1; \
498  } \
499  \
500  KOKKOS_INLINE_FUNCTION \
501  value_type fastAccessDx(int i) const { \
502  USING \
503  return FASTACCESSDX; \
504  } \
505  \
506  protected: \
507  \
508  const T1& expr1; \
509  const T2& expr2; \
510  \
511  }; \
512  \
513  template <typename T1, typename T2> \
514  class OP< T1, T2, false, true, ExprSpecDefault > \
515  : public Expr< OP< T1, T2, false, true, ExprSpecDefault > > { \
516  public: \
517  \
518  typedef typename std::remove_cv<T1>::type ExprT1; \
519  typedef T2 ConstT; \
520  typedef typename ExprT1::value_type value_type; \
521  typedef typename ExprT1::scalar_type scalar_type; \
522  \
523  typedef ExprSpecDefault expr_spec_type; \
524  \
525  KOKKOS_INLINE_FUNCTION \
526  OP(const T1& expr1_, const ConstT& c_) : \
527  expr1(expr1_), c(c_) {} \
528  \
529  KOKKOS_INLINE_FUNCTION \
530  int size() const { \
531  return expr1.size(); \
532  } \
533  \
534  KOKKOS_INLINE_FUNCTION \
535  bool hasFastAccess() const { \
536  return expr1.hasFastAccess(); \
537  } \
538  \
539  KOKKOS_INLINE_FUNCTION \
540  value_type val() const { \
541  USING \
542  return VAL_CONST_DX_2; \
543  } \
544  \
545  KOKKOS_INLINE_FUNCTION \
546  value_type dx(int i) const { \
547  USING \
548  return CONST_DX_2; \
549  } \
550  \
551  KOKKOS_INLINE_FUNCTION \
552  value_type fastAccessDx(int i) const { \
553  USING \
554  return CONST_FASTACCESSDX_2; \
555  } \
556  \
557  protected: \
558  \
559  const T1& expr1; \
560  const ConstT& c; \
561  }; \
562  \
563  template <typename T1, typename T2> \
564  class OP< T1, T2, true, false, ExprSpecDefault > \
565  : public Expr< OP< T1, T2, true, false, ExprSpecDefault > > { \
566  public: \
567  \
568  typedef typename std::remove_cv<T2>::type ExprT2; \
569  typedef T1 ConstT; \
570  typedef typename ExprT2::value_type value_type; \
571  typedef typename ExprT2::scalar_type scalar_type; \
572  \
573  typedef ExprSpecDefault expr_spec_type; \
574  \
575  KOKKOS_INLINE_FUNCTION \
576  OP(const ConstT& c_, const T2& expr2_) : \
577  c(c_), expr2(expr2_) {} \
578  \
579  KOKKOS_INLINE_FUNCTION \
580  int size() const { \
581  return expr2.size(); \
582  } \
583  \
584  KOKKOS_INLINE_FUNCTION \
585  bool hasFastAccess() const { \
586  return expr2.hasFastAccess(); \
587  } \
588  \
589  KOKKOS_INLINE_FUNCTION \
590  value_type val() const { \
591  USING \
592  return VAL_CONST_DX_1; \
593  } \
594  \
595  KOKKOS_INLINE_FUNCTION \
596  value_type dx(int i) const { \
597  USING \
598  return CONST_DX_1; \
599  } \
600  \
601  KOKKOS_INLINE_FUNCTION \
602  value_type fastAccessDx(int i) const { \
603  USING \
604  return CONST_FASTACCESSDX_1; \
605  } \
606  \
607  protected: \
608  \
609  const ConstT& c; \
610  const T2& expr2; \
611  }; \
612  \
613  template <typename T1, typename T2> \
614  KOKKOS_INLINE_FUNCTION \
615  SACADO_FAD_EXP_OP_ENABLE_EXPR_EXPR(OP) \
616  OPNAME (const T1& expr1, const T2& expr2) \
617  { \
618  typedef OP< typename Expr<T1>::derived_type, \
619  typename Expr<T2>::derived_type, \
620  false, false, typename T1::expr_spec_type > expr_t; \
621  \
622  return expr_t(expr1.derived(), expr2.derived()); \
623  } \
624  \
625  template <typename T> \
626  KOKKOS_INLINE_FUNCTION \
627  OP< typename T::value_type, typename Expr<T>::derived_type, \
628  true, false, typename T::expr_spec_type > \
629  OPNAME (const typename T::value_type& c, \
630  const Expr<T>& expr) \
631  { \
632  typedef typename T::value_type ConstT; \
633  typedef OP< ConstT, typename Expr<T>::derived_type, \
634  true, false, typename T::expr_spec_type > expr_t; \
635  \
636  return expr_t(c, expr.derived()); \
637  } \
638  \
639  template <typename T> \
640  KOKKOS_INLINE_FUNCTION \
641  OP< typename Expr<T>::derived_type, typename T::value_type, \
642  false, true, typename T::expr_spec_type > \
643  OPNAME (const Expr<T>& expr, \
644  const typename T::value_type& c) \
645  { \
646  typedef typename T::value_type ConstT; \
647  typedef OP< typename Expr<T>::derived_type, ConstT, \
648  false, true, typename T::expr_spec_type > expr_t; \
649  \
650  return expr_t(expr.derived(), c); \
651  } \
652  \
653  template <typename T> \
654  KOKKOS_INLINE_FUNCTION \
655  SACADO_FAD_EXP_OP_ENABLE_SCALAR_EXPR(OP) \
656  OPNAME (const typename T::scalar_type& c, \
657  const Expr<T>& expr) \
658  { \
659  typedef typename T::scalar_type ConstT; \
660  typedef OP< ConstT, typename Expr<T>::derived_type, \
661  true, false, typename T::expr_spec_type > expr_t; \
662  \
663  return expr_t(c, expr.derived()); \
664  } \
665  \
666  template <typename T> \
667  KOKKOS_INLINE_FUNCTION \
668  SACADO_FAD_EXP_OP_ENABLE_EXPR_SCALAR(OP) \
669  OPNAME (const Expr<T>& expr, \
670  const typename T::scalar_type& c) \
671  { \
672  typedef typename T::scalar_type ConstT; \
673  typedef OP< typename Expr<T>::derived_type, ConstT, \
674  false, true, typename T::expr_spec_type > expr_t; \
675  \
676  return expr_t(expr.derived(), c); \
677  } \
678  \
679  template <typename T1, typename T2, bool c1, bool c2, typename E> \
680  struct ExprLevel< OP< T1, T2, c1, c2, E > > { \
681  static constexpr unsigned value_1 = ExprLevel<T1>::value; \
682  static constexpr unsigned value_2 = ExprLevel<T2>::value; \
683  static constexpr unsigned value = \
684  value_1 >= value_2 ? value_1 : value_2; \
685  }; \
686  \
687  template <typename T1, typename T2, bool c1, bool c2, typename E> \
688  struct IsFadExpr< OP< T1, T2, c1, c2, E > > { \
689  static constexpr unsigned value = true; \
690  }; \
691  \
692  } \
693  } \
694  \
695  template <typename T1, typename T2, bool c1, bool c2, typename E> \
696  struct IsExpr< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
697  static constexpr bool value = true; \
698  }; \
699  \
700  template <typename T1, typename T2, bool c1, bool c2, typename E> \
701  struct BaseExprType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
702  typedef typename BaseExprType<T1>::type base_expr_1; \
703  typedef typename BaseExprType<T2>::type base_expr_2; \
704  typedef typename Sacado::Promote<base_expr_1, \
705  base_expr_2>::type type; \
706  }; \
707  \
708  template <typename T1, typename T2, bool c1, bool c2, typename E> \
709  struct IsSimdType< Fad::Exp::OP< T1, T2, c1, c2, E > > { \
710  static const bool value = \
711  IsSimdType< typename Fad::Exp::OP< T1, T2, c1, c2, E >::value_type >::value; \
712  }; \
713  \
714 }
715 
716 
717 FAD_BINARYOP_MACRO(operator+,
718  AdditionOp,
719  ;,
720  expr1.val() + expr2.val(),
721  expr1.dx(i) + expr2.dx(i),
722  expr2.dx(i),
723  expr1.dx(i),
724  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
725  c + expr2.val(),
726  expr1.val() + c,
727  expr2.dx(i),
728  expr1.dx(i),
729  expr2.fastAccessDx(i),
730  expr1.fastAccessDx(i))
731 FAD_BINARYOP_MACRO(operator-,
733  ;,
734  expr1.val() - expr2.val(),
735  expr1.dx(i) - expr2.dx(i),
736  -expr2.dx(i),
737  expr1.dx(i),
738  expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
739  c - expr2.val(),
740  expr1.val() - c,
741  -expr2.dx(i),
742  expr1.dx(i),
743  -expr2.fastAccessDx(i),
744  expr1.fastAccessDx(i))
745 FAD_BINARYOP_MACRO(operator*,
747  ;,
748  expr1.val() * expr2.val(),
749  expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
750  expr1.val()*expr2.dx(i),
751  expr1.dx(i)*expr2.val(),
752  expr1.val()*expr2.fastAccessDx(i) +
753  expr1.fastAccessDx(i)*expr2.val(),
754  c * expr2.val(),
755  expr1.val() * c,
756  c*expr2.dx(i),
757  expr1.dx(i)*c,
758  c*expr2.fastAccessDx(i),
759  expr1.fastAccessDx(i)*c)
760 FAD_BINARYOP_MACRO(operator/,
762  ;,
763  expr1.val() / expr2.val(),
764  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
765  (expr2.val()*expr2.val()),
766  -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
767  expr1.dx(i)/expr2.val(),
768  (expr1.fastAccessDx(i)*expr2.val() -
769  expr2.fastAccessDx(i)*expr1.val()) /
770  (expr2.val()*expr2.val()),
771  c / expr2.val(),
772  expr1.val() / c,
773  -expr2.dx(i)*c / (expr2.val()*expr2.val()),
774  expr1.dx(i)/c,
775  -expr2.fastAccessDx(i)*c / (expr2.val()*expr2.val()),
776  expr1.fastAccessDx(i)/c)
779  using std::atan2;,
780  atan2(expr1.val(), expr2.val()),
781  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
782  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
783  -expr1.val()*expr2.dx(i)/
784  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
785  expr2.val()*expr1.dx(i)/
786  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
787  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
788  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
789  atan2(c, expr2.val()),
790  atan2(expr1.val(), c),
791  (-c*expr2.dx(i)) / (c*c + expr2.val()*expr2.val()),
792  (c*expr1.dx(i))/ (expr1.val()*expr1.val() + c*c),
793  (-c*expr2.fastAccessDx(i))/ (c*c + expr2.val()*expr2.val()),
794  (c*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c*c))
795 // FAD_BINARYOP_MACRO(pow,
796 // PowerOp,
797 // using std::pow; using std::log;,
798 // pow(expr1.val(), expr2.val()),
799 // 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())) ),
800 // 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())) ),
801 // 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())) ),
802 // 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())) ),
803 // pow(c, expr2.val()),
804 // pow(expr1.val(), c),
805 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) ),
806 // 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)) ),
807 // if_then_else( c == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) ),
808 // if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c))) )
811  ;,
812  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
813  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
814  if_then_else( expr1.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
815  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), value_type(0.0) ),
816  if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
817  if_then_else( c >= expr2.val(), value_type(c), expr2.val() ),
818  if_then_else( expr1.val() >= c, expr1.val(), value_type(c) ),
819  if_then_else( c >= expr2.val(), value_type(0.0), expr2.dx(i) ),
820  if_then_else( expr1.val() >= c, expr1.dx(i), value_type(0.0) ),
821  if_then_else( c >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
822  if_then_else( expr1.val() >= c, expr1.fastAccessDx(i), value_type(0.0) ) )
825  ;,
826  if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
827  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
828  if_then_else( expr1.val() <= expr2.val(), value_type(0.0), expr2.dx(i) ),
829  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), value_type(0.0) ),
830  if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
831  if_then_else( c <= expr2.val(), value_type(c), expr2.val() ),
832  if_then_else( expr1.val() <= c, expr1.val(), value_type(c) ),
833  if_then_else( c <= expr2.val(), value_type(0), expr2.dx(i) ),
834  if_then_else( expr1.val() <= c, expr1.dx(i), value_type(0) ),
835  if_then_else( c <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
836  if_then_else( expr1.val() <= c, expr1.fastAccessDx(i), value_type(0) ) )
837 
838 
839 // Special handling for std::pow() to provide specializations of PowerOp for
840 // "simd" value types that use if_then_else(). The only reason for not using
841 // if_then_else() always is to avoid evaluating the derivative if the value is
842 // zero to avoid throwing FPEs.
843 namespace Sacado {
844  namespace Fad {
845  namespace Exp {
846 
847  template <typename T1, typename T2,
848  bool is_const_T1, bool is_const_T2,
849  typename ExprSpec, bool is_simd >
850  class PowerOp {};
851 
852  //
853  // Implementation for simd type using if_then_else()
854  //
855  template <typename T1, typename T2>
856  class PowerOp< T1, T2, false, false, ExprSpecDefault, true > :
857  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault, true > > {
858  public:
859 
860  typedef typename std::remove_cv<T1>::type ExprT1;
861  typedef typename std::remove_cv<T2>::type ExprT2;
862  typedef typename ExprT1::value_type value_type_1;
863  typedef typename ExprT2::value_type value_type_2;
864  typedef typename Sacado::Promote<value_type_1,
865  value_type_2>::type value_type;
866 
867  typedef typename ExprT1::scalar_type scalar_type_1;
868  typedef typename ExprT2::scalar_type scalar_type_2;
869  typedef typename Sacado::Promote<scalar_type_1,
870  scalar_type_2>::type scalar_type;
871 
872  typedef ExprSpecDefault expr_spec_type;
873 
875  PowerOp(const T1& expr1_, const T2& expr2_) :
876  expr1(expr1_), expr2(expr2_) {}
877 
879  int size() const {
880  const int sz1 = expr1.size(), sz2 = expr2.size();
881  return sz1 > sz2 ? sz1 : sz2;
882  }
883 
885  bool hasFastAccess() const {
886  return expr1.hasFastAccess() && expr2.hasFastAccess();
887  }
888 
890  value_type val() const {
891  using std::pow;
892  return pow(expr1.val(), expr2.val());
893  }
894 
896  value_type dx(int i) const {
897  using std::pow; using std::log;
898  const int sz1 = expr1.size(), sz2 = expr2.size();
899  if (sz1 > 0 && sz2 > 0)
900  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())) );
901  else if (sz1 > 0)
902  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
903  // It seems less accurate and caused convergence problems in some codes
904  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())) );
905  else
906  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())) );
907  }
908 
910  value_type fastAccessDx(int i) const {
911  using std::pow; using std::log;
912  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())) );
913  }
914 
915  protected:
916 
917  const T1& expr1;
918  const T2& expr2;
919 
920  };
921 
922  template <typename T1, typename T2>
923  class PowerOp< T1, T2, false, true, ExprSpecDefault, true >
924  : public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault, true > > {
925  public:
926 
927  typedef typename std::remove_cv<T1>::type ExprT1;
928  typedef T2 ConstT;
929  typedef typename ExprT1::value_type value_type;
930  typedef typename ExprT1::scalar_type scalar_type;
931 
932  typedef ExprSpecDefault expr_spec_type;
933 
935  PowerOp(const T1& expr1_, const ConstT& c_) :
936  expr1(expr1_), c(c_) {}
937 
939  int size() const {
940  return expr1.size();
941  }
942 
944  bool hasFastAccess() const {
945  return expr1.hasFastAccess();
946  }
947 
949  value_type val() const {
950  using std::pow;
951  return pow(expr1.val(), c);
952  }
953 
955  value_type dx(int i) const {
956  using std::pow;
957  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
958  // It seems less accurate and caused convergence problems in some codes
959  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)) );
960  }
961 
963  value_type fastAccessDx(int i) const {
964  using std::pow;
965  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
966  // It seems less accurate and caused convergence problems in some codes
967  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)) );
968  }
969 
970  protected:
971 
972  const T1& expr1;
973  const ConstT& c;
974  };
975 
976  template <typename T1, typename T2>
977  class PowerOp< T1, T2, true, false, ExprSpecDefault, true >
978  : public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault, true > > {
979  public:
980 
981  typedef typename std::remove_cv<T2>::type ExprT2;
982  typedef T1 ConstT;
983  typedef typename ExprT2::value_type value_type;
984  typedef typename ExprT2::scalar_type scalar_type;
985 
986  typedef ExprSpecDefault expr_spec_type;
987 
989  PowerOp(const ConstT& c_, const T2& expr2_) :
990  c(c_), expr2(expr2_) {}
991 
993  int size() const {
994  return expr2.size();
995  }
996 
998  bool hasFastAccess() const {
999  return expr2.hasFastAccess();
1000  }
1001 
1003  value_type val() const {
1004  using std::pow;
1005  return pow(c, expr2.val());
1006  }
1007 
1009  value_type dx(int i) const {
1010  using std::pow; using std::log;
1011  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c)*pow(c,expr2.val())) );
1012  }
1013 
1015  value_type fastAccessDx(int i) const {
1016  using std::pow; using std::log;
1017  return if_then_else( c == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val())) );
1018  }
1019 
1020  protected:
1021 
1022  const ConstT& c;
1023  const T2& expr2;
1024  };
1025 
1026  //
1027  // Specialization for scalar types using ternary operator
1028  //
1029 
1030  template <typename T1, typename T2>
1031  class PowerOp< T1, T2, false, false, ExprSpecDefault, false > :
1032  public Expr< PowerOp< T1, T2, false, false, ExprSpecDefault, false > > {
1033  public:
1034 
1035  typedef typename std::remove_cv<T1>::type ExprT1;
1036  typedef typename std::remove_cv<T2>::type ExprT2;
1037  typedef typename ExprT1::value_type value_type_1;
1038  typedef typename ExprT2::value_type value_type_2;
1039  typedef typename Sacado::Promote<value_type_1,
1040  value_type_2>::type value_type;
1041 
1042  typedef typename ExprT1::scalar_type scalar_type_1;
1043  typedef typename ExprT2::scalar_type scalar_type_2;
1044  typedef typename Sacado::Promote<scalar_type_1,
1045  scalar_type_2>::type scalar_type;
1046 
1047  typedef ExprSpecDefault expr_spec_type;
1048 
1050  PowerOp(const T1& expr1_, const T2& expr2_) :
1051  expr1(expr1_), expr2(expr2_) {}
1052 
1054  int size() const {
1055  const int sz1 = expr1.size(), sz2 = expr2.size();
1056  return sz1 > sz2 ? sz1 : sz2;
1057  }
1058 
1060  bool hasFastAccess() const {
1061  return expr1.hasFastAccess() && expr2.hasFastAccess();
1062  }
1063 
1065  value_type val() const {
1066  using std::pow;
1067  return pow(expr1.val(), expr2.val());
1068  }
1069 
1071  value_type dx(int i) const {
1072  using std::pow; using std::log;
1073  const int sz1 = expr1.size(), sz2 = expr2.size();
1074  if (sz1 > 0 && sz2 > 0)
1075  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()));
1076  else if (sz1 > 0)
1077  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1078  // It seems less accurate and caused convergence problems in some codes
1079  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()));
1080  else
1081  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
1082  }
1083 
1085  value_type fastAccessDx(int i) const {
1086  using std::pow; using std::log;
1087  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()));
1088  }
1089 
1090  protected:
1091 
1092  const T1& expr1;
1093  const T2& expr2;
1094 
1095  };
1096 
1097  template <typename T1, typename T2>
1098  class PowerOp< T1, T2, false, true, ExprSpecDefault, false >
1099  : public Expr< PowerOp< T1, T2, false, true, ExprSpecDefault, false > > {
1100  public:
1101 
1102  typedef typename std::remove_cv<T1>::type ExprT1;
1103  typedef T2 ConstT;
1104  typedef typename ExprT1::value_type value_type;
1105  typedef typename ExprT1::scalar_type scalar_type;
1106 
1107  typedef ExprSpecDefault expr_spec_type;
1108 
1110  PowerOp(const T1& expr1_, const ConstT& c_) :
1111  expr1(expr1_), c(c_) {}
1112 
1114  int size() const {
1115  return expr1.size();
1116  }
1117 
1119  bool hasFastAccess() const {
1120  return expr1.hasFastAccess();
1121  }
1122 
1124  value_type val() const {
1125  using std::pow;
1126  return pow(expr1.val(), c);
1127  }
1128 
1130  value_type dx(int i) const {
1131  using std::pow;
1132  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1133  // It seems less accurate and caused convergence problems in some codes
1134  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.dx(i)/expr1.val()*pow(expr1.val(),c));
1135  }
1136 
1138  value_type fastAccessDx(int i) const {
1139  using std::pow;
1140  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1141  // It seems less accurate and caused convergence problems in some codes
1142  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c));
1143  }
1144 
1145  protected:
1146 
1147  const T1& expr1;
1148  const ConstT& c;
1149  };
1150 
1151  template <typename T1, typename T2>
1152  class PowerOp< T1, T2, true, false, ExprSpecDefault, false >
1153  : public Expr< PowerOp< T1, T2, true, false, ExprSpecDefault, false > > {
1154  public:
1155 
1156  typedef typename std::remove_cv<T2>::type ExprT2;
1157  typedef T1 ConstT;
1158  typedef typename ExprT2::value_type value_type;
1159  typedef typename ExprT2::scalar_type scalar_type;
1160 
1161  typedef ExprSpecDefault expr_spec_type;
1162 
1164  PowerOp(const ConstT& c_, const T2& expr2_) :
1165  c(c_), expr2(expr2_) {}
1166 
1168  int size() const {
1169  return expr2.size();
1170  }
1171 
1173  bool hasFastAccess() const {
1174  return expr2.hasFastAccess();
1175  }
1176 
1178  value_type val() const {
1179  using std::pow;
1180  return pow(c, expr2.val());
1181  }
1182 
1184  value_type dx(int i) const {
1185  using std::pow; using std::log;
1186  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c)*pow(c,expr2.val()));
1187  }
1188 
1190  value_type fastAccessDx(int i) const {
1191  using std::pow; using std::log;
1192  return c == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c)*pow(c,expr2.val()));
1193  }
1194 
1195  protected:
1196 
1197  const ConstT& c;
1198  const T2& expr2;
1199  };
1200 
1201  template <typename T1, typename T2>
1204  pow (const T1& expr1, const T2& expr2)
1205  {
1206  typedef PowerOp< typename Expr<T1>::derived_type,
1207  typename Expr<T2>::derived_type,
1208  false, false, typename T1::expr_spec_type > expr_t;
1209 
1210  return expr_t(expr1.derived(), expr2.derived());
1211  }
1212 
1213  template <typename T>
1215  PowerOp< typename T::value_type, typename Expr<T>::derived_type,
1216  true, false, typename T::expr_spec_type >
1217  pow (const typename T::value_type& c,
1218  const Expr<T>& expr)
1219  {
1220  typedef typename T::value_type ConstT;
1221  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1222  true, false, typename T::expr_spec_type > expr_t;
1223 
1224  return expr_t(c, expr.derived());
1225  }
1226 
1227  template <typename T>
1229  PowerOp< typename Expr<T>::derived_type, typename T::value_type,
1230  false, true, typename T::expr_spec_type >
1231  pow (const Expr<T>& expr,
1232  const typename T::value_type& c)
1233  {
1234  typedef typename T::value_type ConstT;
1235  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1236  false, true, typename T::expr_spec_type > expr_t;
1237 
1238  return expr_t(expr.derived(), c);
1239  }
1240 
1241  template <typename T>
1244  pow (const typename T::scalar_type& c,
1245  const Expr<T>& expr)
1246  {
1247  typedef typename T::scalar_type ConstT;
1248  typedef PowerOp< ConstT, typename Expr<T>::derived_type,
1249  true, false, typename T::expr_spec_type > expr_t;
1250 
1251  return expr_t(c, expr.derived());
1252  }
1253 
1254  template <typename T>
1257  pow (const Expr<T>& expr,
1258  const typename T::scalar_type& c)
1259  {
1260  typedef typename T::scalar_type ConstT;
1261  typedef PowerOp< typename Expr<T>::derived_type, ConstT,
1262  false, true, typename T::expr_spec_type > expr_t;
1263 
1264  return expr_t(expr.derived(), c);
1265  }
1266 
1267  template <typename T1, typename T2, bool c1, bool c2, typename E>
1268  struct ExprLevel< PowerOp< T1, T2, c1, c2, E > > {
1269  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1270  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1271  static constexpr unsigned value =
1272  value_1 >= value_2 ? value_1 : value_2;
1273  };
1274 
1275  template <typename T1, typename T2, bool c1, bool c2, typename E>
1276  struct IsFadExpr< PowerOp< T1, T2, c1, c2, E > > {
1277  static constexpr unsigned value = true;
1278  };
1279 
1280  }
1281  }
1282 
1283  template <typename T1, typename T2, bool c1, bool c2, typename E>
1284  struct IsExpr< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1285  static constexpr bool value = true;
1286  };
1287 
1288  template <typename T1, typename T2, bool c1, bool c2, typename E>
1289  struct BaseExprType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1290  typedef typename BaseExprType<T1>::type base_expr_1;
1291  typedef typename BaseExprType<T2>::type base_expr_2;
1292  typedef typename Sacado::Promote<base_expr_1,
1293  base_expr_2>::type type;
1294  };
1295 
1296  template <typename T1, typename T2, bool c1, bool c2, typename E>
1297  struct IsSimdType< Fad::Exp::PowerOp< T1, T2, c1, c2, E > > {
1298  static const bool value =
1299  IsSimdType< typename Fad::Exp::PowerOp< T1, T2, c1, c2, E >::value_type >::value;
1300  };
1301 
1302 }
1303 
1304 //--------------------------if_then_else operator -----------------------
1305 // Can't use the above macros because it is a ternary operator (sort of).
1306 // Also, relies on C++11
1307 
1308 namespace Sacado {
1309  namespace Fad {
1310  namespace Exp {
1311 
1312  template <typename CondT, typename T1, typename T2,
1313  bool is_const_T1, bool is_const_T2,
1314  typename ExprSpec>
1315  class IfThenElseOp {};
1316 
1317  template <typename CondT, typename T1, typename T2>
1319  public Expr< IfThenElseOp< CondT, T1, T2, false, false, ExprSpecDefault > > {
1320 
1321  public:
1322 
1323  typedef typename std::remove_cv<T1>::type ExprT1;
1324  typedef typename std::remove_cv<T2>::type ExprT2;
1325  typedef typename ExprT1::value_type value_type_1;
1326  typedef typename ExprT2::value_type value_type_2;
1327  typedef typename Sacado::Promote<value_type_1,
1329 
1330  typedef typename ExprT1::scalar_type scalar_type_1;
1331  typedef typename ExprT2::scalar_type scalar_type_2;
1332  typedef typename Sacado::Promote<scalar_type_1,
1334 
1336 
1338  IfThenElseOp(const CondT& cond_, const T1& expr1_, const T2& expr2_) :
1339  cond(cond_), expr1(expr1_), expr2(expr2_) {}
1340 
1342  int size() const {
1343  int sz1 = expr1.size(), sz2 = expr2.size();
1344  return sz1 > sz2 ? sz1 : sz2;
1345  }
1346 
1348  bool hasFastAccess() const {
1349  return expr1.hasFastAccess() && expr2.hasFastAccess();
1350  }
1351 
1353  value_type val() const {
1354  return if_then_else( cond, expr1.val(), expr2.val() );
1355  }
1356 
1358  value_type dx(int i) const {
1359  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
1360  }
1361 
1363  value_type fastAccessDx(int i) const {
1364  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
1365  }
1366 
1367  protected:
1368 
1369  const CondT& cond;
1370  const T1& expr1;
1371  const T2& expr2;
1372 
1373  };
1374 
1375  template <typename CondT, typename T1, typename T2>
1376  class IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > :
1377  public Expr< IfThenElseOp< CondT, T1, T2, false, true, ExprSpecDefault > > {
1378 
1379  public:
1380 
1381  typedef typename std::remove_cv<T1>::type ExprT1;
1382  typedef T2 ConstT;
1383  typedef typename ExprT1::value_type value_type;
1384  typedef typename ExprT1::scalar_type scalar_type;
1385 
1387 
1389  IfThenElseOp(const CondT& cond_, const T1& expr1_, const ConstT& c_) :
1390  cond(cond_), expr1(expr1_), c(c_) {}
1391 
1393  int size() const {
1394  return expr1.size();
1395  }
1396 
1398  bool hasFastAccess() const {
1399  return expr1.hasFastAccess();
1400  }
1401 
1403  value_type val() const {
1404  return if_then_else( cond, expr1.val(), c );
1405  }
1406 
1408  value_type dx(int i) const {
1409  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
1410  }
1411 
1413  value_type fastAccessDx(int i) const {
1414  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
1415  }
1416 
1417  protected:
1418 
1419  const CondT& cond;
1420  const T1& expr1;
1421  const ConstT& c;
1422  };
1423 
1424  template <typename CondT, typename T1, typename T2>
1425  class IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > :
1426  public Expr< IfThenElseOp< CondT, T1, T2, true, false, ExprSpecDefault > > {
1427 
1428  public:
1429 
1430  typedef typename std::remove_cv<T2>::type ExprT2;
1431  typedef T1 ConstT;
1432  typedef typename ExprT2::value_type value_type;
1433  typedef typename ExprT2::scalar_type scalar_type;
1434 
1436 
1438  IfThenElseOp(const CondT& cond_, const ConstT& c_, const T2& expr2_) :
1439  cond(cond_), c(c_), expr2(expr2_) {}
1440 
1442  int size() const {
1443  return expr2.size();
1444  }
1445 
1447  bool hasFastAccess() const {
1448  return expr2.hasFastAccess();
1449  }
1450 
1452  value_type val() const {
1453  return if_then_else( cond, c, expr2.val() );
1454  }
1455 
1457  value_type dx(int i) const {
1458  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
1459  }
1460 
1462  value_type fastAccessDx(int i) const {
1463  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
1464  }
1465 
1466  protected:
1467 
1468  const CondT& cond;
1469  const ConstT& c;
1470  const T2& expr2;
1471  };
1472 
1473  template <typename CondT, typename T1, typename T2>
1477  IfThenElseOp< CondT,
1478  typename Expr<T1>::derived_type,
1479  typename Expr<T2>::derived_type,
1480  false, false,
1481  typename T1::expr_spec_type >
1482  >::type
1483  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
1484  {
1486  typename Expr<T2>::derived_type,
1487  false, false, typename T1::expr_spec_type > expr_t;
1488 
1489  return expr_t(cond, expr1.derived(), expr2.derived());
1490  }
1491 
1492  template <typename CondT, typename T>
1494  IfThenElseOp< CondT, typename T::value_type, typename Expr<T>::derived_type,
1495  true, false, typename T::expr_spec_type >
1496  if_then_else (const CondT& cond, const typename T::value_type& c,
1497  const Expr<T>& expr)
1498  {
1499  typedef typename T::value_type ConstT;
1501  true, false, typename T::expr_spec_type > expr_t;
1502 
1503  return expr_t(cond, c, expr.derived());
1504  }
1505 
1506  template <typename CondT, typename T>
1508  IfThenElseOp< CondT, typename Expr<T>::derived_type, typename T::value_type,
1509  false, true, typename T::expr_spec_type >
1510  if_then_else (const CondT& cond, const Expr<T>& expr,
1511  const typename T::value_type& c)
1512  {
1513  typedef typename T::value_type ConstT;
1515  false, true, typename T::expr_spec_type > expr_t;
1516 
1517  return expr_t(cond, expr.derived(), c);
1518  }
1519 
1520  template <typename CondT, typename T>
1522  typename mpl::disable_if<
1523  mpl::is_same< typename T::value_type,
1524  typename T::scalar_type >,
1525  IfThenElseOp< CondT, typename T::scalar_type,
1526  typename Expr<T>::derived_type,
1527  true, false, typename T::expr_spec_type >
1528  >::type
1529  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
1530  const Expr<T>& expr)
1531  {
1532  typedef typename T::scalar_type ConstT;
1534  true, false, typename T::expr_spec_type > expr_t;
1535 
1536  return expr_t(cond, c, expr.derived());
1537  }
1538 
1539  template <typename CondT, typename T>
1541  typename mpl::disable_if<
1542  mpl::is_same< typename T::value_type,
1543  typename T::scalar_type >,
1544  IfThenElseOp< CondT, typename Expr<T>::derived_type,
1545  typename T::scalar_type,
1546  false, true, typename T::expr_spec_type >
1547  >::type
1548  if_then_else (const CondT& cond, const Expr<T>& expr,
1549  const typename Expr<T>::scalar_type& c)
1550  {
1551  typedef typename T::scalar_type ConstT;
1553  false, true, typename T::expr_spec_type > expr_t;
1554 
1555  return expr_t(cond, expr.derived(), c);
1556  }
1557 
1558  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1559  typename E>
1560  struct ExprLevel< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1561  static constexpr unsigned value_1 = ExprLevel<T1>::value;
1562  static constexpr unsigned value_2 = ExprLevel<T2>::value;
1563  static constexpr unsigned value =
1564  value_1 >= value_2 ? value_1 : value_2;
1565  };
1566 
1567  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1568  typename E>
1569  struct IsFadExpr< IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1570  static constexpr unsigned value = true;
1571  };
1572 
1573  }
1574  }
1575 
1576  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1577  typename E>
1578  struct IsExpr< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1579  static constexpr bool value = true;
1580  };
1581 
1582  template <typename CondT, typename T1, typename T2, bool c1, bool c2,
1583  typename E>
1584  struct BaseExprType< Fad::Exp::IfThenElseOp< CondT, T1, T2, c1, c2, E > > {
1587  typedef typename Sacado::Promote<base_expr_1,
1589  };
1590 }
1591 
1592 #undef FAD_BINARYOP_MACRO
1593 
1594 //-------------------------- Relational Operators -----------------------
1595 
1596 namespace Sacado {
1597  namespace Fad {
1598  namespace Exp {
1599  namespace Impl {
1600  // Helper trait to determine return type of logical comparison operations
1601  // (==, !=, ...), usually bool but maybe something else for SIMD types.
1602  // Need to use SFINAE to restrict to types that define == operator in the
1603  // conditional overloads below, otherwise instantiating ConditionaReturnType
1604  // may fail during overload resolution.
1605  template <typename T1, typename T2 = T1,
1608 
1609  template <typename T1, typename T2>
1611  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
1612  };
1613  }
1614  }
1615  }
1616 }
1617 
1618 #define FAD_RELOP_MACRO(OP) \
1619 namespace Sacado { \
1620  namespace Fad { \
1621  namespace Exp { \
1622  template <typename T1, typename T2> \
1623  KOKKOS_INLINE_FUNCTION \
1624  typename mpl::enable_if_c< \
1625  IsFadExpr<T1>::value && IsFadExpr<T2>::value && \
1626  ExprLevel<T1>::value == ExprLevel<T2>::value, \
1627  typename Impl::ConditionalReturnType<typename T1::value_type, \
1628  typename T2::value_type>::type \
1629  >::type \
1630  operator OP (const T1& expr1, const T2& expr2) \
1631  { \
1632  return expr1.derived().val() OP expr2.derived().val(); \
1633  } \
1634  \
1635  template <typename T2> \
1636  KOKKOS_INLINE_FUNCTION \
1637  typename Impl::ConditionalReturnType<typename T2::value_type>::type \
1638  operator OP (const typename T2::value_type& a, \
1639  const Expr<T2>& expr2) \
1640  { \
1641  return a OP expr2.derived().val(); \
1642  } \
1643  \
1644  template <typename T1> \
1645  KOKKOS_INLINE_FUNCTION \
1646  typename Impl::ConditionalReturnType<typename T1::value_type>::type \
1647  operator OP (const Expr<T1>& expr1, \
1648  const typename T1::value_type& b) \
1649  { \
1650  return expr1.derived().val() OP b; \
1651  } \
1652  } \
1653  } \
1654 } \
1655 
1656 FAD_RELOP_MACRO(==)
1657 FAD_RELOP_MACRO(!=)
1658 FAD_RELOP_MACRO(<)
1659 FAD_RELOP_MACRO(>)
1660 FAD_RELOP_MACRO(<=)
1661 FAD_RELOP_MACRO(>=)
1662 FAD_RELOP_MACRO(<<=)
1663 FAD_RELOP_MACRO(>>=)
1664 FAD_RELOP_MACRO(&)
1665 FAD_RELOP_MACRO(|)
1666 
1667 #undef FAD_RELOP_MACRO
1668 
1669 namespace Sacado {
1670 
1671  namespace Fad {
1672  namespace Exp {
1673 
1674  template <typename ExprT>
1676  bool operator ! (const Expr<ExprT>& expr)
1677  {
1678  return ! expr.derived().val();
1679  }
1680 
1681  } // namespace Exp
1682  } // namespace Fad
1683 
1684 } // namespace Sacado
1685 
1686 //-------------------------- Boolean Operators -----------------------
1687 namespace Sacado {
1688 
1689  namespace Fad {
1690  namespace Exp {
1691 
1692  template <typename T>
1694  bool toBool(const Expr<T>& xx) {
1695  const typename Expr<T>::derived_type& x = xx.derived();
1696  bool is_zero = (x.val() == 0.0);
1697  for (int i=0; i<x.size(); i++)
1698  is_zero = is_zero && (x.dx(i) == 0.0);
1699  return !is_zero;
1700  }
1701 
1702  } // namespace Exp
1703  } // namespace Fad
1704 
1705 } // namespace Sacado
1706 
1707 #define FAD_BOOL_MACRO(OP) \
1708 namespace Sacado { \
1709  namespace Fad { \
1710  namespace Exp { \
1711  template <typename T1, typename T2> \
1712  KOKKOS_INLINE_FUNCTION \
1713  bool \
1714  operator OP (const Expr<T1>& expr1, \
1715  const Expr<T2>& expr2) \
1716  { \
1717  return toBool(expr1) OP toBool(expr2); \
1718  } \
1719  \
1720  template <typename T2> \
1721  KOKKOS_INLINE_FUNCTION \
1722  bool \
1723  operator OP (const typename Expr<T2>::value_type& a, \
1724  const Expr<T2>& expr2) \
1725  { \
1726  return a OP toBool(expr2); \
1727  } \
1728  \
1729  template <typename T1> \
1730  KOKKOS_INLINE_FUNCTION \
1731  bool \
1732  operator OP (const Expr<T1>& expr1, \
1733  const typename Expr<T1>::value_type& b) \
1734  { \
1735  return toBool(expr1) OP b; \
1736  } \
1737  } \
1738  } \
1739 }
1740 
1741 FAD_BOOL_MACRO(&&)
1742 FAD_BOOL_MACRO(||)
1743 
1744 #undef FAD_BOOL_MACRO
1745 
1746 //-------------------------- I/O Operators -----------------------
1747 
1748 namespace Sacado {
1749 
1750  namespace Fad {
1751  namespace Exp {
1752 
1753  template <typename T>
1754  std::ostream& operator << (std::ostream& os, const Expr<T>& xx) {
1755  const typename Expr<T>::derived_type& x = xx.derived();
1756  os << x.val() << " [";
1757 
1758  for (int i=0; i< x.size(); i++) {
1759  os << " " << x.dx(i);
1760  }
1761 
1762  os << " ]";
1763  return os;
1764  }
1765 
1766  } // namespace Exp
1767  } // namespace Fad
1768 
1769 } // namespace Sacado
1770 
1771 #if defined(HAVE_SACADO_KOKKOSCORE)
1772 
1773 //-------------------------- Atomic Operators -----------------------
1774 
1775 namespace Sacado {
1776 
1777  namespace Fad {
1778  namespace Exp {
1779 
1780  // Overload of Kokkos::atomic_add for Fad types.
1781  template <typename S>
1783  void atomic_add(GeneralFad<S>* dst, const GeneralFad<S>& x) {
1784  using Kokkos::atomic_add;
1785 
1786  const int xsz = x.size();
1787  const int sz = dst->size();
1788 
1789  // We currently cannot handle resizing since that would need to be
1790  // done atomically.
1791  if (xsz > sz)
1792  Kokkos::abort(
1793  "Sacado error: Fad resize within atomic_add() not supported!");
1794 
1795  if (xsz != sz && sz > 0 && xsz > 0)
1796  Kokkos::abort(
1797  "Sacado error: Fad assignment of incompatiable sizes!");
1798 
1799 
1800  if (sz > 0 && xsz > 0) {
1801  SACADO_FAD_DERIV_LOOP(i,sz)
1802  atomic_add(&(dst->fastAccessDx(i)), x.fastAccessDx(i));
1803  }
1805  atomic_add(&(dst->val()), x.val());
1806  }
1807 
1808  } // namespace Exp
1809  } // namespace Fad
1810 
1811 } // namespace Sacado
1812 
1813 #endif // HAVE_SACADO_KOKKOSCORE
1814 
1815 #endif // SACADO_FAD_OPS_HPP
Wrapper for a generic expression template.
cbrt(expr.val())
#define FAD_BOOL_MACRO(OP)
expr expr SinOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
asinh(expr.val())
#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 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)
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
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.
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 Atan2Op
expr val()
#define KOKKOS_INLINE_FUNCTION
KOKKOS_INLINE_FUNCTION T safe_sqrt(const T &x)
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())
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
#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 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