Sacado Package Browser (Single Doxygen Collection)  Version of the Day
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  SACADO_INLINE_FUNCTION \
81  explicit Expr(const ExprT& expr_) : expr(expr_) {} \
82  \
83  SACADO_INLINE_FUNCTION \
84  int size() const { return expr.size(); } \
85  \
86  SACADO_INLINE_FUNCTION \
87  bool hasFastAccess() const { return expr.hasFastAccess(); } \
88  \
89  SACADO_INLINE_FUNCTION \
90  bool isPassive() const { return expr.isPassive();} \
91  \
92  SACADO_INLINE_FUNCTION \
93  bool updateValue() const { return expr.updateValue(); } \
94  \
95  SACADO_INLINE_FUNCTION \
96  void cache() const {} \
97  \
98  SACADO_INLINE_FUNCTION \
99  value_type val() const { \
100  USING \
101  return VALUE; \
102  } \
103  \
104  SACADO_INLINE_FUNCTION \
105  value_type dx(int i) const { \
106  USING \
107  return DX; \
108  } \
109  \
110  SACADO_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  SACADO_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()))
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; using Sacado::if_then_else;
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; using Sacado::if_then_else;
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  SACADO_INLINE_FUNCTION \
444  Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
445  expr1(expr1_), expr2(expr2_) {} \
446  \
447  SACADO_INLINE_FUNCTION \
448  int size() const { \
449  int sz1 = expr1.size(), sz2 = expr2.size(); \
450  return sz1 > sz2 ? sz1 : sz2; \
451  } \
452  \
453  SACADO_INLINE_FUNCTION \
454  bool hasFastAccess() const { \
455  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
456  } \
457  \
458  SACADO_INLINE_FUNCTION \
459  bool isPassive() const { \
460  return expr1.isPassive() && expr2.isPassive(); \
461  } \
462  \
463  SACADO_INLINE_FUNCTION \
464  bool updateValue() const { \
465  return expr1.updateValue() && expr2.updateValue(); \
466  } \
467  \
468  SACADO_INLINE_FUNCTION \
469  void cache() const {} \
470  \
471  SACADO_INLINE_FUNCTION \
472  const value_type val() const { \
473  USING \
474  return VALUE; \
475  } \
476  \
477  SACADO_INLINE_FUNCTION \
478  const value_type dx(int i) const { \
479  USING \
480  return DX; \
481  } \
482  \
483  SACADO_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  SACADO_INLINE_FUNCTION \
524  Expr(const ExprT1& expr1_, const ConstT& c_) : \
525  expr1(expr1_), c(c_) {} \
526  \
527  SACADO_INLINE_FUNCTION \
528  int size() const { \
529  return expr1.size(); \
530  } \
531  \
532  SACADO_INLINE_FUNCTION \
533  bool hasFastAccess() const { \
534  return expr1.hasFastAccess(); \
535  } \
536  \
537  SACADO_INLINE_FUNCTION \
538  bool isPassive() const { \
539  return expr1.isPassive(); \
540  } \
541  \
542  SACADO_INLINE_FUNCTION \
543  bool updateValue() const { return expr1.updateValue(); } \
544  \
545  SACADO_INLINE_FUNCTION \
546  void cache() const {} \
547  \
548  SACADO_INLINE_FUNCTION \
549  const value_type val() const { \
550  USING \
551  return VAL_CONST_DX_2; \
552  } \
553  \
554  SACADO_INLINE_FUNCTION \
555  const value_type dx(int i) const { \
556  USING \
557  return CONST_DX_2; \
558  } \
559  \
560  SACADO_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  SACADO_INLINE_FUNCTION \
601  Expr(const ConstT& c_, const ExprT2& expr2_) : \
602  c(c_), expr2(expr2_) {} \
603  \
604  SACADO_INLINE_FUNCTION \
605  int size() const { \
606  return expr2.size(); \
607  } \
608  \
609  SACADO_INLINE_FUNCTION \
610  bool hasFastAccess() const { \
611  return expr2.hasFastAccess(); \
612  } \
613  \
614  SACADO_INLINE_FUNCTION \
615  bool isPassive() const { \
616  return expr2.isPassive(); \
617  } \
618  \
619  SACADO_INLINE_FUNCTION \
620  bool updateValue() const { return expr2.updateValue(); } \
621  \
622  SACADO_INLINE_FUNCTION \
623  void cache() const {} \
624  \
625  SACADO_INLINE_FUNCTION \
626  const value_type val() const { \
627  USING \
628  return VAL_CONST_DX_1; \
629  } \
630  \
631  SACADO_INLINE_FUNCTION \
632  const value_type dx(int i) const { \
633  USING \
634  return CONST_DX_1; \
635  } \
636  \
637  SACADO_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  SACADO_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  SACADO_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  SACADO_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  SACADO_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  SACADO_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  SACADO_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>
1093  operator* (const typename Expr<T>::value_type& c,
1094  const Expr<T>& expr)
1095  {
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  {
1109  typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1110 
1111  return Expr<expr_t>(expr, ConstT(c));
1112  }
1113 
1114  template <typename T>
1117  operator* (const typename Expr<T>::scalar_type& c,
1118  const Expr<T>& expr)
1119  {
1121  typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1122 
1123  return Expr<expr_t>(ConstT(c), expr);
1124  }
1125 
1126  template <typename T>
1129  operator* (const Expr<T>& expr,
1130  const typename Expr<T>::scalar_type& c)
1131  {
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, typename Impl>
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, PowerImpl::Simd >, 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; using Sacado::if_then_else;
1225  const int sz1 = expr1.size(), sz2 = expr2.size();
1226  if (sz1 > 0 && sz2 > 0)
1227  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())) );
1228  else if (sz1 > 0)
1229  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1230  // It seems less accurate and caused convergence problems in some codes
1231  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())) ));
1232  else
1233  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())) );
1234  }
1235 
1237  const value_type fastAccessDx(int i) const {
1238  using std::pow; using std::log; using Sacado::if_then_else;
1239  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())) );
1240  }
1241 
1242  protected:
1243 
1244  const ExprT1& expr1;
1245  const ExprT2& expr2;
1246 
1247  };
1248 
1249  template <typename ExprT1, typename T2>
1250  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Simd >,
1251  ExprSpecDefault > {
1252 
1253  public:
1254 
1257  typedef typename ExprT1::value_type value_type_1;
1259  typedef typename Sacado::Promote<value_type_1,
1261 
1262  typedef typename ExprT1::scalar_type scalar_type_1;
1264  typedef typename Sacado::Promote<scalar_type_1,
1266 
1267  typedef typename ExprT1::base_expr_type base_expr_type_1;
1269  typedef typename Sacado::Promote<base_expr_type_1,
1271 
1273  Expr(const ExprT1& expr1_, const ConstT& c_) :
1274  expr1(expr1_), c(c_) {}
1275 
1277  int size() const {
1278  return expr1.size();
1279  }
1280 
1282  bool hasFastAccess() const {
1283  return expr1.hasFastAccess();
1284  }
1285 
1287  bool isPassive() const {
1288  return expr1.isPassive();
1289  }
1290 
1292  bool updateValue() const { return expr1.updateValue(); }
1293 
1295  void cache() const {}
1296 
1298  const value_type val() const {
1299  using std::pow;
1300  return pow(expr1.val(), c.val());
1301  }
1302 
1304  const value_type dx(int i) const {
1305  using std::pow; using Sacado::if_then_else;
1306  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1307  // It seems less accurate and caused convergence problems in some codes
1308  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())) ));
1309  }
1310 
1312  const value_type fastAccessDx(int i) const {
1313  using std::pow; using Sacado::if_then_else;
1314  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1315  // It seems less accurate and caused convergence problems in some codes
1316  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()))));
1317  }
1318 
1319  protected:
1320 
1321  const ExprT1& expr1;
1323  };
1324 
1325  template <typename T1, typename ExprT2>
1326  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Simd >,
1327  ExprSpecDefault > {
1328 
1329  public:
1330 
1334  typedef typename ExprT2::value_type value_type_2;
1335  typedef typename Sacado::Promote<value_type_1,
1337 
1339  typedef typename ExprT2::scalar_type scalar_type_2;
1340  typedef typename Sacado::Promote<scalar_type_1,
1342 
1344  typedef typename ExprT2::base_expr_type base_expr_type_2;
1345  typedef typename Sacado::Promote<base_expr_type_1,
1347 
1348 
1350  Expr(const ConstT& c_, const ExprT2& expr2_) :
1351  c(c_), expr2(expr2_) {}
1352 
1354  int size() const {
1355  return expr2.size();
1356  }
1357 
1359  bool hasFastAccess() const {
1360  return expr2.hasFastAccess();
1361  }
1362 
1364  bool isPassive() const {
1365  return expr2.isPassive();
1366  }
1367 
1369  bool updateValue() const { return expr2.updateValue(); }
1370 
1372  void cache() const {}
1373 
1375  const value_type val() const {
1376  using std::pow;
1377  return pow(c.val(), expr2.val());
1378  }
1379 
1381  const value_type dx(int i) const {
1382  using std::pow; using std::log; using Sacado::if_then_else;
1383  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())) );
1384  }
1385 
1387  const value_type fastAccessDx(int i) const {
1388  using std::pow; using std::log; using Sacado::if_then_else;
1389  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())) );
1390  }
1391 
1392  protected:
1393 
1395  const ExprT2& expr2;
1396  };
1397 
1398  //
1399  // Specialization for scalar types using ternary operator
1400  //
1401 
1402  template <typename ExprT1, typename ExprT2>
1403  class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Scalar >,
1404  ExprSpecDefault > {
1405 
1406  public:
1407 
1408  typedef typename ExprT1::value_type value_type_1;
1409  typedef typename ExprT2::value_type value_type_2;
1410  typedef typename Sacado::Promote<value_type_1,
1412 
1413  typedef typename ExprT1::scalar_type scalar_type_1;
1414  typedef typename ExprT2::scalar_type scalar_type_2;
1415  typedef typename Sacado::Promote<scalar_type_1,
1417 
1418  typedef typename ExprT1::base_expr_type base_expr_type_1;
1419  typedef typename ExprT2::base_expr_type base_expr_type_2;
1420  typedef typename Sacado::Promote<base_expr_type_1,
1422 
1424  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1425  expr1(expr1_), expr2(expr2_) {}
1426 
1428  int size() const {
1429  int sz1 = expr1.size(), sz2 = expr2.size();
1430  return sz1 > sz2 ? sz1 : sz2;
1431  }
1432 
1434  bool hasFastAccess() const {
1435  return expr1.hasFastAccess() && expr2.hasFastAccess();
1436  }
1437 
1439  bool isPassive() const {
1440  return expr1.isPassive() && expr2.isPassive();
1441  }
1442 
1444  bool updateValue() const {
1445  return expr1.updateValue() && expr2.updateValue();
1446  }
1447 
1449  void cache() const {}
1450 
1452  const value_type val() const {
1453  using std::pow;
1454  return pow(expr1.val(), expr2.val());
1455  }
1456 
1458  const value_type dx(int i) const {
1459  using std::pow; using std::log;
1460  const int sz1 = expr1.size(), sz2 = expr2.size();
1461  if (sz1 > 0 && sz2 > 0)
1462  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()));
1463  else if (sz1 > 0)
1464  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1465  // It seems less accurate and caused convergence problems in some codes
1466  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()));
1467  else
1468  return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
1469  }
1470 
1472  const value_type fastAccessDx(int i) const {
1473  using std::pow; using std::log;
1474  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()));
1475  }
1476 
1477  protected:
1478 
1479  const ExprT1& expr1;
1480  const ExprT2& expr2;
1481 
1482  };
1483 
1484  template <typename ExprT1, typename T2>
1485  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Scalar >,
1486  ExprSpecDefault > {
1487 
1488  public:
1489 
1492  typedef typename ExprT1::value_type value_type_1;
1494  typedef typename Sacado::Promote<value_type_1,
1496 
1497  typedef typename ExprT1::scalar_type scalar_type_1;
1499  typedef typename Sacado::Promote<scalar_type_1,
1501 
1502  typedef typename ExprT1::base_expr_type base_expr_type_1;
1504  typedef typename Sacado::Promote<base_expr_type_1,
1506 
1508  Expr(const ExprT1& expr1_, const ConstT& c_) :
1509  expr1(expr1_), c(c_) {}
1510 
1512  int size() const {
1513  return expr1.size();
1514  }
1515 
1517  bool hasFastAccess() const {
1518  return expr1.hasFastAccess();
1519  }
1520 
1522  bool isPassive() const {
1523  return expr1.isPassive();
1524  }
1525 
1527  bool updateValue() const { return expr1.updateValue(); }
1528 
1530  void cache() const {}
1531 
1533  const value_type val() const {
1534  using std::pow;
1535  return pow(expr1.val(), c.val());
1536  }
1537 
1539  const value_type dx(int i) const {
1540  using std::pow;
1541  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1542  // It seems less accurate and caused convergence problems in some codes
1543  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()));
1544  }
1545 
1547  const value_type fastAccessDx(int i) const {
1548  using std::pow;
1549  // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1550  // It seems less accurate and caused convergence problems in some codes
1551  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()));
1552  }
1553 
1554  protected:
1555 
1556  const ExprT1& expr1;
1558  };
1559 
1560  template <typename T1, typename ExprT2>
1561  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Scalar >,
1562  ExprSpecDefault > {
1563 
1564  public:
1565 
1569  typedef typename ExprT2::value_type value_type_2;
1570  typedef typename Sacado::Promote<value_type_1,
1572 
1574  typedef typename ExprT2::scalar_type scalar_type_2;
1575  typedef typename Sacado::Promote<scalar_type_1,
1577 
1579  typedef typename ExprT2::base_expr_type base_expr_type_2;
1580  typedef typename Sacado::Promote<base_expr_type_1,
1582 
1583 
1585  Expr(const ConstT& c_, const ExprT2& expr2_) :
1586  c(c_), expr2(expr2_) {}
1587 
1589  int size() const {
1590  return expr2.size();
1591  }
1592 
1594  bool hasFastAccess() const {
1595  return expr2.hasFastAccess();
1596  }
1597 
1599  bool isPassive() const {
1600  return expr2.isPassive();
1601  }
1602 
1604  bool updateValue() const { return expr2.updateValue(); }
1605 
1607  void cache() const {}
1608 
1610  const value_type val() const {
1611  using std::pow;
1612  return pow(c.val(), expr2.val());
1613  }
1614 
1616  const value_type dx(int i) const {
1617  using std::pow; using std::log;
1618  return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val()));
1619  }
1620 
1622  const value_type fastAccessDx(int i) const {
1623  using std::pow; using std::log;
1624  return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val()));
1625  }
1626 
1627  protected:
1628 
1630  const ExprT2& expr2;
1631  };
1632 
1633  //
1634  // Specialization for nested derivatives. This version does not use
1635  // if_then_else/ternary-operator on the base so that nested derivatives
1636  // are correct.
1637  //
1638  template <typename ExprT1, typename ExprT2>
1639  class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Nested >,
1640  ExprSpecDefault > {
1641 
1642  public:
1643 
1644  typedef typename ExprT1::value_type value_type_1;
1645  typedef typename ExprT2::value_type value_type_2;
1646  typedef typename Sacado::Promote<value_type_1,
1648 
1649  typedef typename ExprT1::scalar_type scalar_type_1;
1650  typedef typename ExprT2::scalar_type scalar_type_2;
1651  typedef typename Sacado::Promote<scalar_type_1,
1653 
1654  typedef typename ExprT1::base_expr_type base_expr_type_1;
1655  typedef typename ExprT2::base_expr_type base_expr_type_2;
1656  typedef typename Sacado::Promote<base_expr_type_1,
1658 
1660  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1661  expr1(expr1_), expr2(expr2_) {}
1662 
1664  int size() const {
1665  int sz1 = expr1.size(), sz2 = expr2.size();
1666  return sz1 > sz2 ? sz1 : sz2;
1667  }
1668 
1670  bool hasFastAccess() const {
1671  return expr1.hasFastAccess() && expr2.hasFastAccess();
1672  }
1673 
1675  bool isPassive() const {
1676  return expr1.isPassive() && expr2.isPassive();
1677  }
1678 
1680  bool updateValue() const {
1681  return expr1.updateValue() && expr2.updateValue();
1682  }
1683 
1685  void cache() const {}
1686 
1688  const value_type val() const {
1689  using std::pow;
1690  return pow(expr1.val(), expr2.val());
1691  }
1692 
1694  const value_type dx(int i) const {
1695  using std::pow; using std::log;
1696  const int sz1 = expr1.size(), sz2 = expr2.size();
1697  if (sz1 > 0 && sz2 > 0)
1698  return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1699  else if (sz1 > 0)
1700  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)));
1701  else
1702  return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1703  }
1704 
1706  const value_type fastAccessDx(int i) const {
1707  using std::pow; using std::log;
1708  return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1709  }
1710 
1711  protected:
1712 
1713  const ExprT1& expr1;
1714  const ExprT2& expr2;
1715 
1716  };
1717 
1718  template <typename ExprT1, typename T2>
1719  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Nested >,
1720  ExprSpecDefault > {
1721 
1722  public:
1723 
1726  typedef typename ExprT1::value_type value_type_1;
1728  typedef typename Sacado::Promote<value_type_1,
1730 
1731  typedef typename ExprT1::scalar_type scalar_type_1;
1733  typedef typename Sacado::Promote<scalar_type_1,
1735 
1736  typedef typename ExprT1::base_expr_type base_expr_type_1;
1738  typedef typename Sacado::Promote<base_expr_type_1,
1740 
1742  Expr(const ExprT1& expr1_, const ConstT& c_) :
1743  expr1(expr1_), c(c_) {}
1744 
1746  int size() const {
1747  return expr1.size();
1748  }
1749 
1751  bool hasFastAccess() const {
1752  return expr1.hasFastAccess();
1753  }
1754 
1756  bool isPassive() const {
1757  return expr1.isPassive();
1758  }
1759 
1761  bool updateValue() const { return expr1.updateValue(); }
1762 
1764  void cache() const {}
1765 
1767  const value_type val() const {
1768  using std::pow;
1769  return pow(expr1.val(), c.val());
1770  }
1771 
1773  const value_type dx(int i) const {
1774  using std::pow;
1775  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)));
1776  }
1777 
1779  const value_type fastAccessDx(int i) const {
1780  using std::pow;
1781  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)));
1782  }
1783 
1784  protected:
1785 
1786  const ExprT1& expr1;
1788  };
1789 
1790  template <typename T1, typename ExprT2>
1791  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Nested >,
1792  ExprSpecDefault > {
1793 
1794  public:
1795 
1799  typedef typename ExprT2::value_type value_type_2;
1800  typedef typename Sacado::Promote<value_type_1,
1802 
1804  typedef typename ExprT2::scalar_type scalar_type_2;
1805  typedef typename Sacado::Promote<scalar_type_1,
1807 
1809  typedef typename ExprT2::base_expr_type base_expr_type_2;
1810  typedef typename Sacado::Promote<base_expr_type_1,
1812 
1813 
1815  Expr(const ConstT& c_, const ExprT2& expr2_) :
1816  c(c_), expr2(expr2_) {}
1817 
1819  int size() const {
1820  return expr2.size();
1821  }
1822 
1824  bool hasFastAccess() const {
1825  return expr2.hasFastAccess();
1826  }
1827 
1829  bool isPassive() const {
1830  return expr2.isPassive();
1831  }
1832 
1834  bool updateValue() const { return expr2.updateValue(); }
1835 
1837  void cache() const {}
1838 
1840  const value_type val() const {
1841  using std::pow;
1842  return pow(c.val(), expr2.val());
1843  }
1844 
1846  const value_type dx(int i) const {
1847  using std::pow; using std::log;
1848  return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
1849  }
1850 
1852  const value_type fastAccessDx(int i) const {
1853  using std::pow; using std::log;
1854  return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
1855  }
1856 
1857  protected:
1858 
1860  const ExprT2& expr2;
1861  };
1862 
1863  //
1864  // Specialization for nested derivatives. This version does not use
1865  // if_then_else/ternary-operator on the base so that nested derivatives
1866  // are correct.
1867  //
1868  template <typename ExprT1, typename ExprT2>
1869  class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::NestedSimd >,
1870  ExprSpecDefault > {
1871 
1872  public:
1873 
1874  typedef typename ExprT1::value_type value_type_1;
1875  typedef typename ExprT2::value_type value_type_2;
1876  typedef typename Sacado::Promote<value_type_1,
1878 
1879  typedef typename ExprT1::scalar_type scalar_type_1;
1880  typedef typename ExprT2::scalar_type scalar_type_2;
1881  typedef typename Sacado::Promote<scalar_type_1,
1883 
1884  typedef typename ExprT1::base_expr_type base_expr_type_1;
1885  typedef typename ExprT2::base_expr_type base_expr_type_2;
1886  typedef typename Sacado::Promote<base_expr_type_1,
1888 
1890  Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1891  expr1(expr1_), expr2(expr2_) {}
1892 
1894  int size() const {
1895  int sz1 = expr1.size(), sz2 = expr2.size();
1896  return sz1 > sz2 ? sz1 : sz2;
1897  }
1898 
1900  bool hasFastAccess() const {
1901  return expr1.hasFastAccess() && expr2.hasFastAccess();
1902  }
1903 
1905  bool isPassive() const {
1906  return expr1.isPassive() && expr2.isPassive();
1907  }
1908 
1910  bool updateValue() const {
1911  return expr1.updateValue() && expr2.updateValue();
1912  }
1913 
1915  void cache() const {}
1916 
1918  const value_type val() const {
1919  using std::pow;
1920  return pow(expr1.val(), expr2.val());
1921  }
1922 
1924  const value_type dx(int i) const {
1925  using std::pow; using std::log; using Sacado::if_then_else;
1926  const int sz1 = expr1.size(), sz2 = expr2.size();
1927  if (sz1 > 0 && sz2 > 0)
1928  return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1929  else if (sz1 > 0)
1930  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))));
1931  else
1932  return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1933  }
1934 
1936  const value_type fastAccessDx(int i) const {
1937  using std::pow; using std::log;
1938  return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1939  }
1940 
1941  protected:
1942 
1943  const ExprT1& expr1;
1944  const ExprT2& expr2;
1945 
1946  };
1947 
1948  template <typename ExprT1, typename T2>
1949  class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::NestedSimd >,
1950  ExprSpecDefault > {
1951 
1952  public:
1953 
1956  typedef typename ExprT1::value_type value_type_1;
1958  typedef typename Sacado::Promote<value_type_1,
1960 
1961  typedef typename ExprT1::scalar_type scalar_type_1;
1963  typedef typename Sacado::Promote<scalar_type_1,
1965 
1966  typedef typename ExprT1::base_expr_type base_expr_type_1;
1968  typedef typename Sacado::Promote<base_expr_type_1,
1970 
1972  Expr(const ExprT1& expr1_, const ConstT& c_) :
1973  expr1(expr1_), c(c_) {}
1974 
1976  int size() const {
1977  return expr1.size();
1978  }
1979 
1981  bool hasFastAccess() const {
1982  return expr1.hasFastAccess();
1983  }
1984 
1986  bool isPassive() const {
1987  return expr1.isPassive();
1988  }
1989 
1991  bool updateValue() const { return expr1.updateValue(); }
1992 
1994  void cache() const {}
1995 
1997  const value_type val() const {
1998  using std::pow;
1999  return pow(expr1.val(), c.val());
2000  }
2001 
2003  const value_type dx(int i) const {
2004  using std::pow; using Sacado::if_then_else;
2005  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))));
2006  }
2007 
2009  const value_type fastAccessDx(int i) const {
2010  using std::pow; using Sacado::if_then_else;
2011  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))));
2012  }
2013 
2014  protected:
2015 
2016  const ExprT1& expr1;
2018  };
2019 
2020  template <typename T1, typename ExprT2>
2021  class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::NestedSimd >,
2022  ExprSpecDefault > {
2023 
2024  public:
2025 
2029  typedef typename ExprT2::value_type value_type_2;
2030  typedef typename Sacado::Promote<value_type_1,
2032 
2034  typedef typename ExprT2::scalar_type scalar_type_2;
2035  typedef typename Sacado::Promote<scalar_type_1,
2037 
2039  typedef typename ExprT2::base_expr_type base_expr_type_2;
2040  typedef typename Sacado::Promote<base_expr_type_1,
2042 
2043 
2045  Expr(const ConstT& c_, const ExprT2& expr2_) :
2046  c(c_), expr2(expr2_) {}
2047 
2049  int size() const {
2050  return expr2.size();
2051  }
2052 
2054  bool hasFastAccess() const {
2055  return expr2.hasFastAccess();
2056  }
2057 
2059  bool isPassive() const {
2060  return expr2.isPassive();
2061  }
2062 
2064  bool updateValue() const { return expr2.updateValue(); }
2065 
2067  void cache() const {}
2068 
2070  const value_type val() const {
2071  using std::pow;
2072  return pow(c.val(), expr2.val());
2073  }
2074 
2076  const value_type dx(int i) const {
2077  using std::pow; using std::log;
2078  return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
2079  }
2080 
2082  const value_type fastAccessDx(int i) const {
2083  using std::pow; using std::log;
2084  return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
2085  }
2086 
2087  protected:
2088 
2090  const ExprT2& expr2;
2091  };
2092 
2093  template <typename T1, typename T2>
2095  typename mpl::enable_if_c<
2098  >::type
2099  /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(PowerOp)*/
2100  pow (const Expr<T1>& expr1, const Expr<T2>& expr2)
2101  {
2102  typedef PowerOp< Expr<T1>, Expr<T2> > expr_t;
2103 
2104  return Expr<expr_t>(expr1, expr2);
2105  }
2106 
2107  template <typename T>
2109  Expr< PowerOp< Expr<T>, Expr<T> > >
2110  pow (const Expr<T>& expr1, const Expr<T>& expr2)
2111  {
2112  typedef PowerOp< Expr<T>, Expr<T> > expr_t;
2113 
2114  return Expr<expr_t>(expr1, expr2);
2115  }
2116 
2117  template <typename T>
2120  pow (const typename Expr<T>::value_type& c,
2121  const Expr<T>& expr)
2122  {
2124  typedef PowerOp< ConstT, Expr<T> > expr_t;
2125 
2126  return Expr<expr_t>(ConstT(c), expr);
2127  }
2128 
2129  template <typename T>
2131  Expr< PowerOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
2132  pow (const Expr<T>& expr,
2133  const typename Expr<T>::value_type& c)
2134  {
2136  typedef PowerOp< Expr<T>, ConstT > expr_t;
2137 
2138  return Expr<expr_t>(expr, ConstT(c));
2139  }
2140 
2141  template <typename T>
2144  pow (const typename Expr<T>::scalar_type& c,
2145  const Expr<T>& expr)
2146  {
2148  typedef PowerOp< ConstT, Expr<T> > expr_t;
2149 
2150  return Expr<expr_t>(ConstT(c), expr);
2151  }
2152 
2153  template <typename T>
2156  pow (const Expr<T>& expr,
2157  const typename Expr<T>::scalar_type& c)
2158  {
2160  typedef PowerOp< Expr<T>, ConstT > expr_t;
2161 
2162  return Expr<expr_t>(expr, ConstT(c));
2163  }
2164  }
2165 }
2166 
2167 //--------------------------if_then_else operator -----------------------
2168 // Can't use the above macros because it is a ternary operator (sort of).
2169 // Also, relies on C++11
2170 
2171 #ifdef HAVE_SACADO_CXX11
2172 
2173 namespace Sacado {
2174  namespace Fad {
2175 
2176  template <typename CondT, typename ExprT1, typename ExprT2>
2177  class IfThenElseOp {};
2178 
2179  template <typename CondT, typename ExprT1, typename ExprT2>
2180  struct ExprSpec< IfThenElseOp< CondT, ExprT1, ExprT2 > > {
2181  typedef typename ExprSpec<ExprT1>::type type;
2182  };
2183 
2184  template <typename CondT, typename ExprT1, typename ExprT2>
2185  class Expr< IfThenElseOp< CondT, ExprT1, ExprT2 >,ExprSpecDefault > {
2186 
2187  public:
2188 
2189  typedef typename ExprT1::value_type value_type_1;
2190  typedef typename ExprT2::value_type value_type_2;
2191  typedef typename Sacado::Promote<value_type_1,
2192  value_type_2>::type value_type;
2193 
2194  typedef typename ExprT1::scalar_type scalar_type_1;
2195  typedef typename ExprT2::scalar_type scalar_type_2;
2196  typedef typename Sacado::Promote<scalar_type_1,
2197  scalar_type_2>::type scalar_type;
2198 
2199  typedef typename ExprT1::base_expr_type base_expr_type_1;
2200  typedef typename ExprT2::base_expr_type base_expr_type_2;
2201  typedef typename Sacado::Promote<base_expr_type_1,
2202  base_expr_type_2>::type base_expr_type;
2203 
2205  Expr(const CondT& cond_, const ExprT1& expr1_, const ExprT2& expr2_) :
2206  cond(cond_), expr1(expr1_), expr2(expr2_) {}
2207 
2209  int size() const {
2210  int sz1 = expr1.size(), sz2 = expr2.size();
2211  return sz1 > sz2 ? sz1 : sz2;
2212  }
2213 
2215  bool hasFastAccess() const {
2216  return expr1.hasFastAccess() && expr2.hasFastAccess();
2217  }
2218 
2220  bool isPassive() const {
2221  return expr1.isPassive() && expr2.isPassive();
2222  }
2223 
2225  bool updateValue() const {
2226  return expr1.updateValue() && expr2.updateValue();
2227  }
2228 
2230  void cache() const {}
2231 
2233  const value_type val() const {
2234  using Sacado::if_then_else;
2235  return if_then_else( cond, expr1.val(), expr2.val() );
2236  }
2237 
2239  const value_type dx(int i) const {
2240  using Sacado::if_then_else;
2241  return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
2242  }
2243 
2245  const value_type fastAccessDx(int i) const {
2246  using Sacado::if_then_else;
2247  return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
2248  }
2249 
2250  protected:
2251 
2252  const CondT& cond;
2253  const ExprT1& expr1;
2254  const ExprT2& expr2;
2255 
2256  };
2257 
2258  template <typename CondT, typename ExprT1, typename T2>
2259  struct ExprSpec< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> > > {
2260  typedef typename ExprSpec<ExprT1>::type type;
2261  };
2262 
2263  template <typename CondT, typename ExprT1, typename T2>
2264  class Expr< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
2265 
2266  public:
2267 
2268  typedef ConstExpr<T2> ConstT;
2269  typedef ConstExpr<T2> ExprT2;
2270  typedef typename ExprT1::value_type value_type_1;
2271  typedef typename ExprT2::value_type value_type_2;
2272  typedef typename Sacado::Promote<value_type_1,
2273  value_type_2>::type value_type;
2274 
2275  typedef typename ExprT1::scalar_type scalar_type_1;
2276  typedef typename ExprT2::scalar_type scalar_type_2;
2277  typedef typename Sacado::Promote<scalar_type_1,
2278  scalar_type_2>::type scalar_type;
2279 
2280  typedef typename ExprT1::base_expr_type base_expr_type_1;
2281  typedef typename ExprT2::base_expr_type base_expr_type_2;
2282  typedef typename Sacado::Promote<base_expr_type_1,
2283  base_expr_type_2>::type base_expr_type;
2284 
2286  Expr(const CondT& cond_, const ExprT1& expr1_, const ConstT& c_) :
2287  cond(cond_), expr1(expr1_), c(c_) {}
2288 
2290  int size() const {
2291  return expr1.size();
2292  }
2293 
2295  bool hasFastAccess() const {
2296  return expr1.hasFastAccess();
2297  }
2298 
2300  bool isPassive() const {
2301  return expr1.isPassive();
2302  }
2303 
2305  bool updateValue() const { return expr1.updateValue(); }
2306 
2308  void cache() const {}
2309 
2311  const value_type val() const {
2312  using Sacado::if_then_else;
2313  return if_then_else( cond, expr1.val(), c.val() );
2314  }
2315 
2317  const value_type dx(int i) const {
2318  using Sacado::if_then_else;
2319  return if_then_else( cond, expr1.dx(i), value_type(0.0) );
2320  }
2321 
2323  const value_type fastAccessDx(int i) const {
2324  using Sacado::if_then_else;
2325  return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
2326  }
2327 
2328  protected:
2329 
2330  const CondT& cond;
2331  const ExprT1& expr1;
2332  ConstT c;
2333  };
2334 
2335  template <typename CondT, typename T1, typename ExprT2>
2336  struct ExprSpec< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 > > {
2337  typedef typename ExprSpec<ExprT2>::type type;
2338  };
2339 
2340  template <typename CondT, typename T1, typename ExprT2>
2341  class Expr< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
2342 
2343  public:
2344 
2345  typedef ConstExpr<T1> ConstT;
2346  typedef ConstExpr<T1> ExprT1;
2347  typedef typename ExprT1::value_type value_type_1;
2348  typedef typename ExprT2::value_type value_type_2;
2349  typedef typename Sacado::Promote<value_type_1,
2350  value_type_2>::type value_type;
2351 
2352  typedef typename ExprT1::scalar_type scalar_type_1;
2353  typedef typename ExprT2::scalar_type scalar_type_2;
2354  typedef typename Sacado::Promote<scalar_type_1,
2355  scalar_type_2>::type scalar_type;
2356 
2357  typedef typename ExprT1::base_expr_type base_expr_type_1;
2358  typedef typename ExprT2::base_expr_type base_expr_type_2;
2359  typedef typename Sacado::Promote<base_expr_type_1,
2360  base_expr_type_2>::type base_expr_type;
2361 
2363  Expr(const CondT& cond_, const ConstT& c_, const ExprT2& expr2_) :
2364  cond(cond_), c(c_), expr2(expr2_) {}
2365 
2367  int size() const {
2368  return expr2.size();
2369  }
2370 
2372  bool hasFastAccess() const {
2373  return expr2.hasFastAccess();
2374  }
2375 
2377  bool isPassive() const {
2378  return expr2.isPassive();
2379  }
2380 
2382  bool updateValue() const { return expr2.updateValue(); }
2383 
2385  void cache() const {}
2386 
2388  const value_type val() const {
2389  using Sacado::if_then_else;
2390  return if_then_else( cond, c.val(), expr2.val() );
2391  }
2392 
2394  const value_type dx(int i) const {
2395  using Sacado::if_then_else;
2396  return if_then_else( cond, value_type(0.0), expr2.dx(i) );
2397  }
2398 
2400  const value_type fastAccessDx(int i) const {
2401  using Sacado::if_then_else;
2402  return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
2403  }
2404 
2405  protected:
2406 
2407  const CondT& cond;
2408  ConstT c;
2409  const ExprT2& expr2;
2410  };
2411 
2412  template <typename CondT, typename T1, typename T2>
2416  Expr< IfThenElseOp< CondT, T1, T2 > >
2417  >::type
2418  if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
2419  {
2420  typedef IfThenElseOp< CondT, T1, T2 > expr_t;
2421 
2422  return Expr<expr_t>(cond, expr1, expr2);
2423  }
2424 
2425  template <typename CondT, typename T>
2427  Expr< IfThenElseOp< CondT, Expr<T>, Expr<T> > >
2428  if_then_else (const CondT& cond, const Expr<T>& expr1, const Expr<T>& expr2)
2429  {
2430  typedef IfThenElseOp< CondT, Expr<T>, Expr<T> > expr_t;
2431 
2432  return Expr<expr_t>(cond, expr1, expr2);
2433  }
2434 
2435  template <typename CondT, typename T>
2438  Expr<T> > >
2439  if_then_else (const CondT& cond, const typename Expr<T>::value_type& c,
2440  const Expr<T>& expr)
2441  {
2443  typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2444 
2445  return Expr<expr_t>(cond, ConstT(c), expr);
2446  }
2447 
2448  template <typename CondT, typename T>
2450  Expr< IfThenElseOp< CondT, Expr<T>,
2452  if_then_else (const CondT& cond, const Expr<T>& expr,
2453  const typename Expr<T>::value_type& c)
2454  {
2456  typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2457 
2458  return Expr<expr_t>(cond, expr, ConstT(c));
2459  }
2460 
2461  template <typename CondT, typename T>
2463  typename mpl::disable_if<
2464  mpl::is_same< typename Expr<T>::value_type,
2465  typename Expr<T>::scalar_type>,
2466  Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::scalar_type>,
2467  Expr<T> > >
2468  >::type
2469  if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
2470  const Expr<T>& expr)
2471  {
2473  typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2474 
2475  return Expr<expr_t>(cond, ConstT(c), expr);
2476  }
2477 
2478  template <typename CondT, typename T>
2480  typename mpl::disable_if<
2481  mpl::is_same< typename Expr<T>::value_type,
2482  typename Expr<T>::scalar_type>,
2483  Expr< IfThenElseOp< CondT, Expr<T>,
2485  >::type
2486  if_then_else (const CondT& cond, const Expr<T>& expr,
2487  const typename Expr<T>::scalar_type& c)
2488  {
2490  typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2491 
2492  return Expr<expr_t>(cond, expr, ConstT(c));
2493  }
2494  }
2495 }
2496 
2497 #endif
2498 
2499 //-------------------------- Relational Operators -----------------------
2500 
2501 #ifdef HAVE_SACADO_CXX11
2502 
2503 namespace Sacado {
2504  namespace Fad {
2505  template <typename T1, typename T2 = T1>
2506  struct ConditionalReturnType {
2507  typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
2508  };
2509  }
2510 }
2511 
2512 #define FAD_RELOP_MACRO(OP) \
2513 namespace Sacado { \
2514  namespace Fad { \
2515  template <typename ExprT1, typename ExprT2> \
2516  SACADO_INLINE_FUNCTION \
2517  typename ConditionalReturnType<typename Expr<ExprT1>::value_type, \
2518  typename Expr<ExprT2>::value_type>::type \
2519  operator OP (const Expr<ExprT1>& expr1, \
2520  const Expr<ExprT2>& expr2) \
2521  { \
2522  return expr1.val() OP expr2.val(); \
2523  } \
2524  \
2525  template <typename ExprT2> \
2526  SACADO_INLINE_FUNCTION \
2527  typename ConditionalReturnType<typename Expr<ExprT2>::value_type>::type \
2528  operator OP (const typename Expr<ExprT2>::value_type& a, \
2529  const Expr<ExprT2>& expr2) \
2530  { \
2531  return a OP expr2.val(); \
2532  } \
2533  \
2534  template <typename ExprT1> \
2535  SACADO_INLINE_FUNCTION \
2536  typename ConditionalReturnType<typename Expr<ExprT1>::value_type>::type \
2537  operator OP (const Expr<ExprT1>& expr1, \
2538  const typename Expr<ExprT1>::value_type& b) \
2539  { \
2540  return expr1.val() OP b; \
2541  } \
2542  } \
2543 }
2544 
2545 #else
2546 
2547 #define FAD_RELOP_MACRO(OP) \
2548 namespace Sacado { \
2549  namespace Fad { \
2550  template <typename ExprT1, typename ExprT2> \
2551  SACADO_INLINE_FUNCTION \
2552  bool \
2553  operator OP (const Expr<ExprT1>& expr1, \
2554  const Expr<ExprT2>& expr2) \
2555  { \
2556  return expr1.val() OP expr2.val(); \
2557  } \
2558  \
2559  template <typename ExprT2> \
2560  SACADO_INLINE_FUNCTION \
2561  bool \
2562  operator OP (const typename Expr<ExprT2>::value_type& a, \
2563  const Expr<ExprT2>& expr2) \
2564  { \
2565  return a OP expr2.val(); \
2566  } \
2567  \
2568  template <typename ExprT1> \
2569  SACADO_INLINE_FUNCTION \
2570  bool \
2571  operator OP (const Expr<ExprT1>& expr1, \
2572  const typename Expr<ExprT1>::value_type& b) \
2573  { \
2574  return expr1.val() OP b; \
2575  } \
2576  } \
2577 }
2578 
2579 #endif
2580 
2581 FAD_RELOP_MACRO(==)
2582 FAD_RELOP_MACRO(!=)
2583 FAD_RELOP_MACRO(<)
2584 FAD_RELOP_MACRO(>)
2585 FAD_RELOP_MACRO(<=)
2586 FAD_RELOP_MACRO(>=)
2587 FAD_RELOP_MACRO(<<=)
2588 FAD_RELOP_MACRO(>>=)
2589 FAD_RELOP_MACRO(&)
2590 FAD_RELOP_MACRO(|)
2591 
2592 #undef FAD_RELOP_MACRO
2593 
2594 namespace Sacado {
2595 
2596  namespace Fad {
2597 
2598  template <typename ExprT>
2600  bool operator ! (const Expr<ExprT>& expr)
2601  {
2602  return ! expr.val();
2603  }
2604 
2605  } // namespace Fad
2606 
2607 } // namespace Sacado
2608 
2609 //-------------------------- Boolean Operators -----------------------
2610 namespace Sacado {
2611 
2612  namespace Fad {
2613 
2614  template <typename ExprT>
2616  bool toBool(const Expr<ExprT>& x) {
2617  bool is_zero = (x.val() == 0.0);
2618  for (int i=0; i<x.size(); i++)
2619  is_zero = is_zero && (x.dx(i) == 0.0);
2620  return !is_zero;
2621  }
2622 
2623  } // namespace Fad
2624 
2625 } // namespace Sacado
2626 
2627 #define FAD_BOOL_MACRO(OP) \
2628 namespace Sacado { \
2629  namespace Fad { \
2630  template <typename ExprT1, typename ExprT2> \
2631  SACADO_INLINE_FUNCTION \
2632  bool \
2633  operator OP (const Expr<ExprT1>& expr1, \
2634  const Expr<ExprT2>& expr2) \
2635  { \
2636  return toBool(expr1) OP toBool(expr2); \
2637  } \
2638  \
2639  template <typename ExprT2> \
2640  SACADO_INLINE_FUNCTION \
2641  bool \
2642  operator OP (const typename Expr<ExprT2>::value_type& a, \
2643  const Expr<ExprT2>& expr2) \
2644  { \
2645  return a OP toBool(expr2); \
2646  } \
2647  \
2648  template <typename ExprT1> \
2649  SACADO_INLINE_FUNCTION \
2650  bool \
2651  operator OP (const Expr<ExprT1>& expr1, \
2652  const typename Expr<ExprT1>::value_type& b) \
2653  { \
2654  return toBool(expr1) OP b; \
2655  } \
2656  } \
2657 }
2658 
2659 FAD_BOOL_MACRO(&&)
2660 FAD_BOOL_MACRO(||)
2661 
2662 #undef FAD_BOOL_MACRO
2663 
2664 //-------------------------- I/O Operators -----------------------
2665 
2666 namespace Sacado {
2667 
2668  namespace Fad {
2669 
2670  template <typename ExprT>
2671  std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
2672  os << x.val() << " [";
2673 
2674  for (int i=0; i< x.size(); i++) {
2675  os << " " << x.dx(i);
2676  }
2677 
2678  os << " ]";
2679  return os;
2680  }
2681 
2682  } // namespace Fad
2683 
2684 } // namespace Sacado
2685 
2686 #endif // SACADO_FAD_OPS_HPP
Wrapper for a generic expression template.
cbrt(expr.val())
expr expr ASinhOp
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from ConstT)
Constant expression template.
atan2(expr1.val(), expr2.val())
char c_
expr val()
expr expr TanOp
#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)
tanh(expr.val())
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
sqrt(expr.val())
abs(expr.val())
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
expr expr SinOp
expr expr ASinOp
expr expr SinhOp
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ConstT &c_)
#define SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP)
fabs(expr.val())
expr expr Log10Op
expr true
expr expr SqrtOp
exp(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
SimpleFad< ValueT > log(const SimpleFad< ValueT > &a)
expr expr CosOp
log(expr.val())
#define T
Definition: Sacado_rad.hpp:573
SACADO_INLINE_FUNCTION bool toBool(const Expr< ExprT > &x)
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
expr expr AbsOp
#define FAD_BOOL_MACRO(OP)
#define T2(r, f)
Definition: Sacado_rad.hpp:578
atan(expr.val())
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 expr ExpOp
atanh(expr.val())
asin(expr.val())
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
log10(expr.val())
tan(expr.val())
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
expr expr ACosOp
sin(expr.val())
UnaryMinusOp
T2 base_expr_type
Typename of base-expressions.
SACADO_INLINE_FUNCTION pow(const Expr< T > &expr, const typename Expr< T >::scalar_type &c)
#define T1(r, f)
Definition: Sacado_rad.hpp:603
cosh(expr.val())
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
Meta-function for determining nesting with an expression.
expr expr ATanOp
SubtractionOp
int value
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
cos(expr.val())
#define SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP)
SACADO_INLINE_FUNCTION T safe_sqrt(const T &x)
static const unsigned value
expr expr dx(i)
expr1 expr1 expr1 c expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr1 expr2 expr1 expr2 expr1 Atan2Op
SACADO_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
sinh(expr.val())
acosh(expr.val())
expr expr ATanhOp
SimpleFad< ValueT > operator*(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
acos(expr.val())
SACADO_INLINE_FUNCTION Expr(const ConstT &c_, const ExprT2 &expr2_)
SimpleFad< ValueT > max(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
ConstT 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)
if_then_else(expr.val() >=0, expr.dx(i), value_type(-expr.dx(i)))
asinh(expr.val())
expr expr ACoshOp
expr expr expr bar false
#define SACADO_INLINE_FUNCTION
expr2 expr2 c c c c MaxOp
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
Base template specification for Promote.
expr1 expr1 expr1 c expr1 expr2 expr1 expr2 expr1 DivisionOp
expr expr CoshOp
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
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
expr expr TanhOp