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 // Sacado Package
4 //
5 // Copyright 2006 NTESS and the Sacado contributors.
6 // SPDX-License-Identifier: LGPL-2.1-or-later
7 //
8 // ***********************************************************************
9 //
10 // The forward-mode AD classes in Sacado are a derivative work of the
11 // expression template classes in the Fad package by Nicolas Di Cesare.
12 // The following banner is included in the original Fad source code:
13 //
14 // ************ DO NOT REMOVE THIS BANNER ****************
15 //
16 // Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
17 // http://www.ann.jussieu.fr/~dicesare
18 //
19 // CEMRACS 98 : C++ courses,
20 // templates : new C++ techniques
21 // for scientific computing
22 //
23 //********************************************************
24 //
25 // A short implementation ( not all operators and
26 // functions are overloaded ) of 1st order Automatic
27 // Differentiation in forward mode (FAD) using
28 // EXPRESSION TEMPLATES.
29 //
30 //********************************************************
31 // @HEADER
32 
33 #ifndef SACADO_FAD_OPS_HPP
34 #define SACADO_FAD_OPS_HPP
35 
37 #include "Sacado_Fad_Ops_Fwd.hpp"
38 #include "Sacado_cmath.hpp"
39 #include <ostream> // for std::ostream
40 
41 #define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
42 namespace Sacado { \
43  namespace Fad { \
44  \
45  template <typename ExprT> \
46  class OP {}; \
47  \
48  template <typename ExprT> \
49  struct ExprSpec< OP<ExprT> > { \
50  typedef typename ExprSpec<ExprT>::type type; \
51  }; \
52  \
53  template <typename ExprT> \
54  class Expr< OP<ExprT>,ExprSpecDefault > { \
55  public: \
56  \
57  typedef typename ExprT::value_type value_type; \
58  typedef typename ExprT::scalar_type scalar_type; \
59  typedef typename ExprT::base_expr_type base_expr_type; \
60  \
61  SACADO_INLINE_FUNCTION \
62  explicit Expr(const ExprT& expr_) : expr(expr_) {} \
63  \
64  SACADO_INLINE_FUNCTION \
65  int size() const { return expr.size(); } \
66  \
67  SACADO_INLINE_FUNCTION \
68  bool hasFastAccess() const { return expr.hasFastAccess(); } \
69  \
70  SACADO_INLINE_FUNCTION \
71  bool isPassive() const { return expr.isPassive();} \
72  \
73  SACADO_INLINE_FUNCTION \
74  bool updateValue() const { return expr.updateValue(); } \
75  \
76  SACADO_INLINE_FUNCTION \
77  void cache() const {} \
78  \
79  SACADO_INLINE_FUNCTION \
80  value_type val() const { \
81  USING \
82  return VALUE; \
83  } \
84  \
85  SACADO_INLINE_FUNCTION \
86  value_type dx(int i) const { \
87  USING \
88  return DX; \
89  } \
90  \
91  SACADO_INLINE_FUNCTION \
92  value_type fastAccessDx(int i) const { \
93  USING \
94  return FASTACCESSDX; \
95  } \
96  \
97  protected: \
98  \
99  const ExprT& expr; \
100  }; \
101  \
102  template <typename T> \
103  SACADO_INLINE_FUNCTION \
104  Expr< OP< Expr<T> > > \
105  OPNAME (const Expr<T>& expr) \
106  { \
107  typedef OP< Expr<T> > expr_t; \
108  \
109  return Expr<expr_t>(expr); \
110  } \
111  } \
112  \
113 }
114 
115 FAD_UNARYOP_MACRO(operator+,
116  UnaryPlusOp,
117  ;,
118  expr.val(),
119  expr.dx(i),
120  expr.fastAccessDx(i))
121 FAD_UNARYOP_MACRO(operator-,
123  ;,
124  -expr.val(),
125  -expr.dx(i),
126  -expr.fastAccessDx(i))
129  using std::exp;,
130  exp(expr.val()),
131  exp(expr.val())*expr.dx(i),
132  exp(expr.val())*expr.fastAccessDx(i))
135  using std::log;,
136  log(expr.val()),
137  expr.dx(i)/expr.val(),
138  expr.fastAccessDx(i)/expr.val())
141  using std::log10; using std::log;,
142  log10(expr.val()),
143  expr.dx(i)/( log(value_type(10))*expr.val()),
144  expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
147  using std::sqrt;,
148  sqrt(expr.val()),
149  expr.dx(i)/(value_type(2)* sqrt(expr.val())),
150  expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
153  using std::cos; using std::sin;,
154  cos(expr.val()),
155  -expr.dx(i)* sin(expr.val()),
156  -expr.fastAccessDx(i)* sin(expr.val()))
159  using std::cos; using std::sin;,
160  sin(expr.val()),
161  expr.dx(i)* cos(expr.val()),
162  expr.fastAccessDx(i)* cos(expr.val()))
165  using std::tan;,
166  std::tan(expr.val()),
167  expr.dx(i)*
168  (value_type(1)+ tan(expr.val())* tan(expr.val())),
169  expr.fastAccessDx(i)*
170  (value_type(1)+ tan(expr.val())* tan(expr.val())))
173  using std::acos; using std::sqrt;,
174  acos(expr.val()),
175  -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
176  -expr.fastAccessDx(i) /
177  sqrt(value_type(1)-expr.val()*expr.val()))
180  using std::asin; using std::sqrt;,
181  asin(expr.val()),
182  expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
183  expr.fastAccessDx(i) /
184  sqrt(value_type(1)-expr.val()*expr.val()))
187  using std::atan;,
188  atan(expr.val()),
189  expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
190  expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
193  using std::cosh; using std::sinh;,
194  cosh(expr.val()),
195  expr.dx(i)* sinh(expr.val()),
196  expr.fastAccessDx(i)* sinh(expr.val()))
197 FAD_UNARYOP_MACRO(sinh,
199  using std::cosh; using std::sinh;,
200  sinh(expr.val()),
201  expr.dx(i)* cosh(expr.val()),
202  expr.fastAccessDx(i)* cosh(expr.val()))
205  using std::tanh;,
206  tanh(expr.val()),
207  expr.dx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())),
208  expr.fastAccessDx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())))
211  using std::acosh; using std::sqrt;,
212  acosh(expr.val()),
213  expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
214  (expr.val()+value_type(1))),
215  expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
216  (expr.val()+value_type(1))))
219  using std::asinh; using std::sqrt;,
220  asinh(expr.val()),
221  expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
222  expr.fastAccessDx(i)/ sqrt(value_type(1)+
223  expr.val()*expr.val()))
226  using std::atanh;,
227  atanh(expr.val()),
228  expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
229  expr.fastAccessDx(i)/(value_type(1)-
230  expr.val()*expr.val()))
233  using std::abs; using Sacado::if_then_else;,
234  abs(expr.val()),
235  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
236  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
239  using std::fabs; using Sacado::if_then_else;,
240  fabs(expr.val()),
241  if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
242  if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
245  using std::cbrt;,
246  cbrt(expr.val()),
247  expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
248  expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
249 
250 #undef FAD_UNARYOP_MACRO
251 
252 // Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
253 // "simd" value types that use if_then_else(). The only reason for not using
254 // if_then_else() always is to avoid evaluating the derivative if the value is
255 // zero to avoid throwing FPEs.
256 namespace Sacado {
257  namespace Fad {
258 
259  template <typename ExprT, bool is_simd>
260  class SafeSqrtOp {};
261 
262  template <typename ExprT>
263  struct ExprSpec< SafeSqrtOp<ExprT> > {
264  typedef typename ExprSpec<ExprT>::type type;
265  };
266 
267  //
268  // Implementation for simd type using if_then_else()
269  //
270  template <typename ExprT>
271  class Expr< SafeSqrtOp<ExprT,true>,ExprSpecDefault > {
272  public:
273 
274  typedef typename ExprT::value_type value_type;
275  typedef typename ExprT::scalar_type scalar_type;
276  typedef typename ExprT::base_expr_type base_expr_type;
277 
279  explicit Expr(const ExprT& expr_) : expr(expr_) {}
280 
282  int size() const { return expr.size(); }
283 
285  bool hasFastAccess() const { return expr.hasFastAccess(); }
286 
288  bool isPassive() const { return expr.isPassive();}
289 
291  bool updateValue() const { return expr.updateValue(); }
292 
294  void cache() const {}
295 
297  value_type val() const {
298  using std::sqrt;
299  return sqrt(expr.val());
300  }
301 
303  value_type dx(int i) const {
304  using std::sqrt; using Sacado::if_then_else;
305  return if_then_else(
306  expr.val() == value_type(0.0), value_type(0.0),
307  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
308  }
309 
311  value_type fastAccessDx(int i) const {
312  using std::sqrt; using Sacado::if_then_else;
313  return if_then_else(
314  expr.val() == value_type(0.0), value_type(0.0),
315  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
316  }
317 
318  protected:
319 
320  const ExprT& expr;
321  };
322 
323  //
324  // Specialization for scalar types using ternary operator
325  //
326  template <typename ExprT>
327  class Expr< SafeSqrtOp<ExprT,false>,ExprSpecDefault > {
328  public:
329 
330  typedef typename ExprT::value_type value_type;
331  typedef typename ExprT::scalar_type scalar_type;
332  typedef typename ExprT::base_expr_type base_expr_type;
333 
335  explicit Expr(const ExprT& expr_) : expr(expr_) {}
336 
338  int size() const { return expr.size(); }
339 
341  bool hasFastAccess() const { return expr.hasFastAccess(); }
342 
344  bool isPassive() const { return expr.isPassive();}
345 
347  bool updateValue() const { return expr.updateValue(); }
348 
350  void cache() const {}
351 
353  value_type val() const {
354  using std::sqrt;
355  return sqrt(expr.val());
356  }
357 
359  value_type dx(int i) const {
360  using std::sqrt;
361  return expr.val() == value_type(0.0) ? value_type(0.0) :
362  value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
363  }
364 
366  value_type fastAccessDx(int i) const {
367  using std::sqrt;
368  return expr.val() == value_type(0.0) ? value_type(0.0) :
369  value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
370  }
371 
372  protected:
373 
374  const ExprT& expr;
375  };
376 
377  template <typename T>
379  Expr< SafeSqrtOp< Expr<T> > >
380  safe_sqrt (const Expr<T>& expr)
381  {
382  typedef SafeSqrtOp< Expr<T> > expr_t;
383 
384  return Expr<expr_t>(expr);
385  }
386  }
387 
388 }
389 
390 #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) \
391 namespace Sacado { \
392  namespace Fad { \
393  \
394  template <typename ExprT1, typename ExprT2> \
395  class OP {}; \
396  \
397  template <typename ExprT1, typename ExprT2> \
398  struct ExprSpec< OP< ExprT1, ExprT2 > > { \
399  typedef typename ExprSpec<ExprT1>::type type; \
400  }; \
401  \
402  template <typename ExprT1, typename ExprT2> \
403  class Expr< OP< ExprT1, ExprT2 >,ExprSpecDefault > { \
404  \
405  public: \
406  \
407  typedef typename ExprT1::value_type value_type_1; \
408  typedef typename ExprT2::value_type value_type_2; \
409  typedef typename Sacado::Promote<value_type_1, \
410  value_type_2>::type value_type; \
411  \
412  typedef typename ExprT1::scalar_type scalar_type_1; \
413  typedef typename ExprT2::scalar_type scalar_type_2; \
414  typedef typename Sacado::Promote<scalar_type_1, \
415  scalar_type_2>::type scalar_type; \
416  \
417  typedef typename ExprT1::base_expr_type base_expr_type_1; \
418  typedef typename ExprT2::base_expr_type base_expr_type_2; \
419  typedef typename Sacado::Promote<base_expr_type_1, \
420  base_expr_type_2>::type base_expr_type; \
421  \
422  SACADO_INLINE_FUNCTION \
423  Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
424  expr1(expr1_), expr2(expr2_) {} \
425  \
426  SACADO_INLINE_FUNCTION \
427  int size() const { \
428  int sz1 = expr1.size(), sz2 = expr2.size(); \
429  return sz1 > sz2 ? sz1 : sz2; \
430  } \
431  \
432  SACADO_INLINE_FUNCTION \
433  bool hasFastAccess() const { \
434  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
435  } \
436  \
437  SACADO_INLINE_FUNCTION \
438  bool isPassive() const { \
439  return expr1.isPassive() && expr2.isPassive(); \
440  } \
441  \
442  SACADO_INLINE_FUNCTION \
443  bool updateValue() const { \
444  return expr1.updateValue() && expr2.updateValue(); \
445  } \
446  \
447  SACADO_INLINE_FUNCTION \
448  void cache() const {} \
449  \
450  SACADO_INLINE_FUNCTION \
451  const value_type val() const { \
452  USING \
453  return VALUE; \
454  } \
455  \
456  SACADO_INLINE_FUNCTION \
457  const value_type dx(int i) const { \
458  USING \
459  return DX; \
460  } \
461  \
462  SACADO_INLINE_FUNCTION \
463  const value_type fastAccessDx(int i) const { \
464  USING \
465  return FASTACCESSDX; \
466  } \
467  \
468  protected: \
469  \
470  const ExprT1& expr1; \
471  const ExprT2& expr2; \
472  \
473  }; \
474  \
475  template <typename ExprT1, typename T2> \
476  struct ExprSpec< OP< ExprT1, ConstExpr<T2> > > { \
477  typedef typename ExprSpec<ExprT1>::type type; \
478  }; \
479  \
480  template <typename ExprT1, typename T2> \
481  class Expr< OP< ExprT1, ConstExpr<T2> >,ExprSpecDefault > { \
482  \
483  public: \
484  \
485  typedef ConstExpr<T2> ConstT; \
486  typedef ConstExpr<T2> ExprT2; \
487  typedef typename ExprT1::value_type value_type_1; \
488  typedef typename ExprT2::value_type value_type_2; \
489  typedef typename Sacado::Promote<value_type_1, \
490  value_type_2>::type value_type; \
491  \
492  typedef typename ExprT1::scalar_type scalar_type_1; \
493  typedef typename ExprT2::scalar_type scalar_type_2; \
494  typedef typename Sacado::Promote<scalar_type_1, \
495  scalar_type_2>::type scalar_type; \
496  \
497  typedef typename ExprT1::base_expr_type base_expr_type_1; \
498  typedef typename ExprT2::base_expr_type base_expr_type_2; \
499  typedef typename Sacado::Promote<base_expr_type_1, \
500  base_expr_type_2>::type base_expr_type; \
501  \
502  SACADO_INLINE_FUNCTION \
503  Expr(const ExprT1& expr1_, const ConstT& c_) : \
504  expr1(expr1_), c(c_) {} \
505  \
506  SACADO_INLINE_FUNCTION \
507  int size() const { \
508  return expr1.size(); \
509  } \
510  \
511  SACADO_INLINE_FUNCTION \
512  bool hasFastAccess() const { \
513  return expr1.hasFastAccess(); \
514  } \
515  \
516  SACADO_INLINE_FUNCTION \
517  bool isPassive() const { \
518  return expr1.isPassive(); \
519  } \
520  \
521  SACADO_INLINE_FUNCTION \
522  bool updateValue() const { return expr1.updateValue(); } \
523  \
524  SACADO_INLINE_FUNCTION \
525  void cache() const {} \
526  \
527  SACADO_INLINE_FUNCTION \
528  const value_type val() const { \
529  USING \
530  return VAL_CONST_DX_2; \
531  } \
532  \
533  SACADO_INLINE_FUNCTION \
534  const value_type dx(int i) const { \
535  USING \
536  return CONST_DX_2; \
537  } \
538  \
539  SACADO_INLINE_FUNCTION \
540  const value_type fastAccessDx(int i) const { \
541  USING \
542  return CONST_FASTACCESSDX_2; \
543  } \
544  \
545  protected: \
546  \
547  const ExprT1& expr1; \
548  ConstT c; \
549  }; \
550  \
551  template <typename T1, typename ExprT2> \
552  struct ExprSpec< OP< ConstExpr<T1>, ExprT2 > > { \
553  typedef typename ExprSpec<ExprT2>::type type; \
554  }; \
555  \
556  template <typename T1, typename ExprT2> \
557  class Expr< OP< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > { \
558  \
559  public: \
560  \
561  typedef ConstExpr<T1> ConstT; \
562  typedef ConstExpr<T1> ExprT1; \
563  typedef typename ExprT1::value_type value_type_1; \
564  typedef typename ExprT2::value_type value_type_2; \
565  typedef typename Sacado::Promote<value_type_1, \
566  value_type_2>::type value_type; \
567  \
568  typedef typename ExprT1::scalar_type scalar_type_1; \
569  typedef typename ExprT2::scalar_type scalar_type_2; \
570  typedef typename Sacado::Promote<scalar_type_1, \
571  scalar_type_2>::type scalar_type; \
572  \
573  typedef typename ExprT1::base_expr_type base_expr_type_1; \
574  typedef typename ExprT2::base_expr_type base_expr_type_2; \
575  typedef typename Sacado::Promote<base_expr_type_1, \
576  base_expr_type_2>::type base_expr_type; \
577  \
578  \
579  SACADO_INLINE_FUNCTION \
580  Expr(const ConstT& c_, const ExprT2& expr2_) : \
581  c(c_), expr2(expr2_) {} \
582  \
583  SACADO_INLINE_FUNCTION \
584  int size() const { \
585  return expr2.size(); \
586  } \
587  \
588  SACADO_INLINE_FUNCTION \
589  bool hasFastAccess() const { \
590  return expr2.hasFastAccess(); \
591  } \
592  \
593  SACADO_INLINE_FUNCTION \
594  bool isPassive() const { \
595  return expr2.isPassive(); \
596  } \
597  \
598  SACADO_INLINE_FUNCTION \
599  bool updateValue() const { return expr2.updateValue(); } \
600  \
601  SACADO_INLINE_FUNCTION \
602  void cache() const {} \
603  \
604  SACADO_INLINE_FUNCTION \
605  const value_type val() const { \
606  USING \
607  return VAL_CONST_DX_1; \
608  } \
609  \
610  SACADO_INLINE_FUNCTION \
611  const value_type dx(int i) const { \
612  USING \
613  return CONST_DX_1; \
614  } \
615  \
616  SACADO_INLINE_FUNCTION \
617  const value_type fastAccessDx(int i) const { \
618  USING \
619  return CONST_FASTACCESSDX_1; \
620  } \
621  \
622  protected: \
623  \
624  ConstT c; \
625  const ExprT2& expr2; \
626  }; \
627  \
628  template <typename T1, typename T2> \
629  SACADO_INLINE_FUNCTION \
630  typename mpl::enable_if_c< \
631  ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value, \
632  Expr< OP< Expr<T1>, Expr<T2> > > \
633  >::type \
634  /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(OP)*/ \
635  OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2) \
636  { \
637  typedef OP< Expr<T1>, Expr<T2> > expr_t; \
638  \
639  return Expr<expr_t>(expr1, expr2); \
640  } \
641  \
642  template <typename T> \
643  SACADO_INLINE_FUNCTION \
644  Expr< OP< Expr<T>, Expr<T> > > \
645  OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
646  { \
647  typedef OP< Expr<T>, Expr<T> > expr_t; \
648  \
649  return Expr<expr_t>(expr1, expr2); \
650  } \
651  \
652  template <typename T> \
653  SACADO_INLINE_FUNCTION \
654  Expr< OP< ConstExpr<typename Expr<T>::value_type>, \
655  Expr<T> > > \
656  OPNAME (const typename Expr<T>::value_type& c, \
657  const Expr<T>& expr) \
658  { \
659  typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
660  typedef OP< ConstT, Expr<T> > expr_t; \
661  \
662  return Expr<expr_t>(ConstT(c), expr); \
663  } \
664  \
665  template <typename T> \
666  SACADO_INLINE_FUNCTION \
667  Expr< OP< Expr<T>, \
668  ConstExpr<typename Expr<T>::value_type> > > \
669  OPNAME (const Expr<T>& expr, \
670  const typename Expr<T>::value_type& c) \
671  { \
672  typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
673  typedef OP< Expr<T>, ConstT > expr_t; \
674  \
675  return Expr<expr_t>(expr, ConstT(c)); \
676  } \
677  \
678  template <typename T> \
679  SACADO_INLINE_FUNCTION \
680  SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP) \
681  OPNAME (const typename Expr<T>::scalar_type& c, \
682  const Expr<T>& expr) \
683  { \
684  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
685  typedef OP< ConstT, Expr<T> > expr_t; \
686  \
687  return Expr<expr_t>(ConstT(c), expr); \
688  } \
689  \
690  template <typename T> \
691  SACADO_INLINE_FUNCTION \
692  SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP) \
693  OPNAME (const Expr<T>& expr, \
694  const typename Expr<T>::scalar_type& c) \
695  { \
696  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
697  typedef OP< Expr<T>, ConstT > expr_t; \
698  \
699  return Expr<expr_t>(expr, ConstT(c)); \
700  } \
701  } \
702  \
703 }
704 
705 
706 FAD_BINARYOP_MACRO(operator+,
707  AdditionOp,
708  ;,
709  expr1.val() + expr2.val(),
710  expr1.dx(i) + expr2.dx(i),
711  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
712  c.val() + expr2.val(),
713  expr1.val() + c.val(),
714  expr2.dx(i),
715  expr1.dx(i),
716  expr2.fastAccessDx(i),
717  expr1.fastAccessDx(i))
718 FAD_BINARYOP_MACRO(operator-,
720  ;,
721  expr1.val() - expr2.val(),
722  expr1.dx(i) - expr2.dx(i),
723  expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
724  c.val() - expr2.val(),
725  expr1.val() - c.val(),
726  -expr2.dx(i),
727  expr1.dx(i),
728  -expr2.fastAccessDx(i),
729  expr1.fastAccessDx(i))
730 // FAD_BINARYOP_MACRO(operator*,
731 // MultiplicationOp,
732 // ;,
733 // expr1.val() * expr2.val(),
734 // expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
735 // expr1.val()*expr2.fastAccessDx(i) +
736 // expr1.fastAccessDx(i)*expr2.val(),
737 // c.val() * expr2.val(),
738 // expr1.val() * c.val(),
739 // c.val()*expr2.dx(i),
740 // expr1.dx(i)*c.val(),
741 // c.val()*expr2.fastAccessDx(i),
742 // expr1.fastAccessDx(i)*c.val())
743 FAD_BINARYOP_MACRO(operator/,
745  ;,
746  expr1.val() / expr2.val(),
747  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
748  (expr2.val()*expr2.val()),
749  (expr1.fastAccessDx(i)*expr2.val() -
750  expr2.fastAccessDx(i)*expr1.val()) /
751  (expr2.val()*expr2.val()),
752  c.val() / expr2.val(),
753  expr1.val() / c.val(),
754  -expr2.dx(i)*c.val() / (expr2.val()*expr2.val()),
755  expr1.dx(i)/c.val(),
756  -expr2.fastAccessDx(i)*c.val() / (expr2.val()*expr2.val()),
757  expr1.fastAccessDx(i)/c.val())
760  using std::atan2;,
761  atan2(expr1.val(), expr2.val()),
762  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
763  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
764  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
765  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
766  atan2(c.val(), expr2.val()),
767  atan2(expr1.val(), c.val()),
768  (-c.val()*expr2.dx(i)) / (c.val()*c.val() + expr2.val()*expr2.val()),
769  (c.val()*expr1.dx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()),
770  (-c.val()*expr2.fastAccessDx(i))/ (c.val()*c.val() + expr2.val()*expr2.val()),
771  (c.val()*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()))
772 // FAD_BINARYOP_MACRO(pow,
773 // PowerOp,
774 // using std::pow; using std::log; using Sacado::if_then_else;,
775 // pow(expr1.val(), expr2.val()),
776 // 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())) ),
777 // 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())) ),
778 // pow(c.val(), expr2.val()),
779 // pow(expr1.val(), c.val()),
780 // 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())) ),
781 // 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())) ),
782 // 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())) ),
783 // 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()))) )
786  using Sacado::if_then_else;,
787  if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
788  if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
789  if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
790  if_then_else( c.val() >= expr2.val(), value_type(c.val()), expr2.val() ),
791  if_then_else( expr1.val() >= c.val(), expr1.val(), value_type(c.val()) ),
792  if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
793  if_then_else( expr1.val() >= c.val(), expr1.dx(i), value_type(0.0) ),
794  if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
795  if_then_else( expr1.val() >= c.val(), expr1.fastAccessDx(i), value_type(0.0) ) )
798  using Sacado::if_then_else;,
799  if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
800  if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
801  if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
802  if_then_else( c.val() <= expr2.val(), value_type(c.val()), expr2.val() ),
803  if_then_else( expr1.val() <= c.val(), expr1.val(), value_type(c.val()) ),
804  if_then_else( c.val() <= expr2.val(), value_type(0), expr2.dx(i) ),
805  if_then_else( expr1.val() <= c.val(), expr1.dx(i), value_type(0) ),
806  if_then_else( c.val() <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
807  if_then_else( expr1.val() <= c.val(), expr1.fastAccessDx(i), value_type(0) ) )
808 
809 
810 #undef FAD_BINARYOP_MACRO
811 
812 namespace Sacado {
813  namespace Fad {
814 
815  template <typename ExprT1, typename ExprT2>
816  class MultiplicationOp {};
817 
818  template <typename ExprT1, typename ExprT2>
819  struct ExprSpec< MultiplicationOp< ExprT1, ExprT2 > > {
820  typedef typename ExprSpec<ExprT1>::type type;
821  };
822 
823  template <typename ExprT1, typename ExprT2>
824  class Expr< MultiplicationOp< ExprT1, ExprT2 >,ExprSpecDefault > {
825 
826  public:
827 
828  typedef typename ExprT1::value_type value_type_1;
829  typedef typename ExprT2::value_type value_type_2;
830  typedef typename Sacado::Promote<value_type_1,
831  value_type_2>::type value_type;
832 
833  typedef typename ExprT1::scalar_type scalar_type_1;
834  typedef typename ExprT2::scalar_type scalar_type_2;
835  typedef typename Sacado::Promote<scalar_type_1,
836  scalar_type_2>::type scalar_type;
837 
838  typedef typename ExprT1::base_expr_type base_expr_type_1;
839  typedef typename ExprT2::base_expr_type base_expr_type_2;
840  typedef typename Sacado::Promote<base_expr_type_1,
841  base_expr_type_2>::type base_expr_type;
842 
844  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
845  expr1(expr1_), expr2(expr2_) {}
846 
848  int size() const {
849  int sz1 = expr1.size(), sz2 = expr2.size();
850  return sz1 > sz2 ? sz1 : sz2;
851  }
852 
854  bool hasFastAccess() const {
855  return expr1.hasFastAccess() && expr2.hasFastAccess();
856  }
857 
859  bool isPassive() const {
860  return expr1.isPassive() && expr2.isPassive();
861  }
862 
864  bool updateValue() const {
865  return expr1.updateValue() && expr2.updateValue();
866  }
867 
869  void cache() const {}
870 
872  const value_type val() const {
873  return expr1.val()*expr2.val();
874  }
875 
877  const value_type dx(int i) const {
878  if (expr1.size() > 0 && expr2.size() > 0)
879  return expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val();
880  else if (expr1.size() > 0)
881  return expr1.dx(i)*expr2.val();
882  else
883  return expr1.val()*expr2.dx(i);
884  }
885 
887  const value_type fastAccessDx(int i) const {
888  return expr1.val()*expr2.fastAccessDx(i) +
889  expr1.fastAccessDx(i)*expr2.val();
890  }
891 
892  protected:
893 
894  const ExprT1& expr1;
895  const ExprT2& expr2;
896 
897  };
898 
899  template <typename ExprT1, typename T2>
900  struct ExprSpec< MultiplicationOp< ExprT1, ConstExpr<T2> > > {
901  typedef typename ExprSpec<ExprT1>::type type;
902  };
903 
904  template <typename ExprT1, typename T2>
905  class Expr< MultiplicationOp< ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
906 
907  public:
908 
909  typedef ConstExpr<T2> ConstT;
910  typedef ConstExpr<T2> ExprT2;
911  typedef typename ExprT1::value_type value_type_1;
912  typedef typename ExprT2::value_type value_type_2;
913  typedef typename Sacado::Promote<value_type_1,
914  value_type_2>::type value_type;
915 
916  typedef typename ExprT1::scalar_type scalar_type_1;
917  typedef typename ExprT2::scalar_type scalar_type_2;
918  typedef typename Sacado::Promote<scalar_type_1,
919  scalar_type_2>::type scalar_type;
920 
921  typedef typename ExprT1::base_expr_type base_expr_type_1;
922  typedef typename ExprT2::base_expr_type base_expr_type_2;
923  typedef typename Sacado::Promote<base_expr_type_1,
924  base_expr_type_2>::type base_expr_type;
925 
927  Expr(const ExprT1& expr1_, const ConstT& c_) :
928  expr1(expr1_), c(c_) {}
929 
931  int size() const {
932  return expr1.size();
933  }
934 
936  bool hasFastAccess() const {
937  return expr1.hasFastAccess();
938  }
939 
941  bool isPassive() const {
942  return expr1.isPassive();
943  }
944 
946  bool updateValue() const { return expr1.updateValue(); }
947 
949  void cache() const {}
950 
952  const value_type val() const {
953  return expr1.val()*c.val();
954  }
955 
957  const value_type dx(int i) const {
958  return expr1.dx(i)*c.val();
959  }
960 
962  const value_type fastAccessDx(int i) const {
963  return expr1.fastAccessDx(i)*c.val();
964  }
965 
966  protected:
967 
968  const ExprT1& expr1;
969  ConstT c;
970  };
971 
972  template <typename T1, typename ExprT2>
973  struct ExprSpec< MultiplicationOp< ConstExpr<T1>, ExprT2 > > {
974  typedef typename ExprSpec<ExprT2>::type type;
975  };
976 
977  template <typename T1, typename ExprT2>
978  class Expr< MultiplicationOp< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
979 
980  public:
981 
982  typedef ConstExpr<T1> ConstT;
983  typedef ConstExpr<T1> ExprT1;
984  typedef typename ExprT1::value_type value_type_1;
985  typedef typename ExprT2::value_type value_type_2;
986  typedef typename Sacado::Promote<value_type_1,
987  value_type_2>::type value_type;
988 
989  typedef typename ExprT1::scalar_type scalar_type_1;
990  typedef typename ExprT2::scalar_type scalar_type_2;
991  typedef typename Sacado::Promote<scalar_type_1,
992  scalar_type_2>::type scalar_type;
993 
994  typedef typename ExprT1::base_expr_type base_expr_type_1;
995  typedef typename ExprT2::base_expr_type base_expr_type_2;
996  typedef typename Sacado::Promote<base_expr_type_1,
997  base_expr_type_2>::type base_expr_type;
998 
1000  Expr(const ConstT& c_, const ExprT2& expr2_) :
1001  c(c_), expr2(expr2_) {}
1002 
1004  int size() const {
1005  return expr2.size();
1006  }
1007 
1009  bool hasFastAccess() const {
1010  return expr2.hasFastAccess();
1011  }
1012 
1014  bool isPassive() const {
1015  return expr2.isPassive();
1016  }
1017 
1019  bool updateValue() const { return expr2.updateValue(); }
1020 
1022  void cache() const {}
1023 
1025  const value_type val() const {
1026  return c.val()*expr2.val();
1027  }
1028 
1030  const value_type dx(int i) const {
1031  return c.val()*expr2.dx(i);
1032  }
1033 
1035  const value_type fastAccessDx(int i) const {
1036  return c.val()*expr2.fastAccessDx(i);
1037  }
1038 
1039  protected:
1040 
1041  ConstT c;
1042  const ExprT2& expr2;
1043  };
1044 
1045  template <typename T1, typename T2>
1047  typename mpl::enable_if_c<
1048  ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value,
1049  Expr< MultiplicationOp< Expr<T1>, Expr<T2> > >
1050  >::type
1051  /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(MultiplicationOp)*/
1052  operator* (const Expr<T1>& expr1, const Expr<T2>& expr2)
1053  {
1054  typedef MultiplicationOp< Expr<T1>, Expr<T2> > expr_t;
1055 
1056  return Expr<expr_t>(expr1, expr2);
1057  }
1058 
1059  template <typename T>
1061  Expr< MultiplicationOp< Expr<T>, Expr<T> > >
1062  operator* (const Expr<T>& expr1, const Expr<T>& expr2)
1063  {
1064  typedef MultiplicationOp< Expr<T>, Expr<T> > expr_t;
1065 
1066  return Expr<expr_t>(expr1, expr2);
1067  }
1068 
1069  template <typename T>
1071  Expr< MultiplicationOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
1072  operator* (const typename Expr<T>::value_type& c,
1073  const Expr<T>& expr)
1074  {
1075  typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1076  typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1077 
1078  return Expr<expr_t>(ConstT(c), expr);
1079  }
1080 
1081  template <typename T>
1083  Expr< MultiplicationOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
1084  operator* (const Expr<T>& expr,
1085  const typename Expr<T>::value_type& c)
1086  {
1087  typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1088  typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1089 
1090  return Expr<expr_t>(expr, ConstT(c));
1091  }
1092 
1093  template <typename T>
1095  SACADO_FAD_OP_ENABLE_SCALAR_EXPR(MultiplicationOp)
1096  operator* (const typename Expr<T>::scalar_type& c,
1097  const Expr<T>& expr)
1098  {
1099  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1100  typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1101 
1102  return Expr<expr_t>(ConstT(c), expr);
1103  }
1104 
1105  template <typename T>
1107  SACADO_FAD_OP_ENABLE_EXPR_SCALAR(MultiplicationOp)
1108  operator* (const Expr<T>& expr,
1109  const typename Expr<T>::scalar_type& c)
1110  {
1111  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1112  typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1113 
1114  return Expr<expr_t>(expr, ConstT(c));
1115  }
1116  }
1117 }
1118 
1119 // Special handling for std::pow() to provide specializations of PowerOp for
1120 // "simd" value types that use if_then_else(). The only reason for not using
1121 // if_then_else() always is to avoid evaluating the derivative if the value is
1122 // zero to avoid throwing FPEs.
1123 namespace Sacado {
1124  namespace Fad {
1125 
1126  template <typename ExprT1, typename ExprT2, typename Impl>
1127  class PowerOp {};
1128 
1129  template <typename ExprT1, typename ExprT2>
1130  struct ExprSpec< PowerOp< ExprT1, ExprT2 > > {
1131  typedef typename ExprSpec<ExprT1>::type type;
1132  };
1133 
1134  template <typename ExprT1, typename T2>
1135  struct ExprSpec<PowerOp< ExprT1, ConstExpr<T2> > > {
1136  typedef typename ExprSpec<ExprT1>::type type;
1137  };
1138 
1139  template <typename T1, typename ExprT2>
1140  struct ExprSpec< PowerOp< ConstExpr<T1>, ExprT2 > > {
1141  typedef typename ExprSpec<ExprT2>::type type;
1142  };
1143 
1144  //
1145  // Implementation for simd type using if_then_else()
1146  //
1147  template <typename ExprT1, typename ExprT2>
1148  class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Simd >, ExprSpecDefault > {
1149 
1150  public:
1151 
1152  typedef typename ExprT1::value_type value_type_1;
1153  typedef typename ExprT2::value_type value_type_2;
1154  typedef typename Sacado::Promote<value_type_1,
1156 
1157  typedef typename ExprT1::scalar_type scalar_type_1;
1158  typedef typename ExprT2::scalar_type scalar_type_2;
1159  typedef typename Sacado::Promote<scalar_type_1,
1161 
1162  typedef typename ExprT1::base_expr_type base_expr_type_1;
1163  typedef typename ExprT2::base_expr_type base_expr_type_2;
1164  typedef typename Sacado::Promote<base_expr_type_1,
1166 
1168  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1169  expr1(expr1_), expr2(expr2_) {}
1170 
1172  int size() const {
1173  int sz1 = expr1.size(), sz2 = expr2.size();
1174  return sz1 > sz2 ? sz1 : sz2;
1175  }
1176 
1178  bool hasFastAccess() const {
1179  return expr1.hasFastAccess() && expr2.hasFastAccess();
1180  }
1181 
1183  bool isPassive() const {
1184  return expr1.isPassive() && expr2.isPassive();
1185  }
1186 
1188  bool updateValue() const {
1189  return expr1.updateValue() && expr2.updateValue();
1190  }
1191 
1193  void cache() const {}
1194 
1196  const value_type val() const {
1197  using std::pow;
1198  return pow(expr1.val(), expr2.val());
1199  }
1200 
1202  const value_type dx(int i) const {
1203  using std::pow; using std::log; using Sacado::if_then_else;
1204  const int sz1 = expr1.size(), sz2 = expr2.size();
1205  if (sz1 > 0 && sz2 > 0)
1206  return if_then_else( expr1.val() == scalar_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())) );
1207  else if (sz1 > 0)
1208  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1209  // It seems less accurate and caused convergence problems in some codes
1210  return if_then_else( expr2.val() == scalar_type(1.0), expr1.dx(i), if_then_else(expr1.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val())) ));
1211  else
1212  return if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val())) );
1213  }
1214 
1216  const value_type fastAccessDx(int i) const {
1217  using std::pow; using std::log; using Sacado::if_then_else;
1218  return if_then_else( expr1.val() == scalar_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())) );
1219  }
1220 
1221  protected:
1222 
1223  const ExprT1& expr1;
1224  const ExprT2& expr2;
1225 
1226  };
1227 
1228  template <typename ExprT1, typename T2>
1229  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Simd >,
1230  ExprSpecDefault > {
1231 
1232  public:
1233 
1236  typedef typename ExprT1::value_type value_type_1;
1238  typedef typename Sacado::Promote<value_type_1,
1240 
1241  typedef typename ExprT1::scalar_type scalar_type_1;
1243  typedef typename Sacado::Promote<scalar_type_1,
1245 
1246  typedef typename ExprT1::base_expr_type base_expr_type_1;
1248  typedef typename Sacado::Promote<base_expr_type_1,
1250 
1252  Expr(const ExprT1& expr1_, const ConstT& c_) :
1253  expr1(expr1_), c(c_) {}
1254 
1256  int size() const {
1257  return expr1.size();
1258  }
1259 
1261  bool hasFastAccess() const {
1262  return expr1.hasFastAccess();
1263  }
1264 
1266  bool isPassive() const {
1267  return expr1.isPassive();
1268  }
1269 
1271  bool updateValue() const { return expr1.updateValue(); }
1272 
1274  void cache() const {}
1275 
1277  const value_type val() const {
1278  using std::pow;
1279  return pow(expr1.val(), c.val());
1280  }
1281 
1283  const value_type dx(int i) const {
1284  using std::pow; using Sacado::if_then_else;
1285  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1286  // It seems less accurate and caused convergence problems in some codes
1287  return if_then_else( c.val() == scalar_type(1.0), expr1.dx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val())) ));
1288  }
1289 
1291  const value_type fastAccessDx(int i) const {
1292  using std::pow; using Sacado::if_then_else;
1293  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1294  // It seems less accurate and caused convergence problems in some codes
1295  return if_then_else( c.val() == scalar_type(1.0), expr1.fastAccessDx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()))));
1296  }
1297 
1298  protected:
1299 
1300  const ExprT1& expr1;
1302  };
1303 
1304  template <typename T1, typename ExprT2>
1305  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Simd >,
1306  ExprSpecDefault > {
1307 
1308  public:
1309 
1313  typedef typename ExprT2::value_type value_type_2;
1314  typedef typename Sacado::Promote<value_type_1,
1316 
1318  typedef typename ExprT2::scalar_type scalar_type_2;
1319  typedef typename Sacado::Promote<scalar_type_1,
1321 
1323  typedef typename ExprT2::base_expr_type base_expr_type_2;
1324  typedef typename Sacado::Promote<base_expr_type_1,
1326 
1327 
1329  Expr(const ConstT& c_, const ExprT2& expr2_) :
1330  c(c_), expr2(expr2_) {}
1331 
1333  int size() const {
1334  return expr2.size();
1335  }
1336 
1338  bool hasFastAccess() const {
1339  return expr2.hasFastAccess();
1340  }
1341 
1343  bool isPassive() const {
1344  return expr2.isPassive();
1345  }
1346 
1348  bool updateValue() const { return expr2.updateValue(); }
1349 
1351  void cache() const {}
1352 
1354  const value_type val() const {
1355  using std::pow;
1356  return pow(c.val(), expr2.val());
1357  }
1358 
1360  const value_type dx(int i) const {
1361  using std::pow; using std::log; using Sacado::if_then_else;
1362  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())) );
1363  }
1364 
1366  const value_type fastAccessDx(int i) const {
1367  using std::pow; using std::log; using Sacado::if_then_else;
1368  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())) );
1369  }
1370 
1371  protected:
1372 
1374  const ExprT2& expr2;
1375  };
1376 
1377  //
1378  // Specialization for scalar types using ternary operator
1379  //
1380 
1381  template <typename ExprT1, typename ExprT2>
1382  class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Scalar >,
1383  ExprSpecDefault > {
1384 
1385  public:
1386 
1387  typedef typename ExprT1::value_type value_type_1;
1388  typedef typename ExprT2::value_type value_type_2;
1389  typedef typename Sacado::Promote<value_type_1,
1391 
1392  typedef typename ExprT1::scalar_type scalar_type_1;
1393  typedef typename ExprT2::scalar_type scalar_type_2;
1394  typedef typename Sacado::Promote<scalar_type_1,
1396 
1397  typedef typename ExprT1::base_expr_type base_expr_type_1;
1398  typedef typename ExprT2::base_expr_type base_expr_type_2;
1399  typedef typename Sacado::Promote<base_expr_type_1,
1401 
1403  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1404  expr1(expr1_), expr2(expr2_) {}
1405 
1407  int size() const {
1408  int sz1 = expr1.size(), sz2 = expr2.size();
1409  return sz1 > sz2 ? sz1 : sz2;
1410  }
1411 
1413  bool hasFastAccess() const {
1414  return expr1.hasFastAccess() && expr2.hasFastAccess();
1415  }
1416 
1418  bool isPassive() const {
1419  return expr1.isPassive() && expr2.isPassive();
1420  }
1421 
1423  bool updateValue() const {
1424  return expr1.updateValue() && expr2.updateValue();
1425  }
1426 
1428  void cache() const {}
1429 
1431  const value_type val() const {
1432  using std::pow;
1433  return pow(expr1.val(), expr2.val());
1434  }
1435 
1437  const value_type dx(int i) const {
1438  using std::pow; using std::log;
1439  const int sz1 = expr1.size(), sz2 = expr2.size();
1440  if (sz1 > 0 && sz2 > 0)
1441  return expr1.val() == scalar_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()));
1442  else if (sz1 > 0)
1443  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1444  // It seems less accurate and caused convergence problems in some codes
1445  return expr2.val() == scalar_type(1.0) ? expr1.dx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val()));
1446  else
1447  return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
1448  }
1449 
1451  const value_type fastAccessDx(int i) const {
1452  using std::pow; using std::log;
1453  return expr1.val() == scalar_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()));
1454  }
1455 
1456  protected:
1457 
1458  const ExprT1& expr1;
1459  const ExprT2& expr2;
1460 
1461  };
1462 
1463  template <typename ExprT1, typename T2>
1464  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Scalar >,
1465  ExprSpecDefault > {
1466 
1467  public:
1468 
1471  typedef typename ExprT1::value_type value_type_1;
1473  typedef typename Sacado::Promote<value_type_1,
1475 
1476  typedef typename ExprT1::scalar_type scalar_type_1;
1478  typedef typename Sacado::Promote<scalar_type_1,
1480 
1481  typedef typename ExprT1::base_expr_type base_expr_type_1;
1483  typedef typename Sacado::Promote<base_expr_type_1,
1485 
1487  Expr(const ExprT1& expr1_, const ConstT& c_) :
1488  expr1(expr1_), c(c_) {}
1489 
1491  int size() const {
1492  return expr1.size();
1493  }
1494 
1496  bool hasFastAccess() const {
1497  return expr1.hasFastAccess();
1498  }
1499 
1501  bool isPassive() const {
1502  return expr1.isPassive();
1503  }
1504 
1506  bool updateValue() const { return expr1.updateValue(); }
1507 
1509  void cache() const {}
1510 
1512  const value_type val() const {
1513  using std::pow;
1514  return pow(expr1.val(), c.val());
1515  }
1516 
1518  const value_type dx(int i) const {
1519  using std::pow;
1520  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1521  // It seems less accurate and caused convergence problems in some codes
1522  return c.val() == scalar_type(1.0) ? expr1.dx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val()));
1523  }
1524 
1526  const value_type fastAccessDx(int i) const {
1527  using std::pow;
1528  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1529  // It seems less accurate and caused convergence problems in some codes
1530  return c.val() == scalar_type(1.0) ? expr1.fastAccessDx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()));
1531  }
1532 
1533  protected:
1534 
1535  const ExprT1& expr1;
1537  };
1538 
1539  template <typename T1, typename ExprT2>
1540  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Scalar >,
1541  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  //
1613  // Specialization for nested derivatives. This version does not use
1614  // if_then_else/ternary-operator on the base so that nested derivatives
1615  // are correct.
1616  //
1617  template <typename ExprT1, typename ExprT2>
1618  class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Nested >,
1619  ExprSpecDefault > {
1620 
1621  public:
1622 
1623  typedef typename ExprT1::value_type value_type_1;
1624  typedef typename ExprT2::value_type value_type_2;
1625  typedef typename Sacado::Promote<value_type_1,
1627 
1628  typedef typename ExprT1::scalar_type scalar_type_1;
1629  typedef typename ExprT2::scalar_type scalar_type_2;
1630  typedef typename Sacado::Promote<scalar_type_1,
1632 
1633  typedef typename ExprT1::base_expr_type base_expr_type_1;
1634  typedef typename ExprT2::base_expr_type base_expr_type_2;
1635  typedef typename Sacado::Promote<base_expr_type_1,
1637 
1639  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1640  expr1(expr1_), expr2(expr2_) {}
1641 
1643  int size() const {
1644  int sz1 = expr1.size(), sz2 = expr2.size();
1645  return sz1 > sz2 ? sz1 : sz2;
1646  }
1647 
1649  bool hasFastAccess() const {
1650  return expr1.hasFastAccess() && expr2.hasFastAccess();
1651  }
1652 
1654  bool isPassive() const {
1655  return expr1.isPassive() && expr2.isPassive();
1656  }
1657 
1659  bool updateValue() const {
1660  return expr1.updateValue() && expr2.updateValue();
1661  }
1662 
1664  void cache() const {}
1665 
1667  const value_type val() const {
1668  using std::pow;
1669  return pow(expr1.val(), expr2.val());
1670  }
1671 
1673  const value_type dx(int i) const {
1674  using std::pow; using std::log;
1675  const int sz1 = expr1.size(), sz2 = expr2.size();
1676  if (sz1 > 0 && sz2 > 0)
1677  return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1678  else if (sz1 > 0)
1679  return expr2.val() == scalar_type(0.0) ? value_type(0.0) : value_type((expr2.val()*expr1.dx(i))*pow(expr1.val(),expr2.val()-scalar_type(1.0)));
1680  else
1681  return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1682  }
1683 
1685  const value_type fastAccessDx(int i) const {
1686  using std::pow; using std::log;
1687  return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1688  }
1689 
1690  protected:
1691 
1692  const ExprT1& expr1;
1693  const ExprT2& expr2;
1694 
1695  };
1696 
1697  template <typename ExprT1, typename T2>
1698  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Nested >,
1699  ExprSpecDefault > {
1700 
1701  public:
1702 
1705  typedef typename ExprT1::value_type value_type_1;
1707  typedef typename Sacado::Promote<value_type_1,
1709 
1710  typedef typename ExprT1::scalar_type scalar_type_1;
1712  typedef typename Sacado::Promote<scalar_type_1,
1714 
1715  typedef typename ExprT1::base_expr_type base_expr_type_1;
1717  typedef typename Sacado::Promote<base_expr_type_1,
1719 
1721  Expr(const ExprT1& expr1_, const ConstT& c_) :
1722  expr1(expr1_), c(c_) {}
1723 
1725  int size() const {
1726  return expr1.size();
1727  }
1728 
1730  bool hasFastAccess() const {
1731  return expr1.hasFastAccess();
1732  }
1733 
1735  bool isPassive() const {
1736  return expr1.isPassive();
1737  }
1738 
1740  bool updateValue() const { return expr1.updateValue(); }
1741 
1743  void cache() const {}
1744 
1746  const value_type val() const {
1747  using std::pow;
1748  return pow(expr1.val(), c.val());
1749  }
1750 
1752  const value_type dx(int i) const {
1753  using std::pow;
1754  return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
1755  }
1756 
1758  const value_type fastAccessDx(int i) const {
1759  using std::pow;
1760  return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
1761  }
1762 
1763  protected:
1764 
1765  const ExprT1& expr1;
1767  };
1768 
1769  template <typename T1, typename ExprT2>
1770  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Nested >,
1771  ExprSpecDefault > {
1772 
1773  public:
1774 
1778  typedef typename ExprT2::value_type value_type_2;
1779  typedef typename Sacado::Promote<value_type_1,
1781 
1783  typedef typename ExprT2::scalar_type scalar_type_2;
1784  typedef typename Sacado::Promote<scalar_type_1,
1786 
1788  typedef typename ExprT2::base_expr_type base_expr_type_2;
1789  typedef typename Sacado::Promote<base_expr_type_1,
1791 
1792 
1794  Expr(const ConstT& c_, const ExprT2& expr2_) :
1795  c(c_), expr2(expr2_) {}
1796 
1798  int size() const {
1799  return expr2.size();
1800  }
1801 
1803  bool hasFastAccess() const {
1804  return expr2.hasFastAccess();
1805  }
1806 
1808  bool isPassive() const {
1809  return expr2.isPassive();
1810  }
1811 
1813  bool updateValue() const { return expr2.updateValue(); }
1814 
1816  void cache() const {}
1817 
1819  const value_type val() const {
1820  using std::pow;
1821  return pow(c.val(), expr2.val());
1822  }
1823 
1825  const value_type dx(int i) const {
1826  using std::pow; using std::log;
1827  return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
1828  }
1829 
1831  const value_type fastAccessDx(int i) const {
1832  using std::pow; using std::log;
1833  return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
1834  }
1835 
1836  protected:
1837 
1839  const ExprT2& expr2;
1840  };
1841 
1842  //
1843  // Specialization for nested derivatives. This version does not use
1844  // if_then_else/ternary-operator on the base so that nested derivatives
1845  // are correct.
1846  //
1847  template <typename ExprT1, typename ExprT2>
1848  class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::NestedSimd >,
1849  ExprSpecDefault > {
1850 
1851  public:
1852 
1853  typedef typename ExprT1::value_type value_type_1;
1854  typedef typename ExprT2::value_type value_type_2;
1855  typedef typename Sacado::Promote<value_type_1,
1857 
1858  typedef typename ExprT1::scalar_type scalar_type_1;
1859  typedef typename ExprT2::scalar_type scalar_type_2;
1860  typedef typename Sacado::Promote<scalar_type_1,
1862 
1863  typedef typename ExprT1::base_expr_type base_expr_type_1;
1864  typedef typename ExprT2::base_expr_type base_expr_type_2;
1865  typedef typename Sacado::Promote<base_expr_type_1,
1867 
1869  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1870  expr1(expr1_), expr2(expr2_) {}
1871 
1873  int size() const {
1874  int sz1 = expr1.size(), sz2 = expr2.size();
1875  return sz1 > sz2 ? sz1 : sz2;
1876  }
1877 
1879  bool hasFastAccess() const {
1880  return expr1.hasFastAccess() && expr2.hasFastAccess();
1881  }
1882 
1884  bool isPassive() const {
1885  return expr1.isPassive() && expr2.isPassive();
1886  }
1887 
1889  bool updateValue() const {
1890  return expr1.updateValue() && expr2.updateValue();
1891  }
1892 
1894  void cache() const {}
1895 
1897  const value_type val() const {
1898  using std::pow;
1899  return pow(expr1.val(), expr2.val());
1900  }
1901 
1903  const value_type dx(int i) const {
1904  using std::pow; using std::log; using Sacado::if_then_else;
1905  const int sz1 = expr1.size(), sz2 = expr2.size();
1906  if (sz1 > 0 && sz2 > 0)
1907  return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1908  else if (sz1 > 0)
1909  return if_then_else( expr2.val() == scalar_type(0.0), value_type(0.0), value_type((expr2.val()*expr1.dx(i))*pow(expr1.val(),expr2.val()-scalar_type(1.0))));
1910  else
1911  return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1912  }
1913 
1915  const value_type fastAccessDx(int i) const {
1916  using std::pow; using std::log;
1917  return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1918  }
1919 
1920  protected:
1921 
1922  const ExprT1& expr1;
1923  const ExprT2& expr2;
1924 
1925  };
1926 
1927  template <typename ExprT1, typename T2>
1928  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::NestedSimd >,
1929  ExprSpecDefault > {
1930 
1931  public:
1932 
1935  typedef typename ExprT1::value_type value_type_1;
1937  typedef typename Sacado::Promote<value_type_1,
1939 
1940  typedef typename ExprT1::scalar_type scalar_type_1;
1942  typedef typename Sacado::Promote<scalar_type_1,
1944 
1945  typedef typename ExprT1::base_expr_type base_expr_type_1;
1947  typedef typename Sacado::Promote<base_expr_type_1,
1949 
1951  Expr(const ExprT1& expr1_, const ConstT& c_) :
1952  expr1(expr1_), c(c_) {}
1953 
1955  int size() const {
1956  return expr1.size();
1957  }
1958 
1960  bool hasFastAccess() const {
1961  return expr1.hasFastAccess();
1962  }
1963 
1965  bool isPassive() const {
1966  return expr1.isPassive();
1967  }
1968 
1970  bool updateValue() const { return expr1.updateValue(); }
1971 
1973  void cache() const {}
1974 
1976  const value_type val() const {
1977  using std::pow;
1978  return pow(expr1.val(), c.val());
1979  }
1980 
1982  const value_type dx(int i) const {
1983  using std::pow; using Sacado::if_then_else;
1984  return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))));
1985  }
1986 
1988  const value_type fastAccessDx(int i) const {
1989  using std::pow; using Sacado::if_then_else;
1990  return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))));
1991  }
1992 
1993  protected:
1994 
1995  const ExprT1& expr1;
1997  };
1998 
1999  template <typename T1, typename ExprT2>
2000  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::NestedSimd >,
2001  ExprSpecDefault > {
2002 
2003  public:
2004 
2008  typedef typename ExprT2::value_type value_type_2;
2009  typedef typename Sacado::Promote<value_type_1,
2011 
2013  typedef typename ExprT2::scalar_type scalar_type_2;
2014  typedef typename Sacado::Promote<scalar_type_1,
2016 
2018  typedef typename ExprT2::base_expr_type base_expr_type_2;
2019  typedef typename Sacado::Promote<base_expr_type_1,
2021 
2022 
2024  Expr(const ConstT& c_, const ExprT2& expr2_) :
2025  c(c_), expr2(expr2_) {}
2026 
2028  int size() const {
2029  return expr2.size();
2030  }
2031 
2033  bool hasFastAccess() const {
2034  return expr2.hasFastAccess();
2035  }
2036 
2038  bool isPassive() const {
2039  return expr2.isPassive();
2040  }
2041 
2043  bool updateValue() const { return expr2.updateValue(); }
2044 
2046  void cache() const {}
2047 
2049  const value_type val() const {
2050  using std::pow;
2051  return pow(c.val(), expr2.val());
2052  }
2053 
2055  const value_type dx(int i) const {
2056  using std::pow; using std::log;
2057  return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
2058  }
2059 
2061  const value_type fastAccessDx(int i) const {
2062  using std::pow; using std::log;
2063  return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
2064  }
2065 
2066  protected:
2067 
2069  const ExprT2& expr2;
2070  };
2071 
2072  template <typename T1, typename T2>
2074  typename mpl::enable_if_c<
2077  >::type
2078  /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(PowerOp)*/
2079  pow (const Expr<T1>& expr1, const Expr<T2>& expr2)
2080  {
2081  typedef PowerOp< Expr<T1>, Expr<T2> > expr_t;
2082 
2083  return Expr<expr_t>(expr1, expr2);
2084  }
2085 
2086  template <typename T>
2088  Expr< PowerOp< Expr<T>, Expr<T> > >
2089  pow (const Expr<T>& expr1, const Expr<T>& expr2)
2090  {
2091  typedef PowerOp< Expr<T>, Expr<T> > expr_t;
2092 
2093  return Expr<expr_t>(expr1, expr2);
2094  }
2095 
2096  template <typename T>
2098  Expr< PowerOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
2099  pow (const typename Expr<T>::value_type& c,
2100  const Expr<T>& expr)
2101  {
2103  typedef PowerOp< ConstT, Expr<T> > expr_t;
2104 
2105  return Expr<expr_t>(ConstT(c), expr);
2106  }
2107 
2108  template <typename T>
2110  Expr< PowerOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
2111  pow (const Expr<T>& expr,
2112  const typename Expr<T>::value_type& c)
2113  {
2115  typedef PowerOp< Expr<T>, ConstT > expr_t;
2116 
2117  return Expr<expr_t>(expr, ConstT(c));
2118  }
2119 
2120  template <typename T>
2123  pow (const typename Expr<T>::scalar_type& c,
2124  const Expr<T>& expr)
2125  {
2127  typedef PowerOp< ConstT, Expr<T> > expr_t;
2128 
2129  return Expr<expr_t>(ConstT(c), expr);
2130  }
2131 
2132  template <typename T>
2135  pow (const Expr<T>& expr,
2136  const typename Expr<T>::scalar_type& c)
2137  {
2139  typedef PowerOp< Expr<T>, ConstT > expr_t;
2140 
2141  return Expr<expr_t>(expr, ConstT(c));
2142  }
2143  }
2144 }
2145 
2146 //--------------------------if_then_else operator -----------------------
2147 // Can't use the above macros because it is a ternary operator (sort of).
2148 // Also, relies on C++11
2149 
2150 
2151 namespace Sacado {
2152  namespace Fad {
2153 
2154  template <typename CondT, typename ExprT1, typename ExprT2>
2155  class IfThenElseOp {};
2156 
2157  template <typename CondT, typename ExprT1, typename ExprT2>
2158  struct ExprSpec< IfThenElseOp< CondT, ExprT1, ExprT2 > > {
2159  typedef typename ExprSpec<ExprT1>::type type;
2160  };
2161 
2162  template <typename CondT, typename ExprT1, typename ExprT2>
2163  class Expr< IfThenElseOp< CondT, ExprT1, ExprT2 >,ExprSpecDefault > {
2164 
2165  public:
2166 
2167  typedef typename ExprT1::value_type value_type_1;
2168  typedef typename ExprT2::value_type value_type_2;
2169  typedef typename Sacado::Promote<value_type_1,
2171 
2172  typedef typename ExprT1::scalar_type scalar_type_1;
2173  typedef typename ExprT2::scalar_type scalar_type_2;
2174  typedef typename Sacado::Promote<scalar_type_1,
2176 
2177  typedef typename ExprT1::base_expr_type base_expr_type_1;
2178  typedef typename ExprT2::base_expr_type base_expr_type_2;
2179  typedef typename Sacado::Promote<base_expr_type_1,
2181 
2183  Expr(const CondT& cond_, const ExprT1& expr1_, const ExprT2& expr2_) :
2184  cond(cond_), expr1(expr1_), expr2(expr2_) {}
2185 
2187  int size() const {
2188  int sz1 = expr1.size(), sz2 = expr2.size();
2189  return sz1 > sz2 ? sz1 : sz2;
2190  }
2191 
2193  bool hasFastAccess() const {
2194  return expr1.hasFastAccess() && expr2.hasFastAccess();
2195  }
2196 
2198  bool isPassive() const {
2199  return expr1.isPassive() && expr2.isPassive();
2200  }
2201 
2203  bool updateValue() const {
2204  return expr1.updateValue() && expr2.updateValue();
2205  }
2206 
2208  void cache() const {}
2209 
2211  const value_type val() const {
2212  using Sacado::if_then_else;
2213  return if_then_else( cond, expr1.val(), expr2.val() );
2214  }
2215 
2217  const value_type dx(int i) const {
2218  using Sacado::if_then_else;
2219  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
2220  }
2221 
2223  const value_type fastAccessDx(int i) const {
2224  using Sacado::if_then_else;
2225  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
2226  }
2227 
2228  protected:
2229 
2230  const CondT& cond;
2231  const ExprT1& expr1;
2232  const ExprT2& expr2;
2233 
2234  };
2235 
2236  template <typename CondT, typename ExprT1, typename T2>
2237  struct ExprSpec< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> > > {
2238  typedef typename ExprSpec<ExprT1>::type type;
2239  };
2240 
2241  template <typename CondT, typename ExprT1, typename T2>
2242  class Expr< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
2243 
2244  public:
2245 
2248  typedef typename ExprT1::value_type value_type_1;
2250  typedef typename Sacado::Promote<value_type_1,
2252 
2253  typedef typename ExprT1::scalar_type scalar_type_1;
2255  typedef typename Sacado::Promote<scalar_type_1,
2257 
2258  typedef typename ExprT1::base_expr_type base_expr_type_1;
2260  typedef typename Sacado::Promote<base_expr_type_1,
2262 
2264  Expr(const CondT& cond_, const ExprT1& expr1_, const ConstT& c_) :
2265  cond(cond_), expr1(expr1_), c(c_) {}
2266 
2268  int size() const {
2269  return expr1.size();
2270  }
2271 
2273  bool hasFastAccess() const {
2274  return expr1.hasFastAccess();
2275  }
2276 
2278  bool isPassive() const {
2279  return expr1.isPassive();
2280  }
2281 
2283  bool updateValue() const { return expr1.updateValue(); }
2284 
2286  void cache() const {}
2287 
2289  const value_type val() const {
2290  using Sacado::if_then_else;
2291  return if_then_else( cond, expr1.val(), c.val() );
2292  }
2293 
2295  const value_type dx(int i) const {
2296  using Sacado::if_then_else;
2297  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
2298  }
2299 
2301  const value_type fastAccessDx(int i) const {
2302  using Sacado::if_then_else;
2303  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
2304  }
2305 
2306  protected:
2307 
2308  const CondT& cond;
2309  const ExprT1& expr1;
2311  };
2312 
2313  template <typename CondT, typename T1, typename ExprT2>
2314  struct ExprSpec< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 > > {
2315  typedef typename ExprSpec<ExprT2>::type type;
2316  };
2317 
2318  template <typename CondT, typename T1, typename ExprT2>
2319  class Expr< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
2320 
2321  public:
2322 
2326  typedef typename ExprT2::value_type value_type_2;
2327  typedef typename Sacado::Promote<value_type_1,
2329 
2331  typedef typename ExprT2::scalar_type scalar_type_2;
2332  typedef typename Sacado::Promote<scalar_type_1,
2334 
2336  typedef typename ExprT2::base_expr_type base_expr_type_2;
2337  typedef typename Sacado::Promote<base_expr_type_1,
2339 
2341  Expr(const CondT& cond_, const ConstT& c_, const ExprT2& expr2_) :
2342  cond(cond_), c(c_), expr2(expr2_) {}
2343 
2345  int size() const {
2346  return expr2.size();
2347  }
2348 
2350  bool hasFastAccess() const {
2351  return expr2.hasFastAccess();
2352  }
2353 
2355  bool isPassive() const {
2356  return expr2.isPassive();
2357  }
2358 
2360  bool updateValue() const { return expr2.updateValue(); }
2361 
2363  void cache() const {}
2364 
2366  const value_type val() const {
2367  using Sacado::if_then_else;
2368  return if_then_else( cond, c.val(), expr2.val() );
2369  }
2370 
2372  const value_type dx(int i) const {
2373  using Sacado::if_then_else;
2374  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
2375  }
2376 
2378  const value_type fastAccessDx(int i) const {
2379  using Sacado::if_then_else;
2380  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
2381  }
2382 
2383  protected:
2384 
2385  const CondT& cond;
2387  const ExprT2& expr2;
2388  };
2389 
2390  template <typename CondT, typename T1, typename T2>
2395  >::type
2396  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
2397  {
2398  typedef IfThenElseOp< CondT, T1, T2 > expr_t;
2399 
2400  return Expr<expr_t>(cond, expr1, expr2);
2401  }
2402 
2403  template <typename CondT, typename T>
2405  Expr< IfThenElseOp< CondT, Expr<T>, Expr<T> > >
2406  if_then_else (const CondT& cond, const Expr<T>& expr1, const Expr<T>& expr2)
2407  {
2408  typedef IfThenElseOp< CondT, Expr<T>, Expr<T> > expr_t;
2409 
2410  return Expr<expr_t>(cond, expr1, expr2);
2411  }
2412 
2413  template <typename CondT, typename T>
2415  Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::value_type>,
2416  Expr<T> > >
2417  if_then_else (const CondT& cond, const typename Expr<T>::value_type& c,
2418  const Expr<T>& expr)
2419  {
2421  typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2422 
2423  return Expr<expr_t>(cond, ConstT(c), expr);
2424  }
2425 
2426  template <typename CondT, typename T>
2428  Expr< IfThenElseOp< CondT, Expr<T>,
2429  ConstExpr<typename Expr<T>::value_type> > >
2430  if_then_else (const CondT& cond, const Expr<T>& expr,
2431  const typename Expr<T>::value_type& c)
2432  {
2434  typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2435 
2436  return Expr<expr_t>(cond, expr, ConstT(c));
2437  }
2438 
2439  template <typename CondT, typename T>
2441  typename mpl::disable_if<
2442  std::is_same< typename Expr<T>::value_type,
2443  typename Expr<T>::scalar_type>,
2444  Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::scalar_type>,
2445  Expr<T> > >
2446  >::type
2447  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
2448  const Expr<T>& expr)
2449  {
2451  typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2452 
2453  return Expr<expr_t>(cond, ConstT(c), expr);
2454  }
2455 
2456  template <typename CondT, typename T>
2458  typename mpl::disable_if<
2459  std::is_same< typename Expr<T>::value_type,
2460  typename Expr<T>::scalar_type>,
2461  Expr< IfThenElseOp< CondT, Expr<T>,
2462  ConstExpr<typename Expr<T>::scalar_type> > >
2463  >::type
2464  if_then_else (const CondT& cond, const Expr<T>& expr,
2465  const typename Expr<T>::scalar_type& c)
2466  {
2468  typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2469 
2470  return Expr<expr_t>(cond, expr, ConstT(c));
2471  }
2472  }
2473 }
2474 
2475 //-------------------------- Relational Operators -----------------------
2476 
2477 namespace Sacado {
2478  namespace Fad {
2479  template <typename T1, typename T2 = T1>
2481  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
2482  };
2483  }
2484 }
2485 
2486 #define FAD_RELOP_MACRO(OP) \
2487 namespace Sacado { \
2488  namespace Fad { \
2489  template <typename ExprT1, typename ExprT2> \
2490  SACADO_INLINE_FUNCTION \
2491  typename ConditionalReturnType<typename Expr<ExprT1>::value_type, \
2492  typename Expr<ExprT2>::value_type>::type \
2493  operator OP (const Expr<ExprT1>& expr1, \
2494  const Expr<ExprT2>& expr2) \
2495  { \
2496  return expr1.val() OP expr2.val(); \
2497  } \
2498  \
2499  template <typename ExprT2> \
2500  SACADO_INLINE_FUNCTION \
2501  typename ConditionalReturnType<typename Expr<ExprT2>::value_type>::type \
2502  operator OP (const typename Expr<ExprT2>::value_type& a, \
2503  const Expr<ExprT2>& expr2) \
2504  { \
2505  return a OP expr2.val(); \
2506  } \
2507  \
2508  template <typename ExprT1> \
2509  SACADO_INLINE_FUNCTION \
2510  typename ConditionalReturnType<typename Expr<ExprT1>::value_type>::type \
2511  operator OP (const Expr<ExprT1>& expr1, \
2512  const typename Expr<ExprT1>::value_type& b) \
2513  { \
2514  return expr1.val() OP b; \
2515  } \
2516  } \
2517 }
2518 
2519 FAD_RELOP_MACRO(==)
2520 FAD_RELOP_MACRO(!=)
2521 FAD_RELOP_MACRO(<)
2522 FAD_RELOP_MACRO(>)
2523 FAD_RELOP_MACRO(<=)
2524 FAD_RELOP_MACRO(>=)
2525 FAD_RELOP_MACRO(<<=)
2526 FAD_RELOP_MACRO(>>=)
2527 FAD_RELOP_MACRO(&)
2528 FAD_RELOP_MACRO(|)
2529 
2530 #undef FAD_RELOP_MACRO
2531 
2532 namespace Sacado {
2533 
2534  namespace Fad {
2535 
2536  template <typename ExprT>
2538  bool operator ! (const Expr<ExprT>& expr)
2539  {
2540  return ! expr.val();
2541  }
2542 
2543  } // namespace Fad
2544 
2545 } // namespace Sacado
2546 
2547 //-------------------------- Boolean Operators -----------------------
2548 namespace Sacado {
2549 
2550  namespace Fad {
2551 
2552  template <typename ExprT>
2554  bool toBool(const Expr<ExprT>& x) {
2555  bool is_zero = (x.val() == 0.0);
2556  for (int i=0; i<x.size(); i++)
2557  is_zero = is_zero && (x.dx(i) == 0.0);
2558  return !is_zero;
2559  }
2560 
2561  } // namespace Fad
2562 
2563 } // namespace Sacado
2564 
2565 #define FAD_BOOL_MACRO(OP) \
2566 namespace Sacado { \
2567  namespace Fad { \
2568  template <typename ExprT1, typename ExprT2> \
2569  SACADO_INLINE_FUNCTION \
2570  bool \
2571  operator OP (const Expr<ExprT1>& expr1, \
2572  const Expr<ExprT2>& expr2) \
2573  { \
2574  return toBool(expr1) OP toBool(expr2); \
2575  } \
2576  \
2577  template <typename ExprT2> \
2578  SACADO_INLINE_FUNCTION \
2579  bool \
2580  operator OP (const typename Expr<ExprT2>::value_type& a, \
2581  const Expr<ExprT2>& expr2) \
2582  { \
2583  return a OP toBool(expr2); \
2584  } \
2585  \
2586  template <typename ExprT1> \
2587  SACADO_INLINE_FUNCTION \
2588  bool \
2589  operator OP (const Expr<ExprT1>& expr1, \
2590  const typename Expr<ExprT1>::value_type& b) \
2591  { \
2592  return toBool(expr1) OP b; \
2593  } \
2594  } \
2595 }
2596 
2597 FAD_BOOL_MACRO(&&)
2598 FAD_BOOL_MACRO(||)
2599 
2600 #undef FAD_BOOL_MACRO
2601 
2602 //-------------------------- I/O Operators -----------------------
2603 
2604 namespace Sacado {
2605 
2606  namespace Fad {
2607 
2608  template <typename ExprT>
2609  std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
2610  os << x.val() << " [";
2611 
2612  for (int i=0; i< x.size(); i++) {
2613  os << " " << x.dx(i);
2614  }
2615 
2616  os << " ]";
2617  return os;
2618  }
2619 
2620  } // namespace Fad
2621 
2622 } // namespace Sacado
2623 
2624 #endif // SACADO_FAD_OPS_HPP
Wrapper for a generic expression template.
cbrt(expr.val())
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
SACADO_INLINE_FUNCTION Expr(const CondT &cond_, const ConstT &c_, const ExprT2 &expr2_)
asinh(expr.val())
char c_
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
asin(expr.val())
Sacado::Promote< scalar_type_1, scalar_type_2 >::type scalar_type
cosh(expr.val())
expr expr dx(i)
abs(expr.val())
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
SACADO_INLINE_FUNCTION const value_type fastAccessDx(int i) const
SACADO_INLINE_FUNCTION Expr< SafeSqrtOp< Expr< T > > > safe_sqrt(const Expr< T > &)
atanh(expr.val())
expr expr CoshOp
expr expr ATanhOp
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ConstT &c_)
#define SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP)
expr expr TanhOp
SimpleFad< ValueT > sqrt(const SimpleFad< ValueT > &a)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const CondT &cond_, const ExprT1 &expr1_, const ConstT &c_)
expr expr SqrtOp
expr expr ASinhOp
decltype(std::declval< T1 >()==std::declval< T2 >()) typedef type
SACADO_INLINE_FUNCTION mpl::enable_if_c< IsFadExpr< T1 >::value &&IsFadExpr< T2 >::value &&ExprLevel< T1 >::value==ExprLevel< T2 >::value, Expr< IfThenElseOp< CondT, T1, T2 > > >::type if_then_else(const CondT &cond, const T1 &expr1, const T2 &expr2)
atan(expr.val())
SACADO_INLINE_FUNCTION const value_type fastAccessDx(int i) const
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
Sacado::Promote< value_type_1, value_type_2 >::type value_type
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
SimpleFad< ValueT > log(const SimpleFad< ValueT > &a)
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 T
Definition: Sacado_rad.hpp:553
SACADO_INLINE_FUNCTION bool toBool(const Expr< ExprT > &x)
tanh(expr.val())
expr expr CosOp
Determine whether a given type is an expression.
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
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
#define FAD_BINARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, CDX1, CDX2, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
#define FAD_BOOL_MACRO(OP)
SACADO_INLINE_FUNCTION Expr(const CondT &cond_, const ExprT1 &expr1_, const ExprT2 &expr2_)
expr expr ATanOp
#define T2(r, f)
Definition: Sacado_rad.hpp:558
SACADO_INLINE_FUNCTION T if_then_else(const Cond cond, const T &a, const T &b)
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
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
T2 base_expr_type
Typename of base-expressions.
#define T1(r, f)
Definition: Sacado_rad.hpp:583
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
atan2(expr1.val(), expr2.val())
Meta-function for determining nesting with an expression.
sin(expr.val())
int value
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
#define SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP)
SACADO_INLINE_FUNCTION const value_type dx(int i) const
log(expr.val())
SACADO_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
expr expr ACoshOp
expr expr Log10Op
expr expr SinhOp
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
acosh(expr.val())
SimpleFad< ValueT > operator*(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
SACADO_INLINE_FUNCTION Expr(const ConstT &c_, const ExprT2 &expr2_)
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
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
T2 value_type
Typename of argument values.
SACADO_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)
#define FAD_RELOP_MACRO(OP)
#define SACADO_INLINE_FUNCTION
exp(expr.val())
expr expr expr ExpOp
fabs(expr.val())
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
expr expr AbsOp
expr expr TanOp
log10(expr.val())
Base template specification for Promote.
cos(expr.val())
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 PowerOp