Stokhos 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_MP_Vector.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef SACADO_FAD_EXP_MP_VECTOR_HPP
11 #define SACADO_FAD_EXP_MP_VECTOR_HPP
12 
13 #include "Sacado_MP_Vector.hpp"
14 
15 namespace Sacado {
16 
17  namespace Fad {
18  namespace Exp {
19 
21  class ExprSpecMPVector {};
22 
24 
29  template <typename T>
30  class Extender<
31  T,
32  typename std::enable_if<
33  Sacado::is_mp_vector<typename T::value_type>::value >::type
34  > : public T
35  {
36  public:
37 
38  typedef typename T::value_type value_type;
39  typedef typename value_type::value_type val_type;
40 
41  // Define expression template specialization
43 
44  // Bring in constructors
45  using T::T;
46 
47  // Bring in methods we are overloading
48  using T::val;
49  using T::dx;
50  using T::fastAccessDx;
51 
53  KOKKOS_INLINE_FUNCTION
54  const val_type& val(int j) const { return T::val().fastAccessCoeff(j); }
55 
57  KOKKOS_INLINE_FUNCTION
58  val_type& val(int j) { return T::val().fastAccessCoeff(j); }
59 
61  KOKKOS_INLINE_FUNCTION
62  val_type dx(int i, int j) const {
63  return this->size() ? this->dx_[i].fastAccessCoeff(j) : val_type(0.0);
64  }
65 
67  KOKKOS_INLINE_FUNCTION
68  val_type& fastAccessDx(int i, int j) {
69  return this->dx_[i].fastAccessCoeff(j);
70  }
71 
73  KOKKOS_INLINE_FUNCTION
74  const val_type& fastAccessDx(int i, int j) const {
75  return this->dx_[i].fastAccessCoeff(j);
76  }
77 
78  };
79 
81  template <typename DstType>
82  class ExprAssign<
83  DstType,
84  typename std::enable_if<
85  std::is_same< typename DstType::expr_spec_type, ExprSpecMPVector >::value
86  >::type
87  > {
88  public:
89 
91  typedef typename DstType::value_type value_type;
92 
93  // MP::Vector size (assuming static, because that's all we care about)
94  static const int VecNum = Sacado::StaticSize<value_type>::value;
95 
97  template <typename SrcType>
98  KOKKOS_INLINE_FUNCTION
99  static void assign_equal(DstType& dst, const SrcType& x)
100  {
101  const int xsz = x.size();
102 
103  if (xsz != dst.size())
104  dst.resizeAndZero(xsz);
105 
106  const int sz = dst.size();
107 
108  // For ViewStorage, the resize above may not in fact resize the
109  // derivative array, so it is possible that sz != xsz at this point.
110  // The only valid use case here is sz > xsz == 0, so we use sz in the
111  // assignment below
112 
113  if (sz) {
114  if (x.hasFastAccess()) {
115  SACADO_FAD_DERIV_LOOP(i,sz)
116  for (int j=0; j<VecNum; ++j)
117  dst.fastAccessDx(i,j) = x.fastAccessDx(i,j);
118  }
119  else
120  SACADO_FAD_DERIV_LOOP(i,sz)
121  for (int j=0; j<VecNum; ++j)
122  dst.fastAccessDx(i,j) = x.dx(i,j);
123  }
124 
125  for (int j=0; j<VecNum; ++j)
126  dst.val(j) = x.val(j);
127  }
128 
130  template <typename SrcType>
131  KOKKOS_INLINE_FUNCTION
132  static void assign_plus_equal(DstType& dst, const SrcType& x)
133  {
134  const int xsz = x.size(), sz = dst.size();
135 
136 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
137  if ((xsz != sz) && (xsz != 0) && (sz != 0))
138  throw "Fad Error: Attempt to assign with incompatible sizes";
139 #endif
140 
141  if (xsz) {
142  if (sz) {
143  if (x.hasFastAccess())
144  SACADO_FAD_DERIV_LOOP(i,sz)
145  for (int j=0; j<VecNum; ++j)
146  dst.fastAccessDx(i,j) += x.fastAccessDx(i,j);
147  else
148  for (int i=0; i<sz; ++i)
149  for (int j=0; j<VecNum; ++j)
150  dst.fastAccessDx(i,j) += x.dx(i,j);
151  }
152  else {
153  dst.resizeAndZero(xsz);
154  if (x.hasFastAccess())
155  SACADO_FAD_DERIV_LOOP(i,xsz)
156  for (int j=0; j<VecNum; ++j)
157  dst.fastAccessDx(i,j) = x.fastAccessDx(i,j);
158  else
159  SACADO_FAD_DERIV_LOOP(i,xsz)
160  for (int j=0; j<VecNum; ++j)
161  dst.fastAccessDx(i,j) = x.dx(i,j);
162  }
163  }
164 
165  for (int j=0; j<VecNum; ++j)
166  dst.val(j) += x.val(j);
167  }
168 
170  template <typename SrcType>
171  KOKKOS_INLINE_FUNCTION
172  static void assign_minus_equal(DstType& dst, const SrcType& x)
173  {
174  const int xsz = x.size(), sz = dst.size();
175 
176 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
177  if ((xsz != sz) && (xsz != 0) && (sz != 0))
178  throw "Fad Error: Attempt to assign with incompatible sizes";
179 #endif
180 
181  if (xsz) {
182  if (sz) {
183  if (x.hasFastAccess())
184  SACADO_FAD_DERIV_LOOP(i,sz)
185  for (int j=0; j<VecNum; ++j)
186  dst.fastAccessDx(i,j) -= x.fastAccessDx(i,j);
187  else
188  SACADO_FAD_DERIV_LOOP(i,sz)
189  for (int j=0; j<VecNum; ++j)
190  dst.fastAccessDx(i,j) -= x.dx(i,j);
191  }
192  else {
193  dst.resizeAndZero(xsz);
194  if (x.hasFastAccess())
195  SACADO_FAD_DERIV_LOOP(i,xsz)
196  for (int j=0; j<VecNum; ++j)
197  dst.fastAccessDx(i,j) = -x.fastAccessDx(i,j);
198  else
199  SACADO_FAD_DERIV_LOOP(i,xsz)
200  for (int j=0; j<VecNum; ++j)
201  dst.fastAccessDx(i,j) = -x.dx(i,j);
202  }
203  }
204 
205  for (int j=0; j<VecNum; ++j)
206  dst.val(j) -= x.val(j);
207  }
208 
210  template <typename SrcType>
211  KOKKOS_INLINE_FUNCTION
212  static void assign_times_equal(DstType& dst, const SrcType& x)
213  {
214  const int xsz = x.size(), sz = dst.size();
215  const value_type xval = x.val();
216  const value_type v = dst.val();
217 
218 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
219  if ((xsz != sz) && (xsz != 0) && (sz != 0))
220  throw "Fad Error: Attempt to assign with incompatible sizes";
221 #endif
222 
223  if (xsz) {
224  if (sz) {
225  if (x.hasFastAccess())
226  SACADO_FAD_DERIV_LOOP(i,sz)
227  for (int j=0; j<VecNum; ++j)
228  dst.fastAccessDx(i) = v.fastAccessCoeff(j)*x.fastAccessDx(i,j) + dst.fastAccessDx(i,j)*xval.fastAccessCoeff(j);
229  else
230  SACADO_FAD_DERIV_LOOP(i,sz)
231  for (int j=0; j<VecNum; ++j)
232  dst.fastAccessDx(i) = v.fastAccessCoeff(j)*x.dx(i,j) + dst.fastAccessDx(i,j)*xval.fastAccessCoeff(j);
233  }
234  else {
235  dst.resizeAndZero(xsz);
236  if (x.hasFastAccess())
237  SACADO_FAD_DERIV_LOOP(i,xsz)
238  for (int j=0; j<VecNum; ++j)
239  dst.fastAccessDx(i,j) = v.fastAccessCoeff(j)*x.fastAccessDx(i,j);
240  else
241  SACADO_FAD_DERIV_LOOP(i,xsz)
242  for (int j=0; j<VecNum; ++j)
243  dst.fastAccessDx(i,j) = v.fastAccessCoeff(j)*x.dx(i,j);
244  }
245  }
246  else {
247  if (sz) {
248  SACADO_FAD_DERIV_LOOP(i,sz)
249  for (int j=0; j<VecNum; ++j)
250  dst.fastAccessDx(i,j) *= xval.fastAccessCoeff(j);
251  }
252  }
253 
254  for (int j=0; j<VecNum; ++j)
255  dst.val(j) *= xval.fastAccessCoeff(j);
256  }
257 
259  template <typename SrcType>
260  KOKKOS_INLINE_FUNCTION
261  static void assign_divide_equal(DstType& dst, const SrcType& x)
262  {
263  const int xsz = x.size(), sz = dst.size();
264  const value_type xval = x.val();
265  const value_type v = dst.val();
266 
267 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
268  if ((xsz != sz) && (xsz != 0) && (sz != 0))
269  throw "Fad Error: Attempt to assign with incompatible sizes";
270 #endif
271 
272  if (xsz) {
273  const value_type xval2 = xval*xval;
274  if (sz) {
275  if (x.hasFastAccess())
276  SACADO_FAD_DERIV_LOOP(i,sz)
277  for (int j=0; j<VecNum; ++j)
278  dst.fastAccessDx(i,j) =
279  ( dst.fastAccessDx(i,j)*xval.fastAccessCoeff(j) - v.fastAccessCoeff(j)*x.fastAccessDx(i,j) ) / xval2.fastAccessCoeff(j);
280  else
281  SACADO_FAD_DERIV_LOOP(i,sz)
282  for (int j=0; j<VecNum; ++j)
283  dst.fastAccessDx(i,j) =
284  ( dst.fastAccessDx(i,j)*xval.fastAccessCoeff(j) - v.fastAccessCoeff(j)*x.dx(i,j) ) / xval2.fastAccessCoeff(j);
285  }
286  else {
287  dst.resizeAndZero(xsz);
288  if (x.hasFastAccess())
289  SACADO_FAD_DERIV_LOOP(i,xsz)
290  for (int j=0; j<VecNum; ++j)
291  dst.fastAccessDx(i,j) = - v.fastAccessCoeff(j)*x.fastAccessDx(i,j) / xval2.fastAccessCoeff(j);
292  else
293  SACADO_FAD_DERIV_LOOP(i,xsz)
294  for (int j=0; j<VecNum; ++j)
295  dst.fastAccessDx(i,j) = -v.fastAccessCoeff(j)*x.dx(i,j) / xval2.fastAccessCoeff(j);
296  }
297  }
298  else {
299  if (sz) {
300  SACADO_FAD_DERIV_LOOP(i,sz)
301  for (int j=0; j<VecNum; ++j)
302  dst.fastAccessDx(i,j) /= xval.fastAccessCoeff(j);
303  }
304  }
305 
306  for (int j=0; j<VecNum; ++j)
307  dst.val(j) /= xval.fastAccessCoeff(j);
308  }
309 
310  };
311 
316  template <typename DstType>
317  class ExprAssign<
318  DstType,
319  typename std::enable_if<
320  Sacado::IsStaticallySized<DstType>::value &&
321  std::is_same< typename DstType::expr_spec_type, ExprSpecMPVector >::value
322  >::type
323  > {
324  public:
325 
327  typedef typename DstType::value_type value_type;
328 
329  // MP::Vector size (assuming static, because that's all we care about)
330  static const int VecNum = Sacado::StaticSize<value_type>::value;
331 
333  template <typename SrcType>
334  KOKKOS_INLINE_FUNCTION
335  static void assign_equal(DstType& dst, const SrcType& x)
336  {
337  const int sz = dst.size();
338  SACADO_FAD_DERIV_LOOP(i,sz)
339  for (int j=0; j<VecNum; ++j)
340  dst.fastAccessDx(i,j) = x.fastAccessDx(i,j);
341  for (int j=0; j<VecNum; ++j)
342  dst.val(j) = x.val(j);
343  }
344 
346  template <typename SrcType>
347  KOKKOS_INLINE_FUNCTION
348  static void assign_plus_equal(DstType& dst, const SrcType& x)
349  {
350  const int sz = dst.size();
351  SACADO_FAD_DERIV_LOOP(i,sz)
352  for (int j=0; j<VecNum; ++j)
353  dst.fastAccessDx(i,j) += x.fastAccessDx(i,j);
354  for (int j=0; j<VecNum; ++j)
355  dst.val(j) += x.val(j);
356  }
357 
359  template <typename SrcType>
360  KOKKOS_INLINE_FUNCTION
361  static void assign_minus_equal(DstType& dst, const SrcType& x)
362  {
363  const int sz = dst.size();
364  SACADO_FAD_DERIV_LOOP(i,sz)
365  for (int j=0; j<VecNum; ++j)
366  dst.fastAccessDx(i,j) -= x.fastAccessDx(i,j);
367  for (int j=0; j<VecNum; ++j)
368  dst.val(j) -= x.val(j);
369  }
370 
372  template <typename SrcType>
373  KOKKOS_INLINE_FUNCTION
374  static void assign_times_equal(DstType& dst, const SrcType& x)
375  {
376  const int sz = dst.size();
377  const value_type xval = x.val();
378  const value_type v = dst.val();
379  SACADO_FAD_DERIV_LOOP(i,sz)
380  for (int j=0; j<VecNum; ++j)
381  dst.fastAccessDx(i,j) = v.fastAccessCoeff(j)*x.fastAccessDx(i,j) + dst.fastAccessDx(i,j)*xval.fastAccessCoeff(j);
382  for (int j=0; j<VecNum; ++j)
383  dst.val(j) *= xval.fastAccessCoeff(j);
384  }
385 
387  template <typename SrcType>
388  KOKKOS_INLINE_FUNCTION
389  static void assign_divide_equal(DstType& dst, const SrcType& x)
390  {
391  const int sz = dst.size();
392  const value_type xval = x.val();
393  const value_type xval2 = xval*xval;
394  const value_type v = dst.val();
395  SACADO_FAD_DERIV_LOOP(i,sz)
396  for (int j=0; j<VecNum; ++j)
397  dst.fastAccessDx(i,j) =
398  ( dst.fastAccessDx(i,j)*xval.fastAccessCoeff(j) - v.fastAccessCoeff(j)*x.fastAccessDx(i,j) )/ xval2.fastAccessCoeff(j);
399  for (int j=0; j<VecNum; ++j)
400  dst.val(j) /= xval.fastAccessCoeff(j);
401  }
402 
403  };
404 
405  } // namespace Exp
406  } // namespace Fad
407 
408 } // namespace Sacado
409 
410 // Specialize expression template operators to add similar extensions
411 #include "Sacado_Fad_Exp_Ops.hpp"
412 
413 #define FAD_UNARYOP_MACRO(OPNAME,OP,USING,MPVALUE,VALUE,DX,FASTACCESSDX) \
414 namespace Sacado { \
415  namespace Fad { \
416  namespace Exp { \
417  \
418  template <typename T> \
419  class OP< T,ExprSpecMPVector > : \
420  public Expr< OP< T,ExprSpecMPVector > > { \
421  public: \
422  \
423  typedef typename std::remove_cv<T>::type ExprT; \
424  typedef typename ExprT::value_type value_type; \
425  typedef typename ExprT::scalar_type scalar_type; \
426  \
427  typedef typename value_type::value_type val_type; \
428  \
429  typedef ExprSpecMPVector expr_spec_type; \
430  \
431  KOKKOS_INLINE_FUNCTION \
432  OP(const T& expr_) : expr(expr_) {} \
433  \
434  KOKKOS_INLINE_FUNCTION \
435  int size() const { return expr.size(); } \
436  \
437  KOKKOS_INLINE_FUNCTION \
438  bool hasFastAccess() const { \
439  return expr.hasFastAccess(); \
440  } \
441  \
442  KOKKOS_INLINE_FUNCTION \
443  value_type val() const { \
444  USING \
445  return MPVALUE; \
446  } \
447  \
448  KOKKOS_INLINE_FUNCTION \
449  val_type val(int j) const { \
450  USING \
451  return VALUE; \
452  } \
453  \
454  KOKKOS_INLINE_FUNCTION \
455  val_type dx(int i, int j) const { \
456  USING \
457  return DX; \
458  } \
459  \
460  KOKKOS_INLINE_FUNCTION \
461  val_type fastAccessDx(int i, int j) const { \
462  USING \
463  return FASTACCESSDX; \
464  } \
465  \
466  protected: \
467  \
468  const T& expr; \
469  }; \
470  \
471  } \
472  } \
473  \
474 }
475 
476 FAD_UNARYOP_MACRO(operator+,
477  UnaryPlusOp,
478  ;,
479  expr.val(),
480  expr.val(j),
481  expr.dx(i,j),
482  expr.fastAccessDx(i,j))
483 FAD_UNARYOP_MACRO(operator-,
485  ;,
486  -expr.val(),
487  -expr.val(j),
488  -expr.dx(i,j),
489  -expr.fastAccessDx(i,j))
492  using std::exp;,
493  exp(expr.val()),
494  exp(expr.val(j)),
495  exp(expr.val(j))*expr.dx(i,j),
496  exp(expr.val(j))*expr.fastAccessDx(i,j))
498  LogOp,
499  using std::log;,
500  log(expr.val()),
501  log(expr.val(j)),
502  expr.dx(i,j)/expr.val(j),
503  expr.fastAccessDx(i,j)/expr.val(j))
506  using std::log10; using std::log;,
507  log10(expr.val()),
508  log10(expr.val(j)),
509  expr.dx(i,j)/( log(value_type(10))*expr.val()),
510  expr.fastAccessDx(i,j) / ( log(value_type(10))*expr.val()))
513  using std::sqrt;,
514  sqrt(expr.val()),
515  sqrt(expr.val(j)),
516  expr.dx(i,j)/(value_type(2)* sqrt(expr.val())),
517  expr.fastAccessDx(i,j)/(value_type(2)* sqrt(expr.val())))
520  using std::cos; using std::sin;,
521  cos(expr.val()),
522  cos(expr.val(j)),
523  -expr.dx(i,j)* sin(expr.val()),
524  -expr.fastAccessDx(i,j)* sin(expr.val()))
527  using std::cos; using std::sin;,
528  sin(expr.val()),
529  sin(expr.val(j)),
530  expr.dx(i,j)* cos(expr.val()),
531  expr.fastAccessDx(i,j)* cos(expr.val()))
534  using std::tan;,
535  tan(expr.val()),
536  tan(expr.val(j)),
537  expr.dx(i,j)*
538  (value_type(1)+ tan(expr.val())* tan(expr.val())),
539  expr.fastAccessDx(i,j)*
540  (value_type(1)+ tan(expr.val())* tan(expr.val())))
543  using std::acos; using std::sqrt;,
544  acos(expr.val()),
545  acos(expr.val(j)),
546  -expr.dx(i,j)/ sqrt(value_type(1)-expr.val()*expr.val()),
547  -expr.fastAccessDx(i,j) /
548  sqrt(value_type(1)-expr.val()*expr.val()))
551  using std::asin; using std::sqrt;,
552  asin(expr.val()),
553  asin(expr.val(j)),
554  expr.dx(i,j)/ sqrt(value_type(1)-expr.val()*expr.val()),
555  expr.fastAccessDx(i,j) /
556  sqrt(value_type(1)-expr.val()*expr.val()))
559  using std::atan;,
560  atan(expr.val()),
561  atan(expr.val(j)),
562  expr.dx(i,j)/(value_type(1)+expr.val()*expr.val()),
563  expr.fastAccessDx(i,j)/(value_type(1)+expr.val()*expr.val()))
566  using std::cosh; using std::sinh;,
567  cosh(expr.val()),
568  cosh(expr.val(j)),
569  expr.dx(i,j)* sinh(expr.val()),
570  expr.fastAccessDx(i,j)* sinh(expr.val()))
571 FAD_UNARYOP_MACRO(sinh,
573  using std::cosh; using std::sinh;,
574  sinh(expr.val()),
575  sinh(expr.val(j)),
576  expr.dx(i,j)* cosh(expr.val()),
577  expr.fastAccessDx(i,j)* cosh(expr.val()))
580  using std::tanh; using std::cosh;,
581  tanh(expr.val()),
582  tanh(expr.val(j)),
583  expr.dx(i,j)/( cosh(expr.val())* cosh(expr.val())),
584  expr.fastAccessDx(i,j) /
585  ( cosh(expr.val())* cosh(expr.val())))
588  using std::acosh; using std::sqrt;,
589  acosh(expr.val()),
590  acosh(expr.val(j)),
591  expr.dx(i,j)/ sqrt((expr.val()-value_type(1)) *
592  (expr.val()+value_type(1))),
593  expr.fastAccessDx(i,j)/ sqrt((expr.val()-value_type(1)) *
594  (expr.val()+value_type(1))))
597  using std::asinh; using std::sqrt;,
598  asinh(expr.val()),
599  asinh(expr.val(j)),
600  expr.dx(i,j)/ sqrt(value_type(1)+expr.val()*expr.val()),
601  expr.fastAccessDx(i,j)/ sqrt(value_type(1)+
602  expr.val()*expr.val()))
605  using std::atanh;,
606  atanh(expr.val()),
607  atanh(expr.val(j)),
608  expr.dx(i,j)/(value_type(1)-expr.val()*expr.val()),
609  expr.fastAccessDx(i,j)/(value_type(1)-
610  expr.val()*expr.val()))
613  using std::abs;,
614  abs(expr.val()),
615  abs(expr.val(j)),
616  if_then_else( expr.val() >= 0, expr.dx(i,j), value_type(-expr.dx(i,j)) ),
617  if_then_else( expr.val() >= 0, expr.fastAccessDx(i,j), value_type(-expr.fastAccessDx(i,j)) ) )
620  using std::fabs;,
621  fabs(expr.val()),
622  fabs(expr.val(j)),
623  if_then_else( expr.val() >= 0, expr.dx(i,j), value_type(-expr.dx(i,j)) ),
624  if_then_else( expr.val() >= 0, expr.fastAccessDx(i,j), value_type(-expr.fastAccessDx(i,j)) ) )
627  using std::cbrt;,
628  cbrt(expr.val()),
629  cbrt(expr.val(j)),
630  expr.dx(i,j)/(value_type(3)*cbrt(expr.val()*expr.val())),
631  expr.fastAccessDx(i,j)/(value_type(3)*cbrt(expr.val()*expr.val())))
632 
633 #undef FAD_UNARYOP_MACRO
634 
635 namespace Sacado {
636  namespace Fad {
637  namespace Exp {
638 
639  // For MP::Vector scalar type, promote constants up to expression's value
640  // type. If the constant type is the same as the value type, we can store
641  // the constant as a reference. If it isn't, we must copy it into a new
642  // value type object. We do this so that we can always access the constant
643  // as a value type.
644  template <typename ConstType, typename ValueType>
645  struct ConstTypeRef {
646  typedef ValueType type;
647  };
648 
649  template <typename ValueType>
650  struct ConstTypeRef<ValueType, ValueType> {
651  typedef ValueType& type;
652  };
653  }
654  }
655 }
656 
657 #define FAD_BINARYOP_MACRO(OPNAME,OP,USING,MPVALUE,VALUE,DX,CDX1,CDX2,FASTACCESSDX,MPVAL_CONST_DX_1,MPVAL_CONST_DX_2,VAL_CONST_DX_1,VAL_CONST_DX_2,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
658 namespace Sacado { \
659  namespace Fad { \
660  namespace Exp { \
661  \
662  template <typename T1, typename T2 > \
663  class OP< T1, T2, false, false, ExprSpecMPVector > : \
664  public Expr< OP< T1, T2, false, false, ExprSpecMPVector > > { \
665  public: \
666  \
667  typedef typename std::remove_cv<T1>::type ExprT1; \
668  typedef typename std::remove_cv<T2>::type ExprT2; \
669  typedef typename ExprT1::value_type value_type_1; \
670  typedef typename ExprT2::value_type value_type_2; \
671  typedef typename Sacado::Promote<value_type_1, \
672  value_type_2>::type value_type; \
673  \
674  typedef typename ExprT1::scalar_type scalar_type_1; \
675  typedef typename ExprT2::scalar_type scalar_type_2; \
676  typedef typename Sacado::Promote<scalar_type_1, \
677  scalar_type_2>::type scalar_type; \
678  \
679  typedef typename value_type::value_type val_type; \
680  \
681  typedef ExprSpecMPVector expr_spec_type; \
682  \
683  KOKKOS_INLINE_FUNCTION \
684  OP(const T1& expr1_, const T2& expr2_) : \
685  expr1(expr1_), expr2(expr2_) {} \
686  \
687  KOKKOS_INLINE_FUNCTION \
688  int size() const { \
689  const int sz1 = expr1.size(), sz2 = expr2.size(); \
690  return sz1 > sz2 ? sz1 : sz2; \
691  } \
692  \
693  KOKKOS_INLINE_FUNCTION \
694  bool hasFastAccess() const { \
695  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
696  } \
697  \
698  KOKKOS_INLINE_FUNCTION \
699  value_type val() const { \
700  USING \
701  return MPVALUE; \
702  } \
703  \
704  KOKKOS_INLINE_FUNCTION \
705  val_type val(int j) const { \
706  USING \
707  return VALUE; \
708  } \
709  \
710  KOKKOS_INLINE_FUNCTION \
711  val_type dx(int i, int j) const { \
712  USING \
713  const int sz1 = expr1.size(), sz2 = expr2.size(); \
714  if (sz1 > 0 && sz2 > 0) \
715  return DX; \
716  else if (sz1 > 0) \
717  return CDX2; \
718  else \
719  return CDX1; \
720  } \
721  \
722  KOKKOS_INLINE_FUNCTION \
723  val_type fastAccessDx(int i, int j) const { \
724  USING \
725  return FASTACCESSDX; \
726  } \
727  \
728  protected: \
729  \
730  const T1& expr1; \
731  const T2& expr2; \
732  \
733  }; \
734  \
735  template <typename T1, typename T2> \
736  class OP< T1, T2, false, true, ExprSpecMPVector > : \
737  public Expr< OP< T1, T2, false, true, ExprSpecMPVector > > { \
738  public: \
739  \
740  typedef typename std::remove_cv<T1>::type ExprT1; \
741  typedef T2 ConstT; \
742  typedef typename ExprT1::value_type value_type; \
743  typedef typename ExprT1::scalar_type scalar_type; \
744  \
745  typedef typename value_type::value_type val_type; \
746  \
747  typedef ExprSpecMPVector expr_spec_type; \
748  \
749  KOKKOS_INLINE_FUNCTION \
750  OP(const T1& expr1_, const ConstT& c_) : \
751  expr1(expr1_), c(c_) {} \
752  \
753  KOKKOS_INLINE_FUNCTION \
754  int size() const { \
755  return expr1.size(); \
756  } \
757  \
758  KOKKOS_INLINE_FUNCTION \
759  bool hasFastAccess() const { \
760  return expr1.hasFastAccess(); \
761  } \
762  \
763  KOKKOS_INLINE_FUNCTION \
764  value_type val() const { \
765  USING \
766  return MPVAL_CONST_DX_2; \
767  } \
768  \
769  KOKKOS_INLINE_FUNCTION \
770  val_type val(int j) const { \
771  USING \
772  return VAL_CONST_DX_2; \
773  } \
774  \
775  KOKKOS_INLINE_FUNCTION \
776  val_type dx(int i, int j) const { \
777  USING \
778  return CONST_DX_2; \
779  } \
780  \
781  KOKKOS_INLINE_FUNCTION \
782  val_type fastAccessDx(int i, int j) const { \
783  USING \
784  return CONST_FASTACCESSDX_2; \
785  } \
786  \
787  protected: \
788  \
789  const T1& expr1; \
790  const typename ConstTypeRef<ConstT,value_type>::type c; \
791  }; \
792  \
793  template <typename T1, typename T2> \
794  class OP< T1, T2, true, false,ExprSpecMPVector > : \
795  public Expr< OP< T1, T2, true, false, ExprSpecMPVector > > { \
796  public: \
797  \
798  typedef typename std::remove_cv<T2>::type ExprT2; \
799  typedef T1 ConstT; \
800  typedef typename ExprT2::value_type value_type; \
801  typedef typename ExprT2::scalar_type scalar_type; \
802  \
803  typedef typename value_type::value_type val_type; \
804  \
805  typedef ExprSpecMPVector expr_spec_type; \
806  \
807  KOKKOS_INLINE_FUNCTION \
808  OP(const ConstT& c_, const T2& expr2_) : \
809  c(c_), expr2(expr2_) {} \
810  \
811  KOKKOS_INLINE_FUNCTION \
812  int size() const { \
813  return expr2.size(); \
814  } \
815  \
816  KOKKOS_INLINE_FUNCTION \
817  bool hasFastAccess() const { \
818  return expr2.hasFastAccess(); \
819  } \
820  \
821  KOKKOS_INLINE_FUNCTION \
822  value_type val() const { \
823  USING \
824  return MPVAL_CONST_DX_1; \
825  } \
826  \
827  KOKKOS_INLINE_FUNCTION \
828  val_type val(int j) const { \
829  USING \
830  return VAL_CONST_DX_1; \
831  } \
832  \
833  KOKKOS_INLINE_FUNCTION \
834  val_type dx(int i, int j) const { \
835  USING \
836  return CONST_DX_1; \
837  } \
838  \
839  KOKKOS_INLINE_FUNCTION \
840  val_type fastAccessDx(int i, int j) const { \
841  USING \
842  return CONST_FASTACCESSDX_1; \
843  } \
844  \
845  protected: \
846  \
847  const typename ConstTypeRef<ConstT,value_type>::type c; \
848  const T2& expr2; \
849  }; \
850  \
851  } \
852  } \
853  \
854 }
855 
856 
857 FAD_BINARYOP_MACRO(operator+,
858  AdditionOp,
859  ;,
860  expr1.val() + expr2.val(),
861  expr1.val(j) + expr2.val(j),
862  expr1.dx(i,j) + expr2.dx(i,j),
863  expr2.dx(i,j),
864  expr1.dx(i,j),
865  expr1.fastAccessDx(i,j) + expr2.fastAccessDx(i,j),
866  c + expr2.val(),
867  expr1.val() + c,
868  c.fastAccessCoeff(j) + expr2.val(j),
869  expr1.val(j) + c.fastAccessCoeff(j),
870  expr2.dx(i,j),
871  expr1.dx(i,j),
872  expr2.fastAccessDx(i,j),
873  expr1.fastAccessDx(i,j))
874 FAD_BINARYOP_MACRO(operator-,
876  ;,
877  expr1.val() - expr2.val(),
878  expr1.val(j) - expr2.val(j),
879  expr1.dx(i,j) - expr2.dx(i,j),
880  -expr2.dx(i,j),
881  expr1.dx(i,j),
882  expr1.fastAccessDx(i,j) - expr2.fastAccessDx(i,j),
883  c - expr2.val(),
884  expr1.val() - c,
885  c.fastAccessCoeff(j) - expr2.val(j),
886  expr1.val(j) - c.fastAccessCoeff(j),
887  -expr2.dx(i,j),
888  expr1.dx(i,j),
889  -expr2.fastAccessDx(i,j),
890  expr1.fastAccessDx(i,j))
891 FAD_BINARYOP_MACRO(operator*,
893  ;,
894  expr1.val() * expr2.val(),
895  expr1.val(j) * expr2.val(j),
896  expr1.val(j)*expr2.dx(i,j) + expr1.dx(i,j)*expr2.val(j),
897  expr1.val(j)*expr2.dx(i,j),
898  expr1.dx(i,j)*expr2.val(j),
899  expr1.val(j)*expr2.fastAccessDx(i,j) +
900  expr1.fastAccessDx(i,j)*expr2.val(j),
901  c * expr2.val(),
902  expr1.val() * c,
903  c.fastAccessCoeff(j) * expr2.val(j),
904  expr1.val(j) * c.fastAccessCoeff(j),
905  c.fastAccessCoeff(j)*expr2.dx(i,j),
906  expr1.dx(i,j)*c.fastAccessCoeff(j),
907  c.fastAccessCoeff(j)*expr2.fastAccessDx(i,j),
908  expr1.fastAccessDx(i,j)*c.fastAccessCoeff(j))
909 FAD_BINARYOP_MACRO(operator/,
911  ;,
912  expr1.val() / expr2.val(),
913  expr1.val(j) / expr2.val(j),
914  (expr1.dx(i,j)*expr2.val(j) - expr2.dx(i,j)*expr1.val(j)) /
915  (expr2.val(j)*expr2.val(j)),
916  -expr2.dx(i,j)*expr1.val(j) / (expr2.val(j)*expr2.val(j)),
917  expr1.dx(i,j)/expr2.val(j),
918  (expr1.fastAccessDx(i,j)*expr2.val(j) -
919  expr2.fastAccessDx(i,j)*expr1.val(j)) /
920  (expr2.val(j)*expr2.val(j)),
921  c / expr2.val(),
922  expr1.val() / c,
923  c.fastAccessCoeff(j) / expr2.val(j),
924  expr1.val(j) / c.fastAccessCoeff(j),
925  -expr2.dx(i,j)*c.fastAccessCoeff(j) / (expr2.val(j)*expr2.val(j)),
926  expr1.dx(i,j)/c.fastAccessCoeff(j),
927  -expr2.fastAccessDx(i,j)*c.fastAccessCoeff(j) / (expr2.val(j)*expr2.val(j)),
928  expr1.fastAccessDx(i,j)/c.fastAccessCoeff(j))
931  using std::atan2;,
932  atan2(expr1.val(), expr2.val()),
933  atan2(expr1.val(j), expr2.val(j)),
934  (expr2.val(j)*expr1.dx(i,j) - expr1.val(j)*expr2.dx(i,j))/
935  (expr1.val(j)*expr1.val(j) + expr2.val(j)*expr2.val(j)),
936  -expr1.val(j)*expr2.dx(i,j)/
937  (expr1.val(j)*expr1.val(j) + expr2.val(j)*expr2.val(j)),
938  expr2.val(j)*expr1.dx(i,j)/
939  (expr1.val(j)*expr1.val(j) + expr2.val(j)*expr2.val(j)),
940  (expr2.val(j)*expr1.fastAccessDx(i,j) - expr1.val(j)*expr2.fastAccessDx(i,j))/
941  (expr1.val(j)*expr1.val(j) + expr2.val(j)*expr2.val(j)),
942  atan2(c, expr2.val()),
943  atan2(expr1.val(), c),
944  atan2(c.fastAccessCoeff(j), expr2.val(j)),
945  atan2(expr1.val(j), c.fastAccessCoeff(j)),
946  (-c.fastAccessCoeff(j)*expr2.dx(i,j)) / (c.fastAccessCoeff(j)*c.fastAccessCoeff(j) + expr2.val(j)*expr2.val(j)),
947  (c.fastAccessCoeff(j)*expr1.dx(i,j))/ (expr1.val(j)*expr1.val(j) + c.fastAccessCoeff(j)*c.fastAccessCoeff(j)),
948  (-c.fastAccessCoeff(j)*expr2.fastAccessDx(i,j))/ (c.fastAccessCoeff(j)*c.fastAccessCoeff(j) + expr2.val(j)*expr2.val(j)),
949  (c.fastAccessCoeff(j)*expr1.fastAccessDx(i,j))/ (expr1.val(j)*expr1.val(j) + c.fastAccessCoeff(j)*c.fastAccessCoeff(j)))
950 // FAD_BINARYOP_MACRO(pow,
951 // PowerOp,
952 // using std::pow; using std::log;,
953 // pow(expr1.val(), expr2.val()),
954 // pow(expr1.val(j), expr2.val(j)),
955 // if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type((expr2.dx(i,j)*log(expr1.val(j))+expr2.val(j)*expr1.dx(i,j)/expr1.val(j))*pow(expr1.val(j),expr2.val(j))) ),
956 // if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(expr2.dx(i,j)*log(expr1.val(j))*pow(expr1.val(j),expr2.val(j))) ),
957 // if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(expr2.val(j)*expr1.dx(i,j)/expr1.val(j)*pow(expr1.val(j),expr2.val(j))) ),
958 // if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type((expr2.fastAccessDx(i,j)*log(expr1.val(j))+expr2.val(j)*expr1.fastAccessDx(i,j)/expr1.val(j))*pow(expr1.val(j),expr2.val(j))) ),
959 // pow(c, expr2.val()),
960 // pow(expr1.val(), c),
961 // pow(c.fastAccessCoeff(j), expr2.val(j)),
962 // pow(expr1.val(j), c.fastAccessCoeff(j)),
963 // if_then_else( c.fastAccessCoeff(j) == val_type(0.0), val_type(0.0), val_type(expr2.dx(i,j)*log(c.fastAccessCoeff(j))*pow(c.fastAccessCoeff(j),expr2.val(j))) ),
964 // if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(c.fastAccessCoeff(j)*expr1.dx(i,j)/expr1.val(j)*pow(expr1.val(j),c.fastAccessCoeff(j))) ),
965 // if_then_else( c.fastAccessCoeff(j) == val_type(0.0), val_type(0.0), val_type(expr2.fastAccessDx(i,j)*log(c.fastAccessCoeff(j))*pow(c.fastAccessCoeff(j),expr2.val(j))) ),
966 // if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(c.fastAccessCoeff(j)*expr1.fastAccessDx(i,j)/expr1.val(j)*pow(expr1.val(j),c.fastAccessCoeff(j)))) )
969  ;,
970  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
971  if_then_else( expr1.val(j) >= expr2.val(j), expr1.val(j), expr2.val(j) ),
972  if_then_else( expr1.val(j) >= expr2.val(j), expr1.dx(i,j), expr2.dx(i,j) ),
973  if_then_else( expr1.val(j) >= expr2.val(j), val_type(0.0), expr2.dx(i,j) ),
974  if_then_else( expr1.val(j) >= expr2.val(j), expr1.dx(i,j), val_type(0.0) ),
975  if_then_else( expr1.val(j) >= expr2.val(j), expr1.fastAccessDx(i,j), expr2.fastAccessDx(i,j) ),
976  if_then_else( c >= expr2.val(), c, expr2.val() ),
977  if_then_else( expr1.val() >= c, expr1.val(), c ),
978  if_then_else( c.fastAccessCoeff(j) >= expr2.val(j), c.fastAccessCoeff(j), expr2.val(j) ),
979  if_then_else( expr1.val(j) >= c.fastAccessCoeff(j), expr1.val(j), c.fastAccessCoeff(j) ),
980  if_then_else( c.fastAccessCoeff(j) >= expr2.val(j), val_type(0.0), expr2.dx(i,j) ),
981  if_then_else( expr1.val(j) >= c.fastAccessCoeff(j), expr1.dx(i,j), val_type(0.0) ),
982  if_then_else( c.fastAccessCoeff(j) >= expr2.val(j), val_type(0.0), expr2.fastAccessDx(i,j) ),
983  if_then_else( expr1.val(j) >= c.fastAccessCoeff(j), expr1.fastAccessDx(i,j), val_type(0.0) ) )
986  ;,
987  if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
988  if_then_else( expr1.val(j) <= expr2.val(j), expr1.val(j), expr2.val(j) ),
989  if_then_else( expr1.val(j) <= expr2.val(j), expr1.dx(i,j), expr2.dx(i,j) ),
990  if_then_else( expr1.val(j) <= expr2.val(j), val_type(0.0), expr2.dx(i,j) ),
991  if_then_else( expr1.val(j) <= expr2.val(j), expr1.dx(i,j), val_type(0.0) ),
992  if_then_else( expr1.val(j) <= expr2.val(j), expr1.fastAccessDx(i,j), expr2.fastAccessDx(i,j) ),
993  if_then_else( c <= expr2.val(), c, expr2.val() ),
994  if_then_else( expr1.val() <= c, expr1.val(), c ),
995  if_then_else( c.fastAccessCoeff(j) <= expr2.val(j), c.fastAccessCoeff(j), expr2.val(j) ),
996  if_then_else( expr1.val(j) <= c.fastAccessCoeff(j), expr1.val(j), c.fastAccessCoeff(j) ),
997  if_then_else( c.fastAccessCoeff(j) <= expr2.val(j), val_type(0), expr2.dx(i,j) ),
998  if_then_else( expr1.val(j) <= c.fastAccessCoeff(j), expr1.dx(i,j), val_type(0) ),
999  if_then_else( c.fastAccessCoeff(j) <= expr2.val(j), val_type(0), expr2.fastAccessDx(i,j) ),
1000  if_then_else( expr1.val(j) <= c.fastAccessCoeff(j), expr1.fastAccessDx(i,j), val_type(0) ) )
1001 
1002 // Special handling for std::pow() to provide specializations of PowerOp for
1003 // "simd" value types that use if_then_else(). The only reason for not using
1004 // if_then_else() always is to avoid evaluating the derivative if the value is
1005 // zero to avoid throwing FPEs.
1006 namespace Sacado {
1007  namespace Fad {
1008  namespace Exp {
1009 
1010  //
1011  // Implementation for simd type using if_then_else()
1012  //
1013  template <typename T1, typename T2>
1014  class PowerOp< T1, T2, false, false, ExprSpecMPVector, PowerImpl::Simd > :
1015  public Expr< PowerOp< T1, T2, false, false, ExprSpecMPVector,
1016  PowerImpl::Simd > > {
1017  public:
1018 
1019  typedef typename std::remove_cv<T1>::type ExprT1;
1020  typedef typename std::remove_cv<T2>::type ExprT2;
1021  typedef typename ExprT1::value_type value_type_1;
1022  typedef typename ExprT2::value_type value_type_2;
1023  typedef typename Sacado::Promote<value_type_1,
1024  value_type_2>::type value_type;
1025 
1026  typedef typename ExprT1::scalar_type scalar_type_1;
1027  typedef typename ExprT2::scalar_type scalar_type_2;
1028  typedef typename Sacado::Promote<scalar_type_1,
1029  scalar_type_2>::type scalar_type;
1030 
1031  typedef typename value_type::value_type val_type;
1032 
1033  typedef ExprSpecMPVector expr_spec_type;
1034 
1035  KOKKOS_INLINE_FUNCTION
1036  PowerOp(const T1& expr1_, const T2& expr2_) :
1037  expr1(expr1_), expr2(expr2_) {}
1038 
1039  KOKKOS_INLINE_FUNCTION
1040  int size() const {
1041  const int sz1 = expr1.size(), sz2 = expr2.size();
1042  return sz1 > sz2 ? sz1 : sz2;
1043  }
1044 
1045  KOKKOS_INLINE_FUNCTION
1046  bool hasFastAccess() const {
1047  return expr1.hasFastAccess() && expr2.hasFastAccess();
1048  }
1049 
1050  KOKKOS_INLINE_FUNCTION
1051  value_type val() const {
1052  using std::pow;
1053  return pow(expr1.val(), expr2.val());
1054  }
1055 
1056  KOKKOS_INLINE_FUNCTION
1057  val_type val(int j) const {
1058  using std::pow;
1059  return pow(expr1.val(j), expr2.val(j));
1060  }
1061 
1062  KOKKOS_INLINE_FUNCTION
1063  val_type dx(int i, int j) const {
1064  using std::pow; using std::log;
1065  const int sz1 = expr1.size(), sz2 = expr2.size();
1066  if (sz1 > 0 && sz2 > 0)
1067  return if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type((expr2.dx(i,j)*log(expr1.val(j))+expr2.val(j)*expr1.dx(i,j)/expr1.val(j))*pow(expr1.val(j),expr2.val(j))) );
1068  else if (sz1 > 0)
1069  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1070  // It seems less accurate and caused convergence problems in some codes
1071  return if_then_else( expr2.val(j) == scalar_type(1.0), expr1.dx(i,j), if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(expr2.val(j)*expr1.dx(i,j)/expr1.val(j)*pow(expr1.val(j),expr2.val(j))) ));
1072  else
1073  return if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(expr2.dx(i,j)*log(expr1.val(j))*pow(expr1.val(j),expr2.val(j))) );
1074  }
1075 
1076  KOKKOS_INLINE_FUNCTION
1077  val_type fastAccessDx(int i, int j) const {
1078  using std::pow; using std::log;
1079  return if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type((expr2.fastAccessDx(i,j)*log(expr1.val(j))+expr2.val(j)*expr1.fastAccessDx(i,j)/expr1.val(j))*pow(expr1.val(j),expr2.val(j))) );
1080  }
1081 
1082  protected:
1083 
1084  const T1& expr1;
1085  const T2& expr2;
1086 
1087  };
1088 
1089  template <typename T1, typename T2>
1090  class PowerOp< T1, T2, false, true, ExprSpecMPVector, PowerImpl::Simd >
1091  : public Expr< PowerOp< T1, T2, false, true, ExprSpecMPVector,
1092  PowerImpl::Simd > > {
1093  public:
1094 
1095  typedef typename std::remove_cv<T1>::type ExprT1;
1096  typedef T2 ConstT;
1097  typedef typename ExprT1::value_type value_type;
1098  typedef typename ExprT1::scalar_type scalar_type;
1099 
1100  typedef typename value_type::value_type val_type;
1101 
1102  typedef ExprSpecMPVector expr_spec_type;
1103 
1104  KOKKOS_INLINE_FUNCTION
1105  PowerOp(const T1& expr1_, const ConstT& c_) :
1106  expr1(expr1_), c(c_) {}
1107 
1108  KOKKOS_INLINE_FUNCTION
1109  int size() const {
1110  return expr1.size();
1111  }
1112 
1113  KOKKOS_INLINE_FUNCTION
1114  bool hasFastAccess() const {
1115  return expr1.hasFastAccess();
1116  }
1117 
1118  KOKKOS_INLINE_FUNCTION
1119  value_type val() const {
1120  using std::pow;
1121  return pow(expr1.val(), c);
1122  }
1123 
1124  KOKKOS_INLINE_FUNCTION
1125  val_type val(int j) const {
1126  using std::pow;
1127  return pow(expr1.val(j), c.fastAccessCoeff(j));
1128  }
1129 
1130  KOKKOS_INLINE_FUNCTION
1131  val_type dx(int i, int j) const {
1132  using std::pow;
1133  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1134  // It seems less accurate and caused convergence problems in some codes
1135  return if_then_else( c.fastAccessCoeff(j) == scalar_type(1.0), expr1.dx(i,j), if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(c.fastAccessCoeff(j)*expr1.dx(i,j)/expr1.val(j)*pow(expr1.val(j),c.fastAccessCoeff(j))) ));
1136  }
1137 
1138  KOKKOS_INLINE_FUNCTION
1139  val_type fastAccessDx(int i, int j) const {
1140  using std::pow;
1141  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1142  // It seems less accurate and caused convergence problems in some codes
1143  return if_then_else( c.fastAccessCoeff(j) == scalar_type(1.0), expr1.fastAccessDx(i,j), if_then_else( expr1.val(j) == val_type(0.0), val_type(0.0), val_type(c.fastAccessCoeff(j)*expr1.fastAccessDx(i,j)/expr1.val(j)*pow(expr1.val(j),c.fastAccessCoeff(j))) ));
1144  }
1145 
1146  protected:
1147 
1148  const T1& expr1;
1149  const ConstT& c;
1150  };
1151 
1152  template <typename T1, typename T2>
1153  class PowerOp< T1, T2, true, false, ExprSpecMPVector, PowerImpl::Simd >
1154  : public Expr< PowerOp< T1, T2, true, false, ExprSpecMPVector,
1155  PowerImpl::Simd> > {
1156  public:
1157 
1158  typedef typename std::remove_cv<T2>::type ExprT2;
1159  typedef T1 ConstT;
1160  typedef typename ExprT2::value_type value_type;
1161  typedef typename ExprT2::scalar_type scalar_type;
1162 
1163  typedef typename value_type::value_type val_type;
1164 
1165  typedef ExprSpecMPVector expr_spec_type;
1166 
1167  KOKKOS_INLINE_FUNCTION
1168  PowerOp(const ConstT& c_, const T2& expr2_) :
1169  c(c_), expr2(expr2_) {}
1170 
1171  KOKKOS_INLINE_FUNCTION
1172  int size() const {
1173  return expr2.size();
1174  }
1175 
1176  KOKKOS_INLINE_FUNCTION
1177  bool hasFastAccess() const {
1178  return expr2.hasFastAccess();
1179  }
1180 
1181  KOKKOS_INLINE_FUNCTION
1182  value_type val() const {
1183  using std::pow;
1184  return pow(c, expr2.val());
1185  }
1186 
1187  KOKKOS_INLINE_FUNCTION
1188  val_type val(int j) const {
1189  using std::pow;
1190  return pow(c.fastAccessCoeff(j), expr2.val(j));
1191  }
1192 
1193  KOKKOS_INLINE_FUNCTION
1194  val_type dx(int i, int j) const {
1195  using std::pow; using std::log;
1196  return if_then_else( c.fastAccessCoeff(j) == val_type(0.0), val_type(0.0), val_type(expr2.dx(i,j)*log(c.fastAccessCoeff(j))*pow(c.fastAccessCoeff(j),expr2.val(j))) );
1197  }
1198 
1199  KOKKOS_INLINE_FUNCTION
1200  val_type fastAccessDx(int i, int j) const {
1201  using std::pow; using std::log;
1202  return if_then_else( c.fastAccessCoeff(j) == val_type(0.0), val_type(0.0), val_type(expr2.fastAccessDx(i,j)*log(c.fastAccessCoeff(j))*pow(c.fastAccessCoeff(j),expr2.val(j))) );
1203  }
1204 
1205  protected:
1206 
1207  const ConstT& c;
1208  const T2& expr2;
1209  };
1210 
1211  //
1212  // Specialization for nested derivatives. This version does not use
1213  // if_then_else/ternary-operator on the base so that nested derivatives
1214  // are correct.
1215  //
1216  template <typename T1, typename T2>
1217  class PowerOp< T1, T2, false, false, ExprSpecMPVector,
1218  PowerImpl::NestedSimd > :
1219  public Expr< PowerOp< T1, T2, false, false, ExprSpecMPVector,
1220  PowerImpl::NestedSimd > > {
1221  public:
1222 
1223  typedef typename std::remove_cv<T1>::type ExprT1;
1224  typedef typename std::remove_cv<T2>::type ExprT2;
1225  typedef typename ExprT1::value_type value_type_1;
1226  typedef typename ExprT2::value_type value_type_2;
1227  typedef typename Sacado::Promote<value_type_1,
1228  value_type_2>::type value_type;
1229 
1230  typedef typename ExprT1::scalar_type scalar_type_1;
1231  typedef typename ExprT2::scalar_type scalar_type_2;
1232  typedef typename Sacado::Promote<scalar_type_1,
1233  scalar_type_2>::type scalar_type;
1234 
1235  typedef typename value_type::value_type val_type;
1236 
1237  typedef ExprSpecMPVector expr_spec_type;
1238 
1239  KOKKOS_INLINE_FUNCTION
1240  PowerOp(const T1& expr1_, const T2& expr2_) :
1241  expr1(expr1_), expr2(expr2_) {}
1242 
1243  KOKKOS_INLINE_FUNCTION
1244  int size() const {
1245  const int sz1 = expr1.size(), sz2 = expr2.size();
1246  return sz1 > sz2 ? sz1 : sz2;
1247  }
1248 
1249  KOKKOS_INLINE_FUNCTION
1250  bool hasFastAccess() const {
1251  return expr1.hasFastAccess() && expr2.hasFastAccess();
1252  }
1253 
1254  KOKKOS_INLINE_FUNCTION
1255  value_type val() const {
1256  using std::pow;
1257  return pow(expr1.val(), expr2.val());
1258  }
1259 
1260  KOKKOS_INLINE_FUNCTION
1261  val_type val(int j) const {
1262  using std::pow;
1263  return pow(expr1.val(j), expr2.val(j));
1264  }
1265 
1266  KOKKOS_INLINE_FUNCTION
1267  value_type dx(int i, int j) const {
1268  using std::pow; using std::log;
1269  const int sz1 = expr1.size(), sz2 = expr2.size();
1270  if (sz1 > 0 && sz2 > 0)
1271  return (expr2.dx(i,j)*log(expr1.val(j))+expr2.val(j)*expr1.dx(i,j)/expr1.val(j))*pow(expr1.val(j),expr2.val(j));
1272  else if (sz1 > 0)
1273  return if_then_else( expr2.val(j) == scalar_type(0.0), value_type(0.0), value_type((expr2.val(j)*expr1.dx(i,j))*pow(expr1.val(j),expr2.val(j)-scalar_type(1.0))));
1274  else
1275  return expr2.dx(i,j)*log(expr1.val(j))*pow(expr1.val(j),expr2.val(j));
1276  }
1277 
1278  KOKKOS_INLINE_FUNCTION
1279  value_type fastAccessDx(int i, int j) const {
1280  using std::pow; using std::log;
1281  return (expr2.fastAccessDx(i,j)*log(expr1.val(j))+expr2.val(j)*expr1.fastAccessDx(i,j)/expr1.val(j))*pow(expr1.val(j),expr2.val(j));
1282  }
1283 
1284  protected:
1285 
1286  const T1& expr1;
1287  const T2& expr2;
1288 
1289  };
1290 
1291  template <typename T1, typename T2>
1292  class PowerOp< T1, T2, false, true, ExprSpecMPVector,
1293  PowerImpl::NestedSimd > :
1294  public Expr< PowerOp< T1, T2, false, true, ExprSpecMPVector,
1295  PowerImpl::NestedSimd > > {
1296  public:
1297 
1298  typedef typename std::remove_cv<T1>::type ExprT1;
1299  typedef T2 ConstT;
1300  typedef typename ExprT1::value_type value_type;
1301  typedef typename ExprT1::scalar_type scalar_type;
1302 
1303  typedef typename value_type::value_type val_type;
1304 
1305  typedef ExprSpecMPVector expr_spec_type;
1306 
1307  KOKKOS_INLINE_FUNCTION
1308  PowerOp(const T1& expr1_, const ConstT& c_) :
1309  expr1(expr1_), c(c_) {}
1310 
1311  KOKKOS_INLINE_FUNCTION
1312  int size() const {
1313  return expr1.size();
1314  }
1315 
1316  KOKKOS_INLINE_FUNCTION
1317  bool hasFastAccess() const {
1318  return expr1.hasFastAccess();
1319  }
1320 
1321  KOKKOS_INLINE_FUNCTION
1322  value_type val() const {
1323  using std::pow;
1324  return pow(expr1.val(), c);
1325  }
1326 
1327  KOKKOS_INLINE_FUNCTION
1328  val_type val(int j) const {
1329  using std::pow;
1330  return pow(expr1.val(j), c.fastAccessCoeff(j));
1331  }
1332 
1333  KOKKOS_INLINE_FUNCTION
1334  value_type dx(int i, int j) const {
1335  using std::pow;
1336  return if_then_else( c.fastAccessCoeff(j) == scalar_type(0.0), value_type(0.0), value_type(c.fastAccessCoeff(j)*expr1.dx(i,j)*pow(expr1.val(j),c.fastAccessCoeff(j)-scalar_type(1.0))));
1337  }
1338 
1339  KOKKOS_INLINE_FUNCTION
1340  value_type fastAccessDx(int i, int j) const {
1341  using std::pow;
1342  return if_then_else( c.fastAccessCoeff(j) == scalar_type(0.0), value_type(0.0), value_type(c.fastAccessCoeff(j)*expr1.fastAccessDx(i,j)*pow(expr1.val(j),c.fastAccessCoeff(j)-scalar_type(1.0))));
1343  }
1344 
1345  protected:
1346 
1347  const T1& expr1;
1348  const ConstT& c;
1349  };
1350 
1351  template <typename T1, typename T2>
1352  class PowerOp<T1, T2, true, false, ExprSpecMPVector,
1353  PowerImpl::NestedSimd > :
1354  public Expr< PowerOp< T1, T2, true, false, ExprSpecMPVector,
1355  PowerImpl::NestedSimd > > {
1356  public:
1357 
1358  typedef typename std::remove_cv<T2>::type ExprT2;
1359  typedef T1 ConstT;
1360  typedef typename ExprT2::value_type value_type;
1361  typedef typename ExprT2::scalar_type scalar_type;
1362 
1363  typedef typename value_type::value_type val_type;
1364 
1365  typedef ExprSpecMPVector expr_spec_type;
1366 
1367  KOKKOS_INLINE_FUNCTION
1368  PowerOp(const ConstT& c_, const T2& expr2_) :
1369  c(c_), expr2(expr2_) {}
1370 
1371  KOKKOS_INLINE_FUNCTION
1372  int size() const {
1373  return expr2.size();
1374  }
1375 
1376  KOKKOS_INLINE_FUNCTION
1377  bool hasFastAccess() const {
1378  return expr2.hasFastAccess();
1379  }
1380 
1381  KOKKOS_INLINE_FUNCTION
1382  value_type val() const {
1383  using std::pow;
1384  return pow(c, expr2.val());
1385  }
1386 
1387  KOKKOS_INLINE_FUNCTION
1388  val_type val(int j) const {
1389  using std::pow;
1390  return pow(c.fastAccessCoeff(j), expr2.val(j));
1391  }
1392 
1393  KOKKOS_INLINE_FUNCTION
1394  value_type dx(int i, int j) const {
1395  using std::pow; using std::log;
1396  return expr2.dx(i,j)*log(c.fastAccessCoeff(j))*pow(c.fastAccessCoeff(j),expr2.val(j));
1397  }
1398 
1399  KOKKOS_INLINE_FUNCTION
1400  value_type fastAccessDx(int i, int j) const {
1401  using std::pow; using std::log;
1402  return expr2.fastAccessDx(i,j)*log(c.fastAccessCoeff(j))*pow(c.fastAccessCoeff(j),expr2.val(j));
1403  }
1404 
1405  protected:
1406 
1407  const ConstT& c;
1408  const T2& expr2;
1409  };
1410 
1411  }
1412  }
1413 }
1414 
1415 //--------------------------if_then_else operator -----------------------
1416 // Can't use the above macros because it is a ternary operator (sort of).
1417 // Also, relies on C++11
1418 
1419 namespace Sacado {
1420  namespace Fad {
1421  namespace Exp {
1422 
1423  template <typename CondT, typename T1, typename T2>
1424  class IfThenElseOp< CondT,T1,T2,false,false,ExprSpecMPVector > :
1425  public Expr< IfThenElseOp< CondT, T1, T2, false, false, ExprSpecMPVector > > {
1426 
1427  public:
1428 
1429  typedef typename std::remove_cv<T1>::type ExprT1;
1430  typedef typename std::remove_cv<T2>::type ExprT2;
1433  typedef typename Sacado::Promote<value_type_1,
1435 
1438  typedef typename Sacado::Promote<scalar_type_1,
1440 
1442 
1444 
1445  KOKKOS_INLINE_FUNCTION
1446  IfThenElseOp(const CondT& cond_, const T1& expr1_, const T2& expr2_) :
1447  cond(cond_), expr1(expr1_), expr2(expr2_) {}
1448 
1449  KOKKOS_INLINE_FUNCTION
1450  int size() const {
1451  int sz1 = expr1.size(), sz2 = expr2.size();
1452  return sz1 > sz2 ? sz1 : sz2;
1453  }
1454 
1455  KOKKOS_INLINE_FUNCTION
1456  bool hasFastAccess() const {
1457  return expr1.hasFastAccess() && expr2.hasFastAccess();
1458  }
1459 
1460  KOKKOS_INLINE_FUNCTION
1461  value_type val() const {
1462  return if_then_else( cond, expr1.val(), expr2.val() );
1463  }
1464 
1465  KOKKOS_INLINE_FUNCTION
1466  val_type val(int j) const {
1467  return if_then_else( cond, expr1.val(j), expr2.val(j) );
1468  }
1469 
1470  KOKKOS_INLINE_FUNCTION
1471  val_type dx(int i, int j) const {
1472  return if_then_else( cond, expr1.dx(i,j), expr2.dx(i,j) );
1473  }
1474 
1475  KOKKOS_INLINE_FUNCTION
1476  val_type fastAccessDx(int i, int j) const {
1477  return if_then_else( cond, expr1.fastAccessDx(i,j), expr2.fastAccessDx(i,j) );
1478  }
1479 
1480  protected:
1481 
1482  const CondT& cond;
1483  const T1& expr1;
1484  const T2& expr2;
1485 
1486  };
1487 
1488  template <typename CondT, typename T1, typename T2>
1489  class IfThenElseOp< CondT, T1, T2, false, true, ExprSpecMPVector> :
1490  public Expr< IfThenElseOp< CondT, T1, T2, false, true, ExprSpecMPVector > > {
1491 
1492  public:
1493 
1494  typedef typename std::remove_cv<T1>::type ExprT1;
1495  typedef T2 ConstT;
1496  typedef typename ExprT1::value_type value_type;
1498 
1500 
1501  KOKKOS_INLINE_FUNCTION
1502  IfThenElseOp(const CondT& cond_, const T1& expr1_, const ConstT& c_) :
1503  cond(cond_), expr1(expr1_), c(c_) {}
1504 
1505  KOKKOS_INLINE_FUNCTION
1506  int size() const {
1507  return expr1.size();
1508  }
1509 
1510  KOKKOS_INLINE_FUNCTION
1511  bool hasFastAccess() const {
1512  return expr1.hasFastAccess();
1513  }
1514 
1515  KOKKOS_INLINE_FUNCTION
1516  value_type val() const {
1517  return if_then_else( cond, expr1.val(), c );
1518  }
1519 
1520  KOKKOS_INLINE_FUNCTION
1521  val_type val(int j) const {
1522  return if_then_else( cond, expr1.val(j), c.fastAccessCoeff(j) );
1523  }
1524 
1525  KOKKOS_INLINE_FUNCTION
1526  val_type dx(int i, int j) const {
1527  return if_then_else( cond, expr1.dx(i,j), val_type(0.0) );
1528  }
1529 
1530  KOKKOS_INLINE_FUNCTION
1531  val_type fastAccessDx(int i, int j) const {
1532  return if_then_else( cond, expr1.fastAccessDx(i,j), val_type(0.0) );
1533  }
1534 
1535  protected:
1536 
1537  const CondT& cond;
1538  const T1& expr1;
1539  const typename ConstTypeRef<ConstT,value_type>::type c;
1540  };
1541 
1542  template <typename CondT, typename T1, typename T2>
1543  class IfThenElseOp< CondT, T1, T2, true, false, ExprSpecMPVector> :
1544  public Expr< IfThenElseOp< CondT, T1, T2, true, false, ExprSpecMPVector > > {
1545 
1546  public:
1547 
1548  typedef typename std::remove_cv<T2>::type ExprT2;
1549  typedef T1 ConstT;
1550  typedef typename ExprT2::value_type value_type;
1552 
1554 
1556 
1557  KOKKOS_INLINE_FUNCTION
1558  IfThenElseOp(const CondT& cond_, const ConstT& c_, const T2& expr2_) :
1559  cond(cond_), c(c_), expr2(expr2_) {}
1560 
1561  KOKKOS_INLINE_FUNCTION
1562  int size() const {
1563  return expr2.size();
1564  }
1565 
1566  KOKKOS_INLINE_FUNCTION
1567  bool hasFastAccess() const {
1568  return expr2.hasFastAccess();
1569  }
1570 
1571  KOKKOS_INLINE_FUNCTION
1572  value_type val() const {
1573  return if_then_else( cond, c, expr2.val() );
1574  }
1575 
1576  KOKKOS_INLINE_FUNCTION
1577  val_type val(int j) const {
1578  return if_then_else( cond, c.fastAccessCoeff(j), expr2.val(j) );
1579  }
1580 
1581  KOKKOS_INLINE_FUNCTION
1582  val_type dx(int i, int j) const {
1583  return if_then_else( cond, val_type(0.0), expr2.dx(i,j) );
1584  }
1585 
1586  KOKKOS_INLINE_FUNCTION
1587  val_type fastAccessDx(int i, int j) const {
1588  return if_then_else( cond, val_type(0.0), expr2.fastAccessDx(i,j) );
1589  }
1590 
1591  protected:
1592 
1593  const CondT& cond;
1594  const typename ConstTypeRef<ConstT,value_type>::type c;
1595  const T2& expr2;
1596  };
1597 
1598  }
1599  }
1600 }
1601 
1602 #undef FAD_BINARYOP_MACRO
1603 
1604 #endif // SACADO_FAD_EXP_MP_VECTOR_HPP
expr expr ASinhOp
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > fabs(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > tan(const PCE< Storage > &a)
expr expr TanhOp
KOKKOS_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const ConstT &c_, const T2 &expr2_)
KOKKOS_INLINE_FUNCTION PCE< Storage > sinh(const PCE< Storage > &a)
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c expr1 expr2 expr1 expr2 expr1 MultiplicationOp
expr expr ATanhOp
atanh(expr.val())
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 j expr1 expr1 expr1 expr1 j expr1 c *expr2 expr1 c expr1 c expr1 c expr1 DivisionOp
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const ConstT &c_)
expr expr ASinOp
KOKKOS_INLINE_FUNCTION PCE< Storage > tanh(const PCE< Storage > &a)
#define FAD_BINARYOP_MACRO(OPNAME, OP, USING, MPVALUE, VALUE, DX, CDX1, CDX2, FASTACCESSDX, MPVAL_CONST_DX_1, MPVAL_CONST_DX_2, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
KOKKOS_INLINE_FUNCTION PCE< Storage > cbrt(const PCE< Storage > &a)
expr expr TanOp
KOKKOS_INLINE_FUNCTION PCE< Storage > acos(const PCE< Storage > &a)
atan2(expr1.val(), expr2.val())
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
KOKKOS_INLINE_FUNCTION const val_type & fastAccessDx(int i, int j) const
Returns derivative component i without bounds checking.
KOKKOS_INLINE_FUNCTION val_type dx(int i, int j) const
Returns derivative component i with bounds checking.
asinh(expr.val())
expr2 j expr1 expr1 expr2 expr2 j expr1 c c c c MaxOp
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 j expr1 expr1 expr1 expr1 j expr1 c *expr2 expr1 c expr1 c expr1 c expr1 expr1 expr1 expr1 j *expr1 expr2 expr1 expr1 j *expr1 c expr2 expr1 c expr1 expr2 expr1 expr2 expr1 Atan2Op
expr expr CoshOp
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
expr expr expr expr fastAccessDx(i, j)) FAD_UNARYOP_MACRO(exp
KOKKOS_INLINE_FUNCTION PCE< Storage > cosh(const PCE< Storage > &a)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, MPVALUE, VALUE, DX, FASTACCESSDX)
if_then_else(expr.val() >=0, expr.dx(i, j), value_type(-expr.dx(i, j)))
expr expr CosOp
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j)-expr2.val(j)
expr expr expr dx(i, j)
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
expr2 j expr1 expr1 expr2 expr2 j expr1 c c c c MinOp
expr expr SinhOp
expr val()
acosh(expr.val())
expr expr SinOp
KOKKOS_INLINE_FUNCTION val_type & fastAccessDx(int i, int j)
Returns derivative component i without bounds checking.
expr expr expr expr ExpOp
expr expr SqrtOp
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
expr expr AbsOp
expr expr ACoshOp
expr expr ATanOp
KOKKOS_INLINE_FUNCTION PCE< Storage > log(const PCE< Storage > &a)
Expression template specialization tag for Fad&lt; MP::Vector &gt;
KOKKOS_INLINE_FUNCTION PCE< Storage > log10(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > asin(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cos(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION IfThenElseOp(const CondT &cond_, const T1 &expr1_, const T2 &expr2_)
expr expr ACosOp