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_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 // The forward-mode AD classes in Sacado are a derivative work of the
30 // expression template classes in the Fad package by Nicolas Di Cesare.
31 // The following banner is included in the original Fad source code:
32 //
33 // ************ DO NOT REMOVE THIS BANNER ****************
34 //
35 // Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
36 // http://www.ann.jussieu.fr/~dicesare
37 //
38 // CEMRACS 98 : C++ courses,
39 // templates : new C++ techniques
40 // for scientific computing
41 //
42 //********************************************************
43 //
44 // A short implementation ( not all operators and
45 // functions are overloaded ) of 1st order Automatic
46 // Differentiation in forward mode (FAD) using
47 // EXPRESSION TEMPLATES.
48 //
49 //********************************************************
50 // @HEADER
51 
52 #ifndef SACADO_FAD_OPS_HPP
53 #define SACADO_FAD_OPS_HPP
54 
56 #include "Sacado_Fad_Ops_Fwd.hpp"
57 #include "Sacado_cmath.hpp"
58 #include <ostream> // for std::ostream
59 
60 #define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
61 namespace Sacado { \
62  namespace Fad { \
63  \
64  template <typename ExprT> \
65  class OP {}; \
66  \
67  template <typename ExprT> \
68  struct ExprSpec< OP<ExprT> > { \
69  typedef typename ExprSpec<ExprT>::type type; \
70  }; \
71  \
72  template <typename ExprT> \
73  class Expr< OP<ExprT>,ExprSpecDefault > { \
74  public: \
75  \
76  typedef typename ExprT::value_type value_type; \
77  typedef typename ExprT::scalar_type scalar_type; \
78  typedef typename ExprT::base_expr_type base_expr_type; \
79  \
80  KOKKOS_INLINE_FUNCTION \
81  explicit Expr(const ExprT& expr_) : expr(expr_) {} \
82  \
83  KOKKOS_INLINE_FUNCTION \
84  int size() const { return expr.size(); } \
85  \
86  KOKKOS_INLINE_FUNCTION \
87  bool hasFastAccess() const { return expr.hasFastAccess(); } \
88  \
89  KOKKOS_INLINE_FUNCTION \
90  bool isPassive() const { return expr.isPassive();} \
91  \
92  KOKKOS_INLINE_FUNCTION \
93  bool updateValue() const { return expr.updateValue(); } \
94  \
95  KOKKOS_INLINE_FUNCTION \
96  void cache() const {} \
97  \
98  KOKKOS_INLINE_FUNCTION \
99  value_type val() const { \
100  USING \
101  return VALUE; \
102  } \
103  \
104  KOKKOS_INLINE_FUNCTION \
105  value_type dx(int i) const { \
106  USING \
107  return DX; \
108  } \
109  \
110  KOKKOS_INLINE_FUNCTION \
111  value_type fastAccessDx(int i) const { \
112  USING \
113  return FASTACCESSDX; \
114  } \
115  \
116  protected: \
117  \
118  const ExprT& expr; \
119  }; \
120  \
121  template <typename T> \
122  KOKKOS_INLINE_FUNCTION \
123  Expr< OP< Expr<T> > > \
124  OPNAME (const Expr<T>& expr) \
125  { \
126  typedef OP< Expr<T> > expr_t; \
127  \
128  return Expr<expr_t>(expr); \
129  } \
130  } \
131  \
132 }
133 
134 FAD_UNARYOP_MACRO(operator+,
135  UnaryPlusOp,
136  ;,
137  expr.val(),
138  expr.dx(i),
139  expr.fastAccessDx(i))
140 FAD_UNARYOP_MACRO(operator-,
142  ;,
143  -expr.val(),
144  -expr.dx(i),
145  -expr.fastAccessDx(i))
148  using std::exp;,
149  exp(expr.val()),
150  exp(expr.val())*expr.dx(i),
151  exp(expr.val())*expr.fastAccessDx(i))
154  using std::log;,
155  log(expr.val()),
156  expr.dx(i)/expr.val(),
157  expr.fastAccessDx(i)/expr.val())
160  using std::log10; using std::log;,
161  log10(expr.val()),
162  expr.dx(i)/( log(value_type(10))*expr.val()),
163  expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
166  using std::sqrt;,
167  sqrt(expr.val()),
168  expr.dx(i)/(value_type(2)* sqrt(expr.val())),
169  expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
172  using std::cos; using std::sin;,
173  cos(expr.val()),
174  -expr.dx(i)* sin(expr.val()),
175  -expr.fastAccessDx(i)* sin(expr.val()))
178  using std::cos; using std::sin;,
179  sin(expr.val()),
180  expr.dx(i)* cos(expr.val()),
181  expr.fastAccessDx(i)* cos(expr.val()))
184  using std::tan;,
185  std::tan(expr.val()),
186  expr.dx(i)*
187  (value_type(1)+ tan(expr.val())* tan(expr.val())),
188  expr.fastAccessDx(i)*
189  (value_type(1)+ tan(expr.val())* tan(expr.val())))
192  using std::acos; using std::sqrt;,
193  acos(expr.val()),
194  -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
195  -expr.fastAccessDx(i) /
196  sqrt(value_type(1)-expr.val()*expr.val()))
199  using std::asin; using std::sqrt;,
200  asin(expr.val()),
201  expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
202  expr.fastAccessDx(i) /
203  sqrt(value_type(1)-expr.val()*expr.val()))
206  using std::atan;,
207  atan(expr.val()),
208  expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
209  expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
212  using std::cosh; using std::sinh;,
213  cosh(expr.val()),
214  expr.dx(i)* sinh(expr.val()),
215  expr.fastAccessDx(i)* sinh(expr.val()))
216 FAD_UNARYOP_MACRO(sinh,
218  using std::cosh; using std::sinh;,
219  sinh(expr.val()),
220  expr.dx(i)* cosh(expr.val()),
221  expr.fastAccessDx(i)* cosh(expr.val()))
224  using std::tanh;,
225  tanh(expr.val()),
226  expr.dx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())),
227  expr.fastAccessDx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())))
230  using std::acosh; using std::sqrt;,
231  acosh(expr.val()),
232  expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
233  (expr.val()+value_type(1))),
234  expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
235  (expr.val()+value_type(1))))
238  using std::asinh; using std::sqrt;,
239  asinh(expr.val()),
240  expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
241  expr.fastAccessDx(i)/ sqrt(value_type(1)+
242  expr.val()*expr.val()))
245  using std::atanh;,
246  atanh(expr.val()),
247  expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
248  expr.fastAccessDx(i)/(value_type(1)-
249  expr.val()*expr.val()))
252  using std::abs; using Sacado::if_then_else;,
253  abs(expr.val()),
254  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
255  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
258  using std::fabs; using Sacado::if_then_else;,
259  fabs(expr.val()),
260  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
261  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
262 #ifdef HAVE_SACADO_CXX11
264  CbrtOp,
265  using std::cbrt;,
266  cbrt(expr.val()),
267  expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
268  expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
269 #endif
270 
271 #undef FAD_UNARYOP_MACRO
272 
273 // Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
274 // "simd" value types that use if_then_else(). The only reason for not using
275 // if_then_else() always is to avoid evaluating the derivative if the value is
276 // zero to avoid throwing FPEs.
277 namespace Sacado {
278  namespace Fad {
279 
280  template <typename ExprT, bool is_simd>
281  class SafeSqrtOp {};
282 
283  template <typename ExprT>
284  struct ExprSpec< SafeSqrtOp<ExprT> > {
285  typedef typename ExprSpec<ExprT>::type type;
286  };
287 
288  //
289  // Implementation for simd type using if_then_else()
290  //
291  template <typename ExprT>
292  class Expr< SafeSqrtOp<ExprT,true>,ExprSpecDefault > {
293  public:
294 
295  typedef typename ExprT::value_type value_type;
296  typedef typename ExprT::scalar_type scalar_type;
297  typedef typename ExprT::base_expr_type base_expr_type;
298 
300  explicit Expr(const ExprT& expr_) : expr(expr_) {}
301 
303  int size() const { return expr.size(); }
304 
306  bool hasFastAccess() const { return expr.hasFastAccess(); }
307 
309  bool isPassive() const { return expr.isPassive();}
310 
312  bool updateValue() const { return expr.updateValue(); }
313 
315  void cache() const {}
316 
318  value_type val() const {
319  using std::sqrt;
320  return sqrt(expr.val());
321  }
322 
324  value_type dx(int i) const {
325  using std::sqrt;
326  return if_then_else(
327  expr.val() == value_type(0.0), value_type(0.0),
328  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
329  }
330 
332  value_type fastAccessDx(int i) const {
333  using std::sqrt;
334  return if_then_else(
335  expr.val() == value_type(0.0), value_type(0.0),
336  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
337  }
338 
339  protected:
340 
341  const ExprT& expr;
342  };
343 
344  //
345  // Specialization for scalar types using ternary operator
346  //
347  template <typename ExprT>
348  class Expr< SafeSqrtOp<ExprT,false>,ExprSpecDefault > {
349  public:
350 
351  typedef typename ExprT::value_type value_type;
352  typedef typename ExprT::scalar_type scalar_type;
353  typedef typename ExprT::base_expr_type base_expr_type;
354 
356  explicit Expr(const ExprT& expr_) : expr(expr_) {}
357 
359  int size() const { return expr.size(); }
360 
362  bool hasFastAccess() const { return expr.hasFastAccess(); }
363 
365  bool isPassive() const { return expr.isPassive();}
366 
368  bool updateValue() const { return expr.updateValue(); }
369 
371  void cache() const {}
372 
374  value_type val() const {
375  using std::sqrt;
376  return sqrt(expr.val());
377  }
378 
380  value_type dx(int i) const {
381  using std::sqrt;
382  return expr.val() == value_type(0.0) ? value_type(0.0) :
383  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
384  }
385 
387  value_type fastAccessDx(int i) const {
388  using std::sqrt;
389  return expr.val() == value_type(0.0) ? value_type(0.0) :
390  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
391  }
392 
393  protected:
394 
395  const ExprT& expr;
396  };
397 
398  template <typename T>
400  Expr< SafeSqrtOp< Expr<T> > >
401  safe_sqrt (const Expr<T>& expr)
402  {
403  typedef SafeSqrtOp< Expr<T> > expr_t;
404 
405  return Expr<expr_t>(expr);
406  }
407  }
408 
409 }
410 
411 #define FAD_BINARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX,VAL_CONST_DX_1,VAL_CONST_DX_2,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
412 namespace Sacado { \
413  namespace Fad { \
414  \
415  template <typename ExprT1, typename ExprT2> \
416  class OP {}; \
417  \
418  template <typename ExprT1, typename ExprT2> \
419  struct ExprSpec< OP< ExprT1, ExprT2 > > { \
420  typedef typename ExprSpec<ExprT1>::type type; \
421  }; \
422  \
423  template <typename ExprT1, typename ExprT2> \
424  class Expr< OP< ExprT1, ExprT2 >,ExprSpecDefault > { \
425  \
426  public: \
427  \
428  typedef typename ExprT1::value_type value_type_1; \
429  typedef typename ExprT2::value_type value_type_2; \
430  typedef typename Sacado::Promote<value_type_1, \
431  value_type_2>::type value_type; \
432  \
433  typedef typename ExprT1::scalar_type scalar_type_1; \
434  typedef typename ExprT2::scalar_type scalar_type_2; \
435  typedef typename Sacado::Promote<scalar_type_1, \
436  scalar_type_2>::type scalar_type; \
437  \
438  typedef typename ExprT1::base_expr_type base_expr_type_1; \
439  typedef typename ExprT2::base_expr_type base_expr_type_2; \
440  typedef typename Sacado::Promote<base_expr_type_1, \
441  base_expr_type_2>::type base_expr_type; \
442  \
443  KOKKOS_INLINE_FUNCTION \
444  Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
445  expr1(expr1_), expr2(expr2_) {} \
446  \
447  KOKKOS_INLINE_FUNCTION \
448  int size() const { \
449  int sz1 = expr1.size(), sz2 = expr2.size(); \
450  return sz1 > sz2 ? sz1 : sz2; \
451  } \
452  \
453  KOKKOS_INLINE_FUNCTION \
454  bool hasFastAccess() const { \
455  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
456  } \
457  \
458  KOKKOS_INLINE_FUNCTION \
459  bool isPassive() const { \
460  return expr1.isPassive() && expr2.isPassive(); \
461  } \
462  \
463  KOKKOS_INLINE_FUNCTION \
464  bool updateValue() const { \
465  return expr1.updateValue() && expr2.updateValue(); \
466  } \
467  \
468  KOKKOS_INLINE_FUNCTION \
469  void cache() const {} \
470  \
471  KOKKOS_INLINE_FUNCTION \
472  const value_type val() const { \
473  USING \
474  return VALUE; \
475  } \
476  \
477  KOKKOS_INLINE_FUNCTION \
478  const value_type dx(int i) const { \
479  USING \
480  return DX; \
481  } \
482  \
483  KOKKOS_INLINE_FUNCTION \
484  const value_type fastAccessDx(int i) const { \
485  USING \
486  return FASTACCESSDX; \
487  } \
488  \
489  protected: \
490  \
491  const ExprT1& expr1; \
492  const ExprT2& expr2; \
493  \
494  }; \
495  \
496  template <typename ExprT1, typename T2> \
497  struct ExprSpec< OP< ExprT1, ConstExpr<T2> > > { \
498  typedef typename ExprSpec<ExprT1>::type type; \
499  }; \
500  \
501  template <typename ExprT1, typename T2> \
502  class Expr< OP< ExprT1, ConstExpr<T2> >,ExprSpecDefault > { \
503  \
504  public: \
505  \
506  typedef ConstExpr<T2> ConstT; \
507  typedef ConstExpr<T2> ExprT2; \
508  typedef typename ExprT1::value_type value_type_1; \
509  typedef typename ExprT2::value_type value_type_2; \
510  typedef typename Sacado::Promote<value_type_1, \
511  value_type_2>::type value_type; \
512  \
513  typedef typename ExprT1::scalar_type scalar_type_1; \
514  typedef typename ExprT2::scalar_type scalar_type_2; \
515  typedef typename Sacado::Promote<scalar_type_1, \
516  scalar_type_2>::type scalar_type; \
517  \
518  typedef typename ExprT1::base_expr_type base_expr_type_1; \
519  typedef typename ExprT2::base_expr_type base_expr_type_2; \
520  typedef typename Sacado::Promote<base_expr_type_1, \
521  base_expr_type_2>::type base_expr_type; \
522  \
523  KOKKOS_INLINE_FUNCTION \
524  Expr(const ExprT1& expr1_, const ConstT& c_) : \
525  expr1(expr1_), c(c_) {} \
526  \
527  KOKKOS_INLINE_FUNCTION \
528  int size() const { \
529  return expr1.size(); \
530  } \
531  \
532  KOKKOS_INLINE_FUNCTION \
533  bool hasFastAccess() const { \
534  return expr1.hasFastAccess(); \
535  } \
536  \
537  KOKKOS_INLINE_FUNCTION \
538  bool isPassive() const { \
539  return expr1.isPassive(); \
540  } \
541  \
542  KOKKOS_INLINE_FUNCTION \
543  bool updateValue() const { return expr1.updateValue(); } \
544  \
545  KOKKOS_INLINE_FUNCTION \
546  void cache() const {} \
547  \
548  KOKKOS_INLINE_FUNCTION \
549  const value_type val() const { \
550  USING \
551  return VAL_CONST_DX_2; \
552  } \
553  \
554  KOKKOS_INLINE_FUNCTION \
555  const value_type dx(int i) const { \
556  USING \
557  return CONST_DX_2; \
558  } \
559  \
560  KOKKOS_INLINE_FUNCTION \
561  const value_type fastAccessDx(int i) const { \
562  USING \
563  return CONST_FASTACCESSDX_2; \
564  } \
565  \
566  protected: \
567  \
568  const ExprT1& expr1; \
569  ConstT c; \
570  }; \
571  \
572  template <typename T1, typename ExprT2> \
573  struct ExprSpec< OP< ConstExpr<T1>, ExprT2 > > { \
574  typedef typename ExprSpec<ExprT2>::type type; \
575  }; \
576  \
577  template <typename T1, typename ExprT2> \
578  class Expr< OP< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > { \
579  \
580  public: \
581  \
582  typedef ConstExpr<T1> ConstT; \
583  typedef ConstExpr<T1> ExprT1; \
584  typedef typename ExprT1::value_type value_type_1; \
585  typedef typename ExprT2::value_type value_type_2; \
586  typedef typename Sacado::Promote<value_type_1, \
587  value_type_2>::type value_type; \
588  \
589  typedef typename ExprT1::scalar_type scalar_type_1; \
590  typedef typename ExprT2::scalar_type scalar_type_2; \
591  typedef typename Sacado::Promote<scalar_type_1, \
592  scalar_type_2>::type scalar_type; \
593  \
594  typedef typename ExprT1::base_expr_type base_expr_type_1; \
595  typedef typename ExprT2::base_expr_type base_expr_type_2; \
596  typedef typename Sacado::Promote<base_expr_type_1, \
597  base_expr_type_2>::type base_expr_type; \
598  \
599  \
600  KOKKOS_INLINE_FUNCTION \
601  Expr(const ConstT& c_, const ExprT2& expr2_) : \
602  c(c_), expr2(expr2_) {} \
603  \
604  KOKKOS_INLINE_FUNCTION \
605  int size() const { \
606  return expr2.size(); \
607  } \
608  \
609  KOKKOS_INLINE_FUNCTION \
610  bool hasFastAccess() const { \
611  return expr2.hasFastAccess(); \
612  } \
613  \
614  KOKKOS_INLINE_FUNCTION \
615  bool isPassive() const { \
616  return expr2.isPassive(); \
617  } \
618  \
619  KOKKOS_INLINE_FUNCTION \
620  bool updateValue() const { return expr2.updateValue(); } \
621  \
622  KOKKOS_INLINE_FUNCTION \
623  void cache() const {} \
624  \
625  KOKKOS_INLINE_FUNCTION \
626  const value_type val() const { \
627  USING \
628  return VAL_CONST_DX_1; \
629  } \
630  \
631  KOKKOS_INLINE_FUNCTION \
632  const value_type dx(int i) const { \
633  USING \
634  return CONST_DX_1; \
635  } \
636  \
637  KOKKOS_INLINE_FUNCTION \
638  const value_type fastAccessDx(int i) const { \
639  USING \
640  return CONST_FASTACCESSDX_1; \
641  } \
642  \
643  protected: \
644  \
645  ConstT c; \
646  const ExprT2& expr2; \
647  }; \
648  \
649  template <typename T1, typename T2> \
650  KOKKOS_INLINE_FUNCTION \
651  typename mpl::enable_if_c< \
652  ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value, \
653  Expr< OP< Expr<T1>, Expr<T2> > > \
654  >::type \
655  /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(OP)*/ \
656  OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2) \
657  { \
658  typedef OP< Expr<T1>, Expr<T2> > expr_t; \
659  \
660  return Expr<expr_t>(expr1, expr2); \
661  } \
662  \
663  template <typename T> \
664  KOKKOS_INLINE_FUNCTION \
665  Expr< OP< Expr<T>, Expr<T> > > \
666  OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
667  { \
668  typedef OP< Expr<T>, Expr<T> > expr_t; \
669  \
670  return Expr<expr_t>(expr1, expr2); \
671  } \
672  \
673  template <typename T> \
674  KOKKOS_INLINE_FUNCTION \
675  Expr< OP< ConstExpr<typename Expr<T>::value_type>, \
676  Expr<T> > > \
677  OPNAME (const typename Expr<T>::value_type& c, \
678  const Expr<T>& expr) \
679  { \
680  typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
681  typedef OP< ConstT, Expr<T> > expr_t; \
682  \
683  return Expr<expr_t>(ConstT(c), expr); \
684  } \
685  \
686  template <typename T> \
687  KOKKOS_INLINE_FUNCTION \
688  Expr< OP< Expr<T>, \
689  ConstExpr<typename Expr<T>::value_type> > > \
690  OPNAME (const Expr<T>& expr, \
691  const typename Expr<T>::value_type& c) \
692  { \
693  typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
694  typedef OP< Expr<T>, ConstT > expr_t; \
695  \
696  return Expr<expr_t>(expr, ConstT(c)); \
697  } \
698  \
699  template <typename T> \
700  KOKKOS_INLINE_FUNCTION \
701  SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP) \
702  OPNAME (const typename Expr<T>::scalar_type& c, \
703  const Expr<T>& expr) \
704  { \
705  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
706  typedef OP< ConstT, Expr<T> > expr_t; \
707  \
708  return Expr<expr_t>(ConstT(c), expr); \
709  } \
710  \
711  template <typename T> \
712  KOKKOS_INLINE_FUNCTION \
713  SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP) \
714  OPNAME (const Expr<T>& expr, \
715  const typename Expr<T>::scalar_type& c) \
716  { \
717  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
718  typedef OP< Expr<T>, ConstT > expr_t; \
719  \
720  return Expr<expr_t>(expr, ConstT(c)); \
721  } \
722  } \
723  \
724 }
725 
726 
727 FAD_BINARYOP_MACRO(operator+,
728  AdditionOp,
729  ;,
730  expr1.val() + expr2.val(),
731  expr1.dx(i) + expr2.dx(i),
732  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
733  c.val() + expr2.val(),
734  expr1.val() + c.val(),
735  expr2.dx(i),
736  expr1.dx(i),
737  expr2.fastAccessDx(i),
738  expr1.fastAccessDx(i))
739 FAD_BINARYOP_MACRO(operator-,
741  ;,
742  expr1.val() - expr2.val(),
743  expr1.dx(i) - expr2.dx(i),
744  expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
745  c.val() - expr2.val(),
746  expr1.val() - c.val(),
747  -expr2.dx(i),
748  expr1.dx(i),
749  -expr2.fastAccessDx(i),
750  expr1.fastAccessDx(i))
751 // FAD_BINARYOP_MACRO(operator*,
752 // MultiplicationOp,
753 // ;,
754 // expr1.val() * expr2.val(),
755 // expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
756 // expr1.val()*expr2.fastAccessDx(i) +
757 // expr1.fastAccessDx(i)*expr2.val(),
758 // c.val() * expr2.val(),
759 // expr1.val() * c.val(),
760 // c.val()*expr2.dx(i),
761 // expr1.dx(i)*c.val(),
762 // c.val()*expr2.fastAccessDx(i),
763 // expr1.fastAccessDx(i)*c.val())
764 FAD_BINARYOP_MACRO(operator/,
766  ;,
767  expr1.val() / expr2.val(),
768  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
769  (expr2.val()*expr2.val()),
770  (expr1.fastAccessDx(i)*expr2.val() -
771  expr2.fastAccessDx(i)*expr1.val()) /
772  (expr2.val()*expr2.val()),
773  c.val() / expr2.val(),
774  expr1.val() / c.val(),
775  -expr2.dx(i)*c.val() / (expr2.val()*expr2.val()),
776  expr1.dx(i)/c.val(),
777  -expr2.fastAccessDx(i)*c.val() / (expr2.val()*expr2.val()),
778  expr1.fastAccessDx(i)/c.val())
781  using std::atan2;,
782  atan2(expr1.val(), expr2.val()),
783  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
784  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
785  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
786  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
787  atan2(c.val(), expr2.val()),
788  atan2(expr1.val(), c.val()),
789  (-c.val()*expr2.dx(i)) / (c.val()*c.val() + expr2.val()*expr2.val()),
790  (c.val()*expr1.dx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()),
791  (-c.val()*expr2.fastAccessDx(i))/ (c.val()*c.val() + expr2.val()*expr2.val()),
792  (c.val()*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()))
793 // FAD_BINARYOP_MACRO(pow,
794 // PowerOp,
795 // using std::pow; using std::log; using Sacado::if_then_else;,
796 // pow(expr1.val(), expr2.val()),
797 // 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())) ),
798 // 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())) ),
799 // pow(c.val(), expr2.val()),
800 // pow(expr1.val(), c.val()),
801 // if_then_else( c.val() == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val())) ),
802 // if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val())) ),
803 // if_then_else( c.val() == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val())) ),
804 // if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()))) )
807  using Sacado::if_then_else;,
808  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
809  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
810  if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
811  if_then_else( c.val() >= expr2.val(), value_type(c.val()), expr2.val() ),
812  if_then_else( expr1.val() >= c.val(), expr1.val(), value_type(c.val()) ),
813  if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
814  if_then_else( expr1.val() >= c.val(), expr1.dx(i), value_type(0.0) ),
815  if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
816  if_then_else( expr1.val() >= c.val(), expr1.fastAccessDx(i), value_type(0.0) ) )
819  using Sacado::if_then_else;,
820  if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
821  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
822  if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
823  if_then_else( c.val() <= expr2.val(), value_type(c.val()), expr2.val() ),
824  if_then_else( expr1.val() <= c.val(), expr1.val(), value_type(c.val()) ),
825  if_then_else( c.val() <= expr2.val(), value_type(0), expr2.dx(i) ),
826  if_then_else( expr1.val() <= c.val(), expr1.dx(i), value_type(0) ),
827  if_then_else( c.val() <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
828  if_then_else( expr1.val() <= c.val(), expr1.fastAccessDx(i), value_type(0) ) )
829 
830 
831 #undef FAD_BINARYOP_MACRO
832 
833 namespace Sacado {
834  namespace Fad {
835 
836  template <typename ExprT1, typename ExprT2>
837  class MultiplicationOp {};
838 
839  template <typename ExprT1, typename ExprT2>
840  struct ExprSpec< MultiplicationOp< ExprT1, ExprT2 > > {
841  typedef typename ExprSpec<ExprT1>::type type;
842  };
843 
844  template <typename ExprT1, typename ExprT2>
845  class Expr< MultiplicationOp< ExprT1, ExprT2 >,ExprSpecDefault > {
846 
847  public:
848 
849  typedef typename ExprT1::value_type value_type_1;
850  typedef typename ExprT2::value_type value_type_2;
851  typedef typename Sacado::Promote<value_type_1,
852  value_type_2>::type value_type;
853 
854  typedef typename ExprT1::scalar_type scalar_type_1;
855  typedef typename ExprT2::scalar_type scalar_type_2;
856  typedef typename Sacado::Promote<scalar_type_1,
857  scalar_type_2>::type scalar_type;
858 
859  typedef typename ExprT1::base_expr_type base_expr_type_1;
860  typedef typename ExprT2::base_expr_type base_expr_type_2;
861  typedef typename Sacado::Promote<base_expr_type_1,
862  base_expr_type_2>::type base_expr_type;
863 
865  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
866  expr1(expr1_), expr2(expr2_) {}
867 
869  int size() const {
870  int sz1 = expr1.size(), sz2 = expr2.size();
871  return sz1 > sz2 ? sz1 : sz2;
872  }
873 
875  bool hasFastAccess() const {
876  return expr1.hasFastAccess() && expr2.hasFastAccess();
877  }
878 
880  bool isPassive() const {
881  return expr1.isPassive() && expr2.isPassive();
882  }
883 
885  bool updateValue() const {
886  return expr1.updateValue() && expr2.updateValue();
887  }
888 
890  void cache() const {}
891 
893  const value_type val() const {
894  return expr1.val()*expr2.val();
895  }
896 
898  const value_type dx(int i) const {
899  if (expr1.size() > 0 && expr2.size() > 0)
900  return expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val();
901  else if (expr1.size() > 0)
902  return expr1.dx(i)*expr2.val();
903  else
904  return expr1.val()*expr2.dx(i);
905  }
906 
908  const value_type fastAccessDx(int i) const {
909  return expr1.val()*expr2.fastAccessDx(i) +
910  expr1.fastAccessDx(i)*expr2.val();
911  }
912 
913  protected:
914 
915  const ExprT1& expr1;
916  const ExprT2& expr2;
917 
918  };
919 
920  template <typename ExprT1, typename T2>
921  struct ExprSpec< MultiplicationOp< ExprT1, ConstExpr<T2> > > {
922  typedef typename ExprSpec<ExprT1>::type type;
923  };
924 
925  template <typename ExprT1, typename T2>
926  class Expr< MultiplicationOp< ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
927 
928  public:
929 
930  typedef ConstExpr<T2> ConstT;
931  typedef ConstExpr<T2> ExprT2;
932  typedef typename ExprT1::value_type value_type_1;
933  typedef typename ExprT2::value_type value_type_2;
934  typedef typename Sacado::Promote<value_type_1,
935  value_type_2>::type value_type;
936 
937  typedef typename ExprT1::scalar_type scalar_type_1;
938  typedef typename ExprT2::scalar_type scalar_type_2;
939  typedef typename Sacado::Promote<scalar_type_1,
940  scalar_type_2>::type scalar_type;
941 
942  typedef typename ExprT1::base_expr_type base_expr_type_1;
943  typedef typename ExprT2::base_expr_type base_expr_type_2;
944  typedef typename Sacado::Promote<base_expr_type_1,
945  base_expr_type_2>::type base_expr_type;
946 
948  Expr(const ExprT1& expr1_, const ConstT& c_) :
949  expr1(expr1_), c(c_) {}
950 
952  int size() const {
953  return expr1.size();
954  }
955 
957  bool hasFastAccess() const {
958  return expr1.hasFastAccess();
959  }
960 
962  bool isPassive() const {
963  return expr1.isPassive();
964  }
965 
967  bool updateValue() const { return expr1.updateValue(); }
968 
970  void cache() const {}
971 
973  const value_type val() const {
974  return expr1.val()*c.val();
975  }
976 
978  const value_type dx(int i) const {
979  return expr1.dx(i)*c.val();
980  }
981 
983  const value_type fastAccessDx(int i) const {
984  return expr1.fastAccessDx(i)*c.val();
985  }
986 
987  protected:
988 
989  const ExprT1& expr1;
990  ConstT c;
991  };
992 
993  template <typename T1, typename ExprT2>
994  struct ExprSpec< MultiplicationOp< ConstExpr<T1>, ExprT2 > > {
995  typedef typename ExprSpec<ExprT2>::type type;
996  };
997 
998  template <typename T1, typename ExprT2>
999  class Expr< MultiplicationOp< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
1000 
1001  public:
1002 
1003  typedef ConstExpr<T1> ConstT;
1004  typedef ConstExpr<T1> ExprT1;
1005  typedef typename ExprT1::value_type value_type_1;
1006  typedef typename ExprT2::value_type value_type_2;
1007  typedef typename Sacado::Promote<value_type_1,
1008  value_type_2>::type value_type;
1009 
1010  typedef typename ExprT1::scalar_type scalar_type_1;
1011  typedef typename ExprT2::scalar_type scalar_type_2;
1012  typedef typename Sacado::Promote<scalar_type_1,
1013  scalar_type_2>::type scalar_type;
1014 
1015  typedef typename ExprT1::base_expr_type base_expr_type_1;
1016  typedef typename ExprT2::base_expr_type base_expr_type_2;
1017  typedef typename Sacado::Promote<base_expr_type_1,
1018  base_expr_type_2>::type base_expr_type;
1019 
1021  Expr(const ConstT& c_, const ExprT2& expr2_) :
1022  c(c_), expr2(expr2_) {}
1023 
1025  int size() const {
1026  return expr2.size();
1027  }
1028 
1030  bool hasFastAccess() const {
1031  return expr2.hasFastAccess();
1032  }
1033 
1035  bool isPassive() const {
1036  return expr2.isPassive();
1037  }
1038 
1040  bool updateValue() const { return expr2.updateValue(); }
1041 
1043  void cache() const {}
1044 
1046  const value_type val() const {
1047  return c.val()*expr2.val();
1048  }
1049 
1051  const value_type dx(int i) const {
1052  return c.val()*expr2.dx(i);
1053  }
1054 
1056  const value_type fastAccessDx(int i) const {
1057  return c.val()*expr2.fastAccessDx(i);
1058  }
1059 
1060  protected:
1061 
1062  ConstT c;
1063  const ExprT2& expr2;
1064  };
1065 
1066  template <typename T1, typename T2>
1068  typename mpl::enable_if_c<
1069  ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value,
1070  Expr< MultiplicationOp< Expr<T1>, Expr<T2> > >
1071  >::type
1072  /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(MultiplicationOp)*/
1073  operator* (const Expr<T1>& expr1, const Expr<T2>& expr2)
1074  {
1075  typedef MultiplicationOp< Expr<T1>, Expr<T2> > expr_t;
1076 
1077  return Expr<expr_t>(expr1, expr2);
1078  }
1079 
1080  template <typename T>
1082  Expr< MultiplicationOp< Expr<T>, Expr<T> > >
1083  operator* (const Expr<T>& expr1, const Expr<T>& expr2)
1084  {
1085  typedef MultiplicationOp< Expr<T>, Expr<T> > expr_t;
1086 
1087  return Expr<expr_t>(expr1, expr2);
1088  }
1089 
1090  template <typename T>
1092  Expr< MultiplicationOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
1093  operator* (const typename Expr<T>::value_type& c,
1094  const Expr<T>& expr)
1095  {
1096  typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1097  typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1098 
1099  return Expr<expr_t>(ConstT(c), expr);
1100  }
1101 
1102  template <typename T>
1104  Expr< MultiplicationOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
1105  operator* (const Expr<T>& expr,
1106  const typename Expr<T>::value_type& c)
1107  {
1108  typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1109  typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1110 
1111  return Expr<expr_t>(expr, ConstT(c));
1112  }
1113 
1114  template <typename T>
1116  SACADO_FAD_OP_ENABLE_SCALAR_EXPR(MultiplicationOp)
1117  operator* (const typename Expr<T>::scalar_type& c,
1118  const Expr<T>& expr)
1119  {
1120  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1121  typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1122 
1123  return Expr<expr_t>(ConstT(c), expr);
1124  }
1125 
1126  template <typename T>
1128  SACADO_FAD_OP_ENABLE_EXPR_SCALAR(MultiplicationOp)
1129  operator* (const Expr<T>& expr,
1130  const typename Expr<T>::scalar_type& c)
1131  {
1132  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1133  typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1134 
1135  return Expr<expr_t>(expr, ConstT(c));
1136  }
1137  }
1138 }
1139 
1140 // Special handling for std::pow() to provide specializations of PowerOp for
1141 // "simd" value types that use if_then_else(). The only reason for not using
1142 // if_then_else() always is to avoid evaluating the derivative if the value is
1143 // zero to avoid throwing FPEs.
1144 namespace Sacado {
1145  namespace Fad {
1146 
1147  template <typename ExprT1, typename ExprT2, bool is_simd>
1148  class PowerOp {};
1149 
1150  template <typename ExprT1, typename ExprT2>
1151  struct ExprSpec< PowerOp< ExprT1, ExprT2 > > {
1152  typedef typename ExprSpec<ExprT1>::type type;
1153  };
1154 
1155  template <typename ExprT1, typename T2>
1156  struct ExprSpec<PowerOp< ExprT1, ConstExpr<T2> > > {
1157  typedef typename ExprSpec<ExprT1>::type type;
1158  };
1159 
1160  template <typename T1, typename ExprT2>
1161  struct ExprSpec< PowerOp< ConstExpr<T1>, ExprT2 > > {
1162  typedef typename ExprSpec<ExprT2>::type type;
1163  };
1164 
1165  //
1166  // Implementation for simd type using if_then_else()
1167  //
1168  template <typename ExprT1, typename ExprT2>
1169  class Expr< PowerOp< ExprT1, ExprT2, true >, ExprSpecDefault > {
1170 
1171  public:
1172 
1173  typedef typename ExprT1::value_type value_type_1;
1174  typedef typename ExprT2::value_type value_type_2;
1175  typedef typename Sacado::Promote<value_type_1,
1177 
1178  typedef typename ExprT1::scalar_type scalar_type_1;
1179  typedef typename ExprT2::scalar_type scalar_type_2;
1180  typedef typename Sacado::Promote<scalar_type_1,
1182 
1183  typedef typename ExprT1::base_expr_type base_expr_type_1;
1184  typedef typename ExprT2::base_expr_type base_expr_type_2;
1185  typedef typename Sacado::Promote<base_expr_type_1,
1187 
1189  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1190  expr1(expr1_), expr2(expr2_) {}
1191 
1193  int size() const {
1194  int sz1 = expr1.size(), sz2 = expr2.size();
1195  return sz1 > sz2 ? sz1 : sz2;
1196  }
1197 
1199  bool hasFastAccess() const {
1200  return expr1.hasFastAccess() && expr2.hasFastAccess();
1201  }
1202 
1204  bool isPassive() const {
1205  return expr1.isPassive() && expr2.isPassive();
1206  }
1207 
1209  bool updateValue() const {
1210  return expr1.updateValue() && expr2.updateValue();
1211  }
1212 
1214  void cache() const {}
1215 
1217  const value_type val() const {
1218  using std::pow;
1219  return pow(expr1.val(), expr2.val());
1220  }
1221 
1223  const value_type dx(int i) const {
1224  using std::pow; using std::log;
1225  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())) );
1226  }
1227 
1229  const value_type fastAccessDx(int i) const {
1230  using std::pow; using std::log;
1231  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())) );
1232  }
1233 
1234  protected:
1235 
1236  const ExprT1& expr1;
1237  const ExprT2& expr2;
1238 
1239  };
1240 
1241  template <typename ExprT1, typename T2>
1242  class Expr< PowerOp< ExprT1, ConstExpr<T2>, true >, ExprSpecDefault > {
1243 
1244  public:
1245 
1248  typedef typename ExprT1::value_type value_type_1;
1250  typedef typename Sacado::Promote<value_type_1,
1252 
1253  typedef typename ExprT1::scalar_type scalar_type_1;
1255  typedef typename Sacado::Promote<scalar_type_1,
1257 
1258  typedef typename ExprT1::base_expr_type base_expr_type_1;
1260  typedef typename Sacado::Promote<base_expr_type_1,
1262 
1264  Expr(const ExprT1& expr1_, const ConstT& c_) :
1265  expr1(expr1_), c(c_) {}
1266 
1268  int size() const {
1269  return expr1.size();
1270  }
1271 
1273  bool hasFastAccess() const {
1274  return expr1.hasFastAccess();
1275  }
1276 
1278  bool isPassive() const {
1279  return expr1.isPassive();
1280  }
1281 
1283  bool updateValue() const { return expr1.updateValue(); }
1284 
1286  void cache() const {}
1287 
1289  const value_type val() const {
1290  using std::pow;
1291  return pow(expr1.val(), c.val());
1292  }
1293 
1295  const value_type dx(int i) const {
1296  using std::pow;
1297  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1298  // It seems less accurate and caused convergence problems in some codes
1299  return if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val())) );
1300  }
1301 
1303  const value_type fastAccessDx(int i) const {
1304  using std::pow;
1305  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1306  // It seems less accurate and caused convergence problems in some codes
1307  return if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val())));
1308  }
1309 
1310  protected:
1311 
1312  const ExprT1& expr1;
1314  };
1315 
1316  template <typename T1, typename ExprT2>
1317  class Expr< PowerOp< ConstExpr<T1>, ExprT2, true >, ExprSpecDefault > {
1318 
1319  public:
1320 
1324  typedef typename ExprT2::value_type value_type_2;
1325  typedef typename Sacado::Promote<value_type_1,
1327 
1329  typedef typename ExprT2::scalar_type scalar_type_2;
1330  typedef typename Sacado::Promote<scalar_type_1,
1332 
1334  typedef typename ExprT2::base_expr_type base_expr_type_2;
1335  typedef typename Sacado::Promote<base_expr_type_1,
1337 
1338 
1340  Expr(const ConstT& c_, const ExprT2& expr2_) :
1341  c(c_), expr2(expr2_) {}
1342 
1344  int size() const {
1345  return expr2.size();
1346  }
1347 
1349  bool hasFastAccess() const {
1350  return expr2.hasFastAccess();
1351  }
1352 
1354  bool isPassive() const {
1355  return expr2.isPassive();
1356  }
1357 
1359  bool updateValue() const { return expr2.updateValue(); }
1360 
1362  void cache() const {}
1363 
1365  const value_type val() const {
1366  using std::pow;
1367  return pow(c.val(), expr2.val());
1368  }
1369 
1371  const value_type dx(int i) const {
1372  using std::pow; using std::log;
1373  return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val())) );
1374  }
1375 
1377  const value_type fastAccessDx(int i) const {
1378  using std::pow; using std::log;
1379  return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val())) );
1380  }
1381 
1382  protected:
1383 
1385  const ExprT2& expr2;
1386  };
1387 
1388  //
1389  // Specialization for scalar types using ternary operator
1390  //
1391 
1392  template <typename ExprT1, typename ExprT2>
1393  class Expr< PowerOp< ExprT1, ExprT2, false >, ExprSpecDefault > {
1394 
1395  public:
1396 
1397  typedef typename ExprT1::value_type value_type_1;
1398  typedef typename ExprT2::value_type value_type_2;
1399  typedef typename Sacado::Promote<value_type_1,
1401 
1402  typedef typename ExprT1::scalar_type scalar_type_1;
1403  typedef typename ExprT2::scalar_type scalar_type_2;
1404  typedef typename Sacado::Promote<scalar_type_1,
1406 
1407  typedef typename ExprT1::base_expr_type base_expr_type_1;
1408  typedef typename ExprT2::base_expr_type base_expr_type_2;
1409  typedef typename Sacado::Promote<base_expr_type_1,
1411 
1413  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1414  expr1(expr1_), expr2(expr2_) {}
1415 
1417  int size() const {
1418  int sz1 = expr1.size(), sz2 = expr2.size();
1419  return sz1 > sz2 ? sz1 : sz2;
1420  }
1421 
1423  bool hasFastAccess() const {
1424  return expr1.hasFastAccess() && expr2.hasFastAccess();
1425  }
1426 
1428  bool isPassive() const {
1429  return expr1.isPassive() && expr2.isPassive();
1430  }
1431 
1433  bool updateValue() const {
1434  return expr1.updateValue() && expr2.updateValue();
1435  }
1436 
1438  void cache() const {}
1439 
1441  const value_type val() const {
1442  using std::pow;
1443  return pow(expr1.val(), expr2.val());
1444  }
1445 
1447  const value_type dx(int i) const {
1448  using std::pow; using std::log;
1449  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()));
1450  }
1451 
1453  const value_type fastAccessDx(int i) const {
1454  using std::pow; using std::log;
1455  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()));
1456  }
1457 
1458  protected:
1459 
1460  const ExprT1& expr1;
1461  const ExprT2& expr2;
1462 
1463  };
1464 
1465  template <typename ExprT1, typename T2>
1466  class Expr< PowerOp< ExprT1, ConstExpr<T2>, false >, ExprSpecDefault > {
1467 
1468  public:
1469 
1472  typedef typename ExprT1::value_type value_type_1;
1474  typedef typename Sacado::Promote<value_type_1,
1476 
1477  typedef typename ExprT1::scalar_type scalar_type_1;
1479  typedef typename Sacado::Promote<scalar_type_1,
1481 
1482  typedef typename ExprT1::base_expr_type base_expr_type_1;
1484  typedef typename Sacado::Promote<base_expr_type_1,
1486 
1488  Expr(const ExprT1& expr1_, const ConstT& c_) :
1489  expr1(expr1_), c(c_) {}
1490 
1492  int size() const {
1493  return expr1.size();
1494  }
1495 
1497  bool hasFastAccess() const {
1498  return expr1.hasFastAccess();
1499  }
1500 
1502  bool isPassive() const {
1503  return expr1.isPassive();
1504  }
1505 
1507  bool updateValue() const { return expr1.updateValue(); }
1508 
1510  void cache() const {}
1511 
1513  const value_type val() const {
1514  using std::pow;
1515  return pow(expr1.val(), c.val());
1516  }
1517 
1519  const value_type dx(int i) const {
1520  using std::pow;
1521  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1522  // It seems less accurate and caused convergence problems in some codes
1523  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val()));
1524  }
1525 
1527  const value_type fastAccessDx(int i) const {
1528  using std::pow;
1529  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1530  // It seems less accurate and caused convergence problems in some codes
1531  return expr1.val() == value_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()));
1532  }
1533 
1534  protected:
1535 
1536  const ExprT1& expr1;
1538  };
1539 
1540  template <typename T1, typename ExprT2>
1541  class Expr< PowerOp< ConstExpr<T1>, ExprT2, false >, ExprSpecDefault > {
1542 
1543  public:
1544 
1548  typedef typename ExprT2::value_type value_type_2;
1549  typedef typename Sacado::Promote<value_type_1,
1551 
1553  typedef typename ExprT2::scalar_type scalar_type_2;
1554  typedef typename Sacado::Promote<scalar_type_1,
1556 
1558  typedef typename ExprT2::base_expr_type base_expr_type_2;
1559  typedef typename Sacado::Promote<base_expr_type_1,
1561 
1562 
1564  Expr(const ConstT& c_, const ExprT2& expr2_) :
1565  c(c_), expr2(expr2_) {}
1566 
1568  int size() const {
1569  return expr2.size();
1570  }
1571 
1573  bool hasFastAccess() const {
1574  return expr2.hasFastAccess();
1575  }
1576 
1578  bool isPassive() const {
1579  return expr2.isPassive();
1580  }
1581 
1583  bool updateValue() const { return expr2.updateValue(); }
1584 
1586  void cache() const {}
1587 
1589  const value_type val() const {
1590  using std::pow;
1591  return pow(c.val(), expr2.val());
1592  }
1593 
1595  const value_type dx(int i) const {
1596  using std::pow; using std::log;
1597  return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val()));
1598  }
1599 
1601  const value_type fastAccessDx(int i) const {
1602  using std::pow; using std::log;
1603  return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val()));
1604  }
1605 
1606  protected:
1607 
1609  const ExprT2& expr2;
1610  };
1611 
1612  template <typename T1, typename T2>
1614  typename mpl::enable_if_c<
1615  ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value,
1617  >::type
1618  /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(PowerOp)*/
1619  pow (const Expr<T1>& expr1, const Expr<T2>& expr2)
1620  {
1621  typedef PowerOp< Expr<T1>, Expr<T2> > expr_t;
1622 
1623  return Expr<expr_t>(expr1, expr2);
1624  }
1625 
1626  template <typename T>
1628  Expr< PowerOp< Expr<T>, Expr<T> > >
1629  pow (const Expr<T>& expr1, const Expr<T>& expr2)
1630  {
1631  typedef PowerOp< Expr<T>, Expr<T> > expr_t;
1632 
1633  return Expr<expr_t>(expr1, expr2);
1634  }
1635 
1636  template <typename T>
1638  Expr< PowerOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
1639  pow (const typename Expr<T>::value_type& c,
1640  const Expr<T>& expr)
1641  {
1643  typedef PowerOp< ConstT, Expr<T> > expr_t;
1644 
1645  return Expr<expr_t>(ConstT(c), expr);
1646  }
1647 
1648  template <typename T>
1650  Expr< PowerOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
1651  pow (const Expr<T>& expr,
1652  const typename Expr<T>::value_type& c)
1653  {
1655  typedef PowerOp< Expr<T>, ConstT > expr_t;
1656 
1657  return Expr<expr_t>(expr, ConstT(c));
1658  }
1659 
1660  template <typename T>
1663  pow (const typename Expr<T>::scalar_type& c,
1664  const Expr<T>& expr)
1665  {
1667  typedef PowerOp< ConstT, Expr<T> > expr_t;
1668 
1669  return Expr<expr_t>(ConstT(c), expr);
1670  }
1671 
1672  template <typename T>
1675  pow (const Expr<T>& expr,
1676  const typename Expr<T>::scalar_type& c)
1677  {
1679  typedef PowerOp< Expr<T>, ConstT > expr_t;
1680 
1681  return Expr<expr_t>(expr, ConstT(c));
1682  }
1683  }
1684 }
1685 
1686 //--------------------------if_then_else operator -----------------------
1687 // Can't use the above macros because it is a ternary operator (sort of).
1688 // Also, relies on C++11
1689 
1690 #ifdef HAVE_SACADO_CXX11
1691 
1692 namespace Sacado {
1693  namespace Fad {
1694 
1695  template <typename CondT, typename ExprT1, typename ExprT2>
1696  class IfThenElseOp {};
1697 
1698  template <typename CondT, typename ExprT1, typename ExprT2>
1699  struct ExprSpec< IfThenElseOp< CondT, ExprT1, ExprT2 > > {
1700  typedef typename ExprSpec<ExprT1>::type type;
1701  };
1702 
1703  template <typename CondT, typename ExprT1, typename ExprT2>
1704  class Expr< IfThenElseOp< CondT, ExprT1, ExprT2 >,ExprSpecDefault > {
1705 
1706  public:
1707 
1708  typedef typename ExprT1::value_type value_type_1;
1709  typedef typename ExprT2::value_type value_type_2;
1710  typedef typename Sacado::Promote<value_type_1,
1711  value_type_2>::type value_type;
1712 
1713  typedef typename ExprT1::scalar_type scalar_type_1;
1714  typedef typename ExprT2::scalar_type scalar_type_2;
1715  typedef typename Sacado::Promote<scalar_type_1,
1716  scalar_type_2>::type scalar_type;
1717 
1718  typedef typename ExprT1::base_expr_type base_expr_type_1;
1719  typedef typename ExprT2::base_expr_type base_expr_type_2;
1720  typedef typename Sacado::Promote<base_expr_type_1,
1721  base_expr_type_2>::type base_expr_type;
1722 
1724  Expr(const CondT& cond_, const ExprT1& expr1_, const ExprT2& expr2_) :
1725  cond(cond_), expr1(expr1_), expr2(expr2_) {}
1726 
1728  int size() const {
1729  int sz1 = expr1.size(), sz2 = expr2.size();
1730  return sz1 > sz2 ? sz1 : sz2;
1731  }
1732 
1734  bool hasFastAccess() const {
1735  return expr1.hasFastAccess() && expr2.hasFastAccess();
1736  }
1737 
1739  bool isPassive() const {
1740  return expr1.isPassive() && expr2.isPassive();
1741  }
1742 
1744  bool updateValue() const {
1745  return expr1.updateValue() && expr2.updateValue();
1746  }
1747 
1749  void cache() const {}
1750 
1752  const value_type val() const {
1753  using Sacado::if_then_else;
1754  return if_then_else( cond, expr1.val(), expr2.val() );
1755  }
1756 
1758  const value_type dx(int i) const {
1759  using Sacado::if_then_else;
1760  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
1761  }
1762 
1764  const value_type fastAccessDx(int i) const {
1765  using Sacado::if_then_else;
1766  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
1767  }
1768 
1769  protected:
1770 
1771  const CondT& cond;
1772  const ExprT1& expr1;
1773  const ExprT2& expr2;
1774 
1775  };
1776 
1777  template <typename CondT, typename ExprT1, typename T2>
1778  struct ExprSpec< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> > > {
1779  typedef typename ExprSpec<ExprT1>::type type;
1780  };
1781 
1782  template <typename CondT, typename ExprT1, typename T2>
1783  class Expr< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
1784 
1785  public:
1786 
1787  typedef ConstExpr<T2> ConstT;
1788  typedef ConstExpr<T2> ExprT2;
1789  typedef typename ExprT1::value_type value_type_1;
1790  typedef typename ExprT2::value_type value_type_2;
1791  typedef typename Sacado::Promote<value_type_1,
1792  value_type_2>::type value_type;
1793 
1794  typedef typename ExprT1::scalar_type scalar_type_1;
1795  typedef typename ExprT2::scalar_type scalar_type_2;
1796  typedef typename Sacado::Promote<scalar_type_1,
1797  scalar_type_2>::type scalar_type;
1798 
1799  typedef typename ExprT1::base_expr_type base_expr_type_1;
1800  typedef typename ExprT2::base_expr_type base_expr_type_2;
1801  typedef typename Sacado::Promote<base_expr_type_1,
1802  base_expr_type_2>::type base_expr_type;
1803 
1805  Expr(const CondT& cond_, const ExprT1& expr1_, const ConstT& c_) :
1806  cond(cond_), expr1(expr1_), c(c_) {}
1807 
1809  int size() const {
1810  return expr1.size();
1811  }
1812 
1814  bool hasFastAccess() const {
1815  return expr1.hasFastAccess();
1816  }
1817 
1819  bool isPassive() const {
1820  return expr1.isPassive();
1821  }
1822 
1824  bool updateValue() const { return expr1.updateValue(); }
1825 
1827  void cache() const {}
1828 
1830  const value_type val() const {
1831  using Sacado::if_then_else;
1832  return if_then_else( cond, expr1.val(), c.val() );
1833  }
1834 
1836  const value_type dx(int i) const {
1837  using Sacado::if_then_else;
1838  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
1839  }
1840 
1842  const value_type fastAccessDx(int i) const {
1843  using Sacado::if_then_else;
1844  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
1845  }
1846 
1847  protected:
1848 
1849  const CondT& cond;
1850  const ExprT1& expr1;
1851  ConstT c;
1852  };
1853 
1854  template <typename CondT, typename T1, typename ExprT2>
1855  struct ExprSpec< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 > > {
1856  typedef typename ExprSpec<ExprT2>::type type;
1857  };
1858 
1859  template <typename CondT, typename T1, typename ExprT2>
1860  class Expr< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
1861 
1862  public:
1863 
1864  typedef ConstExpr<T1> ConstT;
1865  typedef ConstExpr<T1> ExprT1;
1866  typedef typename ExprT1::value_type value_type_1;
1867  typedef typename ExprT2::value_type value_type_2;
1868  typedef typename Sacado::Promote<value_type_1,
1869  value_type_2>::type value_type;
1870 
1871  typedef typename ExprT1::scalar_type scalar_type_1;
1872  typedef typename ExprT2::scalar_type scalar_type_2;
1873  typedef typename Sacado::Promote<scalar_type_1,
1874  scalar_type_2>::type scalar_type;
1875 
1876  typedef typename ExprT1::base_expr_type base_expr_type_1;
1877  typedef typename ExprT2::base_expr_type base_expr_type_2;
1878  typedef typename Sacado::Promote<base_expr_type_1,
1879  base_expr_type_2>::type base_expr_type;
1880 
1882  Expr(const CondT& cond_, const ConstT& c_, const ExprT2& expr2_) :
1883  cond(cond_), c(c_), expr2(expr2_) {}
1884 
1886  int size() const {
1887  return expr2.size();
1888  }
1889 
1891  bool hasFastAccess() const {
1892  return expr2.hasFastAccess();
1893  }
1894 
1896  bool isPassive() const {
1897  return expr2.isPassive();
1898  }
1899 
1901  bool updateValue() const { return expr2.updateValue(); }
1902 
1904  void cache() const {}
1905 
1907  const value_type val() const {
1908  using Sacado::if_then_else;
1909  return if_then_else( cond, c.val(), expr2.val() );
1910  }
1911 
1913  const value_type dx(int i) const {
1914  using Sacado::if_then_else;
1915  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
1916  }
1917 
1919  const value_type fastAccessDx(int i) const {
1920  using Sacado::if_then_else;
1921  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
1922  }
1923 
1924  protected:
1925 
1926  const CondT& cond;
1927  ConstT c;
1928  const ExprT2& expr2;
1929  };
1930 
1931  template <typename CondT, typename T1, typename T2>
1933  typename mpl::enable_if_c< IsFadExpr<T1>::value && IsFadExpr<T2>::value &&
1934  ExprLevel<T1>::value == ExprLevel<T2>::value,
1935  Expr< IfThenElseOp< CondT, T1, T2 > >
1936  >::type
1937  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
1938  {
1939  typedef IfThenElseOp< CondT, T1, T2 > expr_t;
1940 
1941  return Expr<expr_t>(cond, expr1, expr2);
1942  }
1943 
1944  template <typename CondT, typename T>
1946  Expr< IfThenElseOp< CondT, Expr<T>, Expr<T> > >
1947  if_then_else (const CondT& cond, const Expr<T>& expr1, const Expr<T>& expr2)
1948  {
1949  typedef IfThenElseOp< CondT, Expr<T>, Expr<T> > expr_t;
1950 
1951  return Expr<expr_t>(cond, expr1, expr2);
1952  }
1953 
1954  template <typename CondT, typename T>
1956  Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::value_type>,
1957  Expr<T> > >
1958  if_then_else (const CondT& cond, const typename Expr<T>::value_type& c,
1959  const Expr<T>& expr)
1960  {
1961  typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1962  typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
1963 
1964  return Expr<expr_t>(cond, ConstT(c), expr);
1965  }
1966 
1967  template <typename CondT, typename T>
1969  Expr< IfThenElseOp< CondT, Expr<T>,
1970  ConstExpr<typename Expr<T>::value_type> > >
1971  if_then_else (const CondT& cond, const Expr<T>& expr,
1972  const typename Expr<T>::value_type& c)
1973  {
1974  typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1975  typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
1976 
1977  return Expr<expr_t>(cond, expr, ConstT(c));
1978  }
1979 
1980  template <typename CondT, typename T>
1982  typename mpl::disable_if<
1983  mpl::is_same< typename Expr<T>::value_type,
1984  typename Expr<T>::scalar_type>,
1985  Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::scalar_type>,
1986  Expr<T> > >
1987  >::type
1988  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
1989  const Expr<T>& expr)
1990  {
1991  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1992  typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
1993 
1994  return Expr<expr_t>(cond, ConstT(c), expr);
1995  }
1996 
1997  template <typename CondT, typename T>
1999  typename mpl::disable_if<
2000  mpl::is_same< typename Expr<T>::value_type,
2001  typename Expr<T>::scalar_type>,
2002  Expr< IfThenElseOp< CondT, Expr<T>,
2003  ConstExpr<typename Expr<T>::scalar_type> > >
2004  >::type
2005  if_then_else (const CondT& cond, const Expr<T>& expr,
2006  const typename Expr<T>::scalar_type& c)
2007  {
2008  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
2009  typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2010 
2011  return Expr<expr_t>(cond, expr, ConstT(c));
2012  }
2013  }
2014 }
2015 
2016 #endif
2017 
2018 //-------------------------- Relational Operators -----------------------
2019 
2020 #ifdef HAVE_SACADO_CXX11
2021 
2022 namespace Sacado {
2023  namespace Fad {
2024  template <typename T1, typename T2 = T1>
2025  struct ConditionalReturnType {
2026  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
2027  };
2028  }
2029 }
2030 
2031 #define FAD_RELOP_MACRO(OP) \
2032 namespace Sacado { \
2033  namespace Fad { \
2034  template <typename ExprT1, typename ExprT2> \
2035  KOKKOS_INLINE_FUNCTION \
2036  typename ConditionalReturnType<typename Expr<ExprT1>::value_type, \
2037  typename Expr<ExprT2>::value_type>::type \
2038  operator OP (const Expr<ExprT1>& expr1, \
2039  const Expr<ExprT2>& expr2) \
2040  { \
2041  return expr1.val() OP expr2.val(); \
2042  } \
2043  \
2044  template <typename ExprT2> \
2045  KOKKOS_INLINE_FUNCTION \
2046  typename ConditionalReturnType<typename Expr<ExprT2>::value_type>::type \
2047  operator OP (const typename Expr<ExprT2>::value_type& a, \
2048  const Expr<ExprT2>& expr2) \
2049  { \
2050  return a OP expr2.val(); \
2051  } \
2052  \
2053  template <typename ExprT1> \
2054  KOKKOS_INLINE_FUNCTION \
2055  typename ConditionalReturnType<typename Expr<ExprT1>::value_type>::type \
2056  operator OP (const Expr<ExprT1>& expr1, \
2057  const typename Expr<ExprT1>::value_type& b) \
2058  { \
2059  return expr1.val() OP b; \
2060  } \
2061  } \
2062 }
2063 
2064 #else
2065 
2066 #define FAD_RELOP_MACRO(OP) \
2067 namespace Sacado { \
2068  namespace Fad { \
2069  template <typename ExprT1, typename ExprT2> \
2070  KOKKOS_INLINE_FUNCTION \
2071  bool \
2072  operator OP (const Expr<ExprT1>& expr1, \
2073  const Expr<ExprT2>& expr2) \
2074  { \
2075  return expr1.val() OP expr2.val(); \
2076  } \
2077  \
2078  template <typename ExprT2> \
2079  KOKKOS_INLINE_FUNCTION \
2080  bool \
2081  operator OP (const typename Expr<ExprT2>::value_type& a, \
2082  const Expr<ExprT2>& expr2) \
2083  { \
2084  return a OP expr2.val(); \
2085  } \
2086  \
2087  template <typename ExprT1> \
2088  KOKKOS_INLINE_FUNCTION \
2089  bool \
2090  operator OP (const Expr<ExprT1>& expr1, \
2091  const typename Expr<ExprT1>::value_type& b) \
2092  { \
2093  return expr1.val() OP b; \
2094  } \
2095  } \
2096 }
2097 
2098 #endif
2099 
2100 FAD_RELOP_MACRO(==)
2101 FAD_RELOP_MACRO(!=)
2102 FAD_RELOP_MACRO(<)
2103 FAD_RELOP_MACRO(>)
2104 FAD_RELOP_MACRO(<=)
2105 FAD_RELOP_MACRO(>=)
2106 FAD_RELOP_MACRO(<<=)
2107 FAD_RELOP_MACRO(>>=)
2108 FAD_RELOP_MACRO(&)
2109 FAD_RELOP_MACRO(|)
2110 
2111 #undef FAD_RELOP_MACRO
2112 
2113 namespace Sacado {
2114 
2115  namespace Fad {
2116 
2117  template <typename ExprT>
2119  bool operator ! (const Expr<ExprT>& expr)
2120  {
2121  return ! expr.val();
2122  }
2123 
2124  } // namespace Fad
2125 
2126 } // namespace Sacado
2127 
2128 //-------------------------- Boolean Operators -----------------------
2129 namespace Sacado {
2130 
2131  namespace Fad {
2132 
2133  template <typename ExprT>
2135  bool toBool(const Expr<ExprT>& x) {
2136  bool is_zero = (x.val() == 0.0);
2137  for (int i=0; i<x.size(); i++)
2138  is_zero = is_zero && (x.dx(i) == 0.0);
2139  return !is_zero;
2140  }
2141 
2142  } // namespace Fad
2143 
2144 } // namespace Sacado
2145 
2146 #define FAD_BOOL_MACRO(OP) \
2147 namespace Sacado { \
2148  namespace Fad { \
2149  template <typename ExprT1, typename ExprT2> \
2150  KOKKOS_INLINE_FUNCTION \
2151  bool \
2152  operator OP (const Expr<ExprT1>& expr1, \
2153  const Expr<ExprT2>& expr2) \
2154  { \
2155  return toBool(expr1) OP toBool(expr2); \
2156  } \
2157  \
2158  template <typename ExprT2> \
2159  KOKKOS_INLINE_FUNCTION \
2160  bool \
2161  operator OP (const typename Expr<ExprT2>::value_type& a, \
2162  const Expr<ExprT2>& expr2) \
2163  { \
2164  return a OP toBool(expr2); \
2165  } \
2166  \
2167  template <typename ExprT1> \
2168  KOKKOS_INLINE_FUNCTION \
2169  bool \
2170  operator OP (const Expr<ExprT1>& expr1, \
2171  const typename Expr<ExprT1>::value_type& b) \
2172  { \
2173  return toBool(expr1) OP b; \
2174  } \
2175  } \
2176 }
2177 
2178 FAD_BOOL_MACRO(&&)
2179 FAD_BOOL_MACRO(||)
2180 
2181 #undef FAD_BOOL_MACRO
2182 
2183 //-------------------------- I/O Operators -----------------------
2184 
2185 namespace Sacado {
2186 
2187  namespace Fad {
2188 
2189  template <typename ExprT>
2190  std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
2191  os << x.val() << " [";
2192 
2193  for (int i=0; i< x.size(); i++) {
2194  os << " " << x.dx(i);
2195  }
2196 
2197  os << " ]";
2198  return os;
2199  }
2200 
2201  } // namespace Fad
2202 
2203 } // namespace Sacado
2204 
2205 #endif // SACADO_FAD_OPS_HPP
Wrapper for a generic expression template.
cbrt(expr.val())
KOKKOS_INLINE_FUNCTION Expr(const ConstT &c_, const ExprT2 &expr2_)
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
expr expr SinOp
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from ConstT)
Constant expression template.
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())
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
expr expr dx(i)
abs(expr.val())
KOKKOS_INLINE_FUNCTION const value_type dx(int i) const
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MinOp
atanh(expr.val())
KOKKOS_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ConstT &c_)
expr expr CoshOp
KOKKOS_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
KOKKOS_INLINE_FUNCTION const value_type fastAccessDx(int i) const
expr expr ATanhOp
#define SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP)
expr expr TanhOp
KOKKOS_INLINE_FUNCTION const value_type fastAccessDx(int i) const
KOKKOS_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
expr expr SqrtOp
expr expr ASinhOp
expr true
atan(expr.val())
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
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)
SimpleFad< ValueT > log(const SimpleFad< ValueT > &a)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
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
Sacado::Promote< value_type_1, value_type_2 >::type value_type
#define T
Definition: Sacado_rad.hpp:573
KOKKOS_INLINE_FUNCTION T safe_sqrt(const T &x)
KOKKOS_INLINE_FUNCTION Expr(const ConstT &c_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
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)
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
#define FAD_BOOL_MACRO(OP)
expr expr ATanOp
#define T2(r, f)
Definition: Sacado_rad.hpp:578
KOKKOS_INLINE_FUNCTION const value_type dx(int i) const
SimpleFad< ValueT > min(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ACosOp
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
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
T2 base_expr_type
Typename of base-expressions.
#define T1(r, f)
Definition: Sacado_rad.hpp:603
atan2(expr1.val(), expr2.val())
Meta-function for determining nesting with an expression.
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
KOKKOS_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
sin(expr.val())
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
KOKKOS_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ConstT &c_)
KOKKOS_INLINE_FUNCTION const value_type fastAccessDx(int i) const
#define SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP)
log(expr.val())
KOKKOS_INLINE_FUNCTION const value_type fastAccessDx(int i) const
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
expr expr ACoshOp
expr expr Log10Op
expr expr SinhOp
acosh(expr.val())
SimpleFad< ValueT > operator*(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
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
T2 value_type
Typename of argument values.
KOKKOS_INLINE_FUNCTION const value_type fastAccessDx(int i) const
#define FAD_RELOP_MACRO(OP)
KOKKOS_INLINE_FUNCTION const value_type fastAccessDx(int i) const
KOKKOS_INLINE_FUNCTION bool toBool(const Expr< ExprT > &x)
KOKKOS_INLINE_FUNCTION T if_then_else(const Cond cond, const T &a, const T &b)
Sacado::Promote< value_type_1, value_type_2 >::type value_type
expr expr expr bar false
exp(expr.val())
expr expr expr ExpOp
fabs(expr.val())
expr expr AbsOp
expr expr TanOp
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
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