Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_ELRFad_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_ELRFAD_OPS_HPP
34 #define SACADO_ELRFAD_OPS_HPP
35 
37 #include "Sacado_cmath.hpp"
38 #include <ostream> // for std::ostream
39 
40 #define FAD_UNARYOP_MACRO(OPNAME,OP,VALUE,ADJOINT, \
41  LINEAR,DX,FASTACCESSDX) \
42 namespace Sacado { \
43  namespace ELRFad { \
44  \
45  template <typename ExprT> \
46  class OP {}; \
47  \
48  template <typename ExprT> \
49  class Expr< OP<ExprT> > { \
50  public: \
51  \
52  typedef typename ExprT::value_type value_type; \
53  typedef typename ExprT::scalar_type scalar_type; \
54  typedef typename ExprT::base_expr_type base_expr_type; \
55  \
56  static const int num_args = ExprT::num_args; \
57  \
58  static const bool is_linear = LINEAR; \
59  \
60  SACADO_INLINE_FUNCTION \
61  explicit Expr(const ExprT& expr_) : expr(expr_) {} \
62  \
63  SACADO_INLINE_FUNCTION \
64  int size() const { return expr.size(); } \
65  \
66  template <int Arg> \
67  SACADO_INLINE_FUNCTION \
68  bool isActive() const { return expr.template isActive<Arg>(); } \
69  \
70  SACADO_INLINE_FUNCTION \
71  bool isActive2(int j) const { return expr.isActive2(j); } \
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  return VALUE; \
82  } \
83  \
84  SACADO_INLINE_FUNCTION \
85  void computePartials(const value_type& bar, \
86  value_type partials[]) const { \
87  expr.computePartials(ADJOINT, partials); \
88  } \
89  \
90  SACADO_INLINE_FUNCTION \
91  void getTangents(int i, value_type dots[]) const { \
92  expr.getTangents(i, dots); } \
93  \
94  template <int Arg> \
95  SACADO_INLINE_FUNCTION \
96  const value_type& getTangent(int i) const { \
97  return expr.template getTangent<Arg>(i); \
98  } \
99  \
100  SACADO_INLINE_FUNCTION \
101  bool isLinear() const { \
102  return LINEAR; \
103  } \
104  \
105  SACADO_INLINE_FUNCTION \
106  bool hasFastAccess() const { \
107  return expr.hasFastAccess(); \
108  } \
109  \
110  SACADO_INLINE_FUNCTION \
111  const value_type dx(int i) const { \
112  return DX; \
113  } \
114  \
115  SACADO_INLINE_FUNCTION \
116  const value_type fastAccessDx(int i) const { \
117  return FASTACCESSDX; \
118  } \
119  \
120  SACADO_INLINE_FUNCTION \
121  const value_type* getDx(int j) const { \
122  return expr.getDx(j); \
123  } \
124  \
125  SACADO_INLINE_FUNCTION \
126  int numActiveArgs() const { \
127  return expr.numActiveArgs(); \
128  } \
129  \
130  SACADO_INLINE_FUNCTION \
131  void computeActivePartials(const value_type& bar, \
132  value_type *partials) const { \
133  expr.computePartials(ADJOINT, partials); \
134  } \
135  \
136  protected: \
137  \
138  const ExprT& expr; \
139  }; \
140  \
141  template <typename T> \
142  SACADO_INLINE_FUNCTION \
143  Expr< OP< Expr<T> > > \
144  OPNAME (const Expr<T>& expr) \
145  { \
146  typedef OP< Expr<T> > expr_t; \
147  \
148  return Expr<expr_t>(expr); \
149  } \
150  } \
151 }
152 
153 FAD_UNARYOP_MACRO(operator+,
154  UnaryPlusOp,
155  expr.val(),
156  bar,
157  true,
158  expr.dx(i),
159  expr.fastAccessDx(i))
160 FAD_UNARYOP_MACRO(operator-,
162  -expr.val(),
163  -bar,
164  true,
165  -expr.dx(i),
166  -expr.fastAccessDx(i))
169  std::exp(expr.val()),
170  bar*std::exp(expr.val()),
171  false,
172  std::exp(expr.val())*expr.dx(i),
173  std::exp(expr.val())*expr.fastAccessDx(i))
176  std::log(expr.val()),
177  bar/expr.val(),
178  false,
179  expr.dx(i)/expr.val(),
180  expr.fastAccessDx(i)/expr.val())
183  std::log10(expr.val()),
184  bar/( std::log(value_type(10.))*expr.val() ),
185  false,
186  expr.dx(i)/( std::log(value_type(10))*expr.val()),
187  expr.fastAccessDx(i) / ( std::log(value_type(10))*expr.val()))
189  SqrtOp,
190  std::sqrt(expr.val()),
191  value_type(0.5)*bar/std::sqrt(expr.val()),
192  false,
193  expr.dx(i)/(value_type(2)* std::sqrt(expr.val())),
194  expr.fastAccessDx(i)/(value_type(2)* std::sqrt(expr.val())))
196  SafeSqrtOp,
197  std::sqrt(expr.val()),
198  expr.val() == value_type(0.0) ? value_type(0.0) : value_type(value_type(0.5)*bar/std::sqrt(expr.val())),
199  false,
200  expr.val() == value_type(0.0) ? value_type(0.0) : value_type(expr.dx(i)/(value_type(2)*std::sqrt(expr.val()))),
201  expr.val() == value_type(0.0) ? value_type(0.0) : value_type(expr.fastAccessDx(i)/(value_type(2)*std::sqrt(expr.val()))))
203  CosOp,
204  std::cos(expr.val()),
205  -bar*std::sin(expr.val()),
206  false,
207  -expr.dx(i)* std::sin(expr.val()),
208  -expr.fastAccessDx(i)* std::sin(expr.val()))
210  SinOp,
211  std::sin(expr.val()),
212  bar*std::cos(expr.val()),
213  false,
214  expr.dx(i)* std::cos(expr.val()),
215  expr.fastAccessDx(i)* std::cos(expr.val()))
217  TanOp,
218  std::tan(expr.val()),
219  bar*(value_type(1.)+ std::tan(expr.val())*std::tan(expr.val())),
220  false,
221  expr.dx(i)*
222  (value_type(1)+ std::tan(expr.val())* std::tan(expr.val())),
223  expr.fastAccessDx(i)*
224  (value_type(1)+ std::tan(expr.val())* std::tan(expr.val())))
226  ACosOp,
227  std::acos(expr.val()),
228  -bar/std::sqrt(value_type(1.)-expr.val()*expr.val()),
229  false,
230  -expr.dx(i)/ std::sqrt(value_type(1)-expr.val()*expr.val()),
231  -expr.fastAccessDx(i) /
232  std::sqrt(value_type(1)-expr.val()*expr.val()))
234  ASinOp,
235  std::asin(expr.val()),
236  bar/std::sqrt(value_type(1.)-expr.val()*expr.val()),
237  false,
238  expr.dx(i)/ std::sqrt(value_type(1)-expr.val()*expr.val()),
239  expr.fastAccessDx(i) /
240  std::sqrt(value_type(1)-expr.val()*expr.val()))
242  ATanOp,
243  std::atan(expr.val()),
244  bar/(value_type(1.)+expr.val()*expr.val()),
245  false,
246  expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
247  expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
249  CoshOp,
250  std::cosh(expr.val()),
251  bar*std::sinh(expr.val()),
252  false,
253  expr.dx(i)* std::sinh(expr.val()),
254  expr.fastAccessDx(i)* std::sinh(expr.val()))
256  SinhOp,
257  std::sinh(expr.val()),
258  bar*std::cosh(expr.val()),
259  false,
260  expr.dx(i)* std::cosh(expr.val()),
261  expr.fastAccessDx(i)* std::cosh(expr.val()))
263  TanhOp,
264  std::tanh(expr.val()),
265  bar*(value_type(1)-std::tanh(expr.val())*std::tanh(expr.val())),
266  false,
267  expr.dx(i)*(value_type(1)-std::tanh(expr.val())*std::tanh(expr.val())),
268  expr.fastAccessDx(i)*(value_type(1)-std::tanh(expr.val())*std::tanh(expr.val())))
270  ACoshOp,
271  std::acosh(expr.val()),
272  bar/std::sqrt((expr.val()-value_type(1.)) *
273  (expr.val()+value_type(1.))),
274  false,
275  expr.dx(i)/ std::sqrt((expr.val()-value_type(1)) *
276  (expr.val()+value_type(1))),
277  expr.fastAccessDx(i)/ std::sqrt((expr.val()-value_type(1)) *
278  (expr.val()+value_type(1))))
280  ASinhOp,
281  std::asinh(expr.val()),
282  bar/std::sqrt(value_type(1.)+expr.val()*expr.val()),
283  false,
284  expr.dx(i)/ std::sqrt(value_type(1)+expr.val()*expr.val()),
285  expr.fastAccessDx(i)/ std::sqrt(value_type(1)+
286  expr.val()*expr.val()))
288  ATanhOp,
289  std::atanh(expr.val()),
290  bar/(value_type(1.)-expr.val()*expr.val()),
291  false,
292  expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
293  expr.fastAccessDx(i)/(value_type(1)-
294  expr.val()*expr.val()))
296  AbsOp,
297  std::abs(expr.val()),
298  (expr.val() >= value_type(0.)) ? bar : value_type(-bar),
299  false,
300  expr.val() >= 0 ? value_type(+expr.dx(i)) :
301  value_type(-expr.dx(i)),
302  expr.val() >= 0 ? value_type(+expr.fastAccessDx(i)) :
303  value_type(-expr.fastAccessDx(i)))
305  FAbsOp,
306  std::fabs(expr.val()),
307  (expr.val() >= value_type(0.)) ? bar : value_type(-bar),
308  false,
309  expr.val() >= 0 ? value_type(+expr.dx(i)) :
310  value_type(-expr.dx(i)),
311  expr.val() >= 0 ? value_type(+expr.fastAccessDx(i)) :
312  value_type(-expr.fastAccessDx(i)))
314  CbrtOp,
315  std::cbrt(expr.val()),
316  bar/(value_type(3)*std::cbrt(expr.val()*expr.val())),
317  false,
318  expr.dx(i)/(value_type(3)*std::cbrt(expr.val()*expr.val())),
319  expr.fastAccessDx(i)/(value_type(3)*std::cbrt(expr.val()*expr.val())))
320 
321 #undef FAD_UNARYOP_MACRO
322 
323 #define FAD_BINARYOP_MACRO( \
324  OPNAME,OP,VALUE,LADJOINT,RADJOINT, \
325  LINEAR,CONST_LINEAR_1, CONST_LINEAR_2, \
326  LINEAR_2,CONST_LINEAR_1_2, CONST_LINEAR_2_2, \
327  DX,FASTACCESSDX,CONST_DX_1,CONST_DX_2, \
328  CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
329 namespace Sacado { \
330  namespace ELRFad { \
331  \
332  template <typename ExprT1, typename ExprT2> \
333  class OP {}; \
334  \
335  template <typename ExprT1, typename ExprT2> \
336  class Expr< OP<ExprT1,ExprT2> > { \
337  \
338  public: \
339  \
340  typedef typename ExprT1::value_type value_type_1; \
341  typedef typename ExprT2::value_type value_type_2; \
342  typedef typename Sacado::Promote<value_type_1, \
343  value_type_2>::type value_type; \
344  \
345  typedef typename ExprT1::scalar_type scalar_type_1; \
346  typedef typename ExprT2::scalar_type scalar_type_2; \
347  typedef typename Sacado::Promote<scalar_type_1, \
348  scalar_type_2>::type scalar_type; \
349  \
350  typedef typename ExprT1::base_expr_type base_expr_type_1; \
351  typedef typename ExprT2::base_expr_type base_expr_type_2; \
352  typedef typename Sacado::Promote<base_expr_type_1, \
353  base_expr_type_2>::type base_expr_type; \
354  \
355  static const int num_args1 = ExprT1::num_args; \
356  static const int num_args2 = ExprT2::num_args; \
357  static const int num_args = num_args1 + num_args2; \
358  \
359  static const bool is_linear = LINEAR_2; \
360  \
361  SACADO_INLINE_FUNCTION \
362  Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
363  expr1(expr1_), expr2(expr2_) {} \
364  \
365  SACADO_INLINE_FUNCTION \
366  int size() const { \
367  int sz1 = expr1.size(), sz2 = expr2.size(); \
368  return sz1 > sz2 ? sz1 : sz2; \
369  } \
370  \
371  template <int Arg> \
372  SACADO_INLINE_FUNCTION \
373  bool isActive() const { \
374  if (Arg < num_args1) \
375  return expr1.template isActive<Arg>(); \
376  else \
377  return expr2.template isActive<Arg-num_args1>(); \
378  } \
379  \
380  SACADO_INLINE_FUNCTION \
381  bool isActive2(int j) const { \
382  if (j < num_args1) \
383  return expr1.isActive2(j); \
384  else \
385  return expr2.isActive2(j); \
386  } \
387  \
388  SACADO_INLINE_FUNCTION \
389  bool updateValue() const { \
390  return expr1.updateValue() && expr2.updateValue(); \
391  } \
392  \
393  SACADO_INLINE_FUNCTION \
394  void cache() const {} \
395  \
396  SACADO_INLINE_FUNCTION \
397  value_type val() const { \
398  return VALUE; \
399  } \
400  \
401  SACADO_INLINE_FUNCTION \
402  void computePartials(const value_type& bar, \
403  value_type partials[]) const { \
404  if (num_args1 > 0) \
405  expr1.computePartials(LADJOINT, partials); \
406  if (num_args2 > 0) \
407  expr2.computePartials(RADJOINT, partials+num_args1); \
408  } \
409  \
410  SACADO_INLINE_FUNCTION \
411  void getTangents(int i, value_type dots[]) const { \
412  expr1.getTangents(i, dots); \
413  expr2.getTangents(i, dots+num_args1); \
414  } \
415  \
416  template <int Arg> \
417  SACADO_INLINE_FUNCTION \
418  const value_type& getTangent(int i) const { \
419  if (Arg < num_args1) \
420  return expr1.template getTangent<Arg>(i); \
421  else \
422  return expr2.template getTangent<Arg-num_args1>(i); \
423  } \
424  \
425  SACADO_INLINE_FUNCTION \
426  bool isLinear() const { \
427  return LINEAR; \
428  } \
429  \
430  SACADO_INLINE_FUNCTION \
431  bool hasFastAccess() const { \
432  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
433  } \
434  \
435  SACADO_INLINE_FUNCTION \
436  const value_type dx(int i) const { \
437  return DX; \
438  } \
439  \
440  SACADO_INLINE_FUNCTION \
441  const value_type fastAccessDx(int i) const { \
442  return FASTACCESSDX; \
443  } \
444  \
445  SACADO_INLINE_FUNCTION \
446  const value_type* getDx(int j) const { \
447  if (j < num_args1) \
448  return expr1.getDx(j); \
449  else \
450  return expr2.getDx(j-num_args1); \
451  } \
452  \
453  SACADO_INLINE_FUNCTION \
454  int numActiveArgs() const { \
455  return expr1.numActiveArgs() + expr2.numActiveArgs(); \
456  } \
457  \
458  SACADO_INLINE_FUNCTION \
459  void computeActivePartials(const value_type& bar, \
460  value_type *partials) const { \
461  if (expr1.numActiveArgs() > 0) \
462  expr1.computePartials(LADJOINT, partials); \
463  if (expr2.numActiveArgs() > 0) \
464  expr2.computePartials(RADJOINT, partials+expr2.numActiveArgs()); \
465  } \
466  protected: \
467  \
468  typename ExprConstRef<ExprT1>::type expr1; \
469  typename ExprConstRef<ExprT2>::type expr2; \
470  \
471  }; \
472  \
473  template <typename ExprT1, typename T2> \
474  class Expr< OP<ExprT1, ConstExpr<T2> > > { \
475  \
476  public: \
477  \
478  typedef ConstExpr<T2> ExprT2; \
479  typedef typename ExprT1::value_type value_type_1; \
480  typedef typename ExprT2::value_type value_type_2; \
481  typedef typename Sacado::Promote<value_type_1, \
482  value_type_2>::type value_type; \
483  \
484  typedef typename ExprT1::scalar_type scalar_type_1; \
485  typedef typename ExprT2::scalar_type scalar_type_2; \
486  typedef typename Sacado::Promote<scalar_type_1, \
487  scalar_type_2>::type scalar_type; \
488  \
489  typedef typename ExprT1::base_expr_type base_expr_type_1; \
490  typedef typename ExprT2::base_expr_type base_expr_type_2; \
491  typedef typename Sacado::Promote<base_expr_type_1, \
492  base_expr_type_2>::type base_expr_type; \
493  \
494  static const int num_args = ExprT1::num_args; \
495  \
496  static const bool is_linear = CONST_LINEAR_2_2; \
497  \
498  SACADO_INLINE_FUNCTION \
499  Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
500  expr1(expr1_), expr2(expr2_) {} \
501  \
502  SACADO_INLINE_FUNCTION \
503  int size() const { \
504  return expr1.size(); \
505  } \
506  \
507  template <int Arg> \
508  SACADO_INLINE_FUNCTION \
509  bool isActive() const { \
510  return expr1.template isActive<Arg>(); \
511  } \
512  \
513  SACADO_INLINE_FUNCTION \
514  bool isActive2(int j) const { return expr1.isActive2(j); } \
515  \
516  SACADO_INLINE_FUNCTION \
517  bool updateValue() const { \
518  return expr1.updateValue(); \
519  } \
520  \
521  SACADO_INLINE_FUNCTION \
522  void cache() const {} \
523  \
524  SACADO_INLINE_FUNCTION \
525  value_type val() const { \
526  return VALUE; \
527  } \
528  \
529  SACADO_INLINE_FUNCTION \
530  void computePartials(const value_type& bar, \
531  value_type partials[]) const { \
532  expr1.computePartials(LADJOINT, partials); \
533  } \
534  \
535  SACADO_INLINE_FUNCTION \
536  void getTangents(int i, value_type dots[]) const { \
537  expr1.getTangents(i, dots); \
538  } \
539  \
540  template <int Arg> \
541  SACADO_INLINE_FUNCTION \
542  const value_type& getTangent(int i) const { \
543  return expr1.template getTangent<Arg>(i); \
544  } \
545  \
546  SACADO_INLINE_FUNCTION \
547  bool isLinear() const { \
548  return CONST_LINEAR_2; \
549  } \
550  \
551  SACADO_INLINE_FUNCTION \
552  bool hasFastAccess() const { \
553  return expr1.hasFastAccess(); \
554  } \
555  \
556  SACADO_INLINE_FUNCTION \
557  const value_type dx(int i) const { \
558  return CONST_DX_2; \
559  } \
560  \
561  SACADO_INLINE_FUNCTION \
562  const value_type fastAccessDx(int i) const { \
563  return CONST_FASTACCESSDX_2; \
564  } \
565  \
566  SACADO_INLINE_FUNCTION \
567  const value_type* getDx(int j) const { \
568  return expr1.getDx(j); \
569  } \
570  \
571  SACADO_INLINE_FUNCTION \
572  int numActiveArgs() const { \
573  return expr1.numActiveArgs(); \
574  } \
575  \
576  SACADO_INLINE_FUNCTION \
577  void computeActivePartials(const value_type& bar, \
578  value_type *partials) const { \
579  expr1.computePartials(LADJOINT, partials); \
580  } \
581  \
582  protected: \
583  \
584  typename ExprConstRef<ExprT1>::type expr1; \
585  typename ExprConstRef<ExprT2>::type expr2; \
586  \
587  }; \
588  \
589  template <typename T1, typename ExprT2> \
590  class Expr< OP<ConstExpr<T1>,ExprT2> > { \
591  \
592  public: \
593  \
594  typedef ConstExpr<T1> ExprT1; \
595  typedef typename ExprT1::value_type value_type_1; \
596  typedef typename ExprT2::value_type value_type_2; \
597  typedef typename Sacado::Promote<value_type_1, \
598  value_type_2>::type value_type; \
599  \
600  typedef typename ExprT1::scalar_type scalar_type_1; \
601  typedef typename ExprT2::scalar_type scalar_type_2; \
602  typedef typename Sacado::Promote<scalar_type_1, \
603  scalar_type_2>::type scalar_type; \
604  \
605  typedef typename ExprT1::base_expr_type base_expr_type_1; \
606  typedef typename ExprT2::base_expr_type base_expr_type_2; \
607  typedef typename Sacado::Promote<base_expr_type_1, \
608  base_expr_type_2>::type base_expr_type; \
609  \
610  static const int num_args = ExprT2::num_args; \
611  \
612  static const bool is_linear = CONST_LINEAR_1_2; \
613  \
614  SACADO_INLINE_FUNCTION \
615  Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
616  expr1(expr1_), expr2(expr2_) {} \
617  \
618  SACADO_INLINE_FUNCTION \
619  int size() const { \
620  return expr2.size(); \
621  } \
622  \
623  template <int Arg> \
624  SACADO_INLINE_FUNCTION \
625  bool isActive() const { \
626  return expr2.template isActive<Arg>(); \
627  } \
628  \
629  SACADO_INLINE_FUNCTION \
630  bool isActive2(int j) const { return expr2.isActive2(j); } \
631  \
632  SACADO_INLINE_FUNCTION \
633  bool updateValue() const { \
634  return expr2.updateValue(); \
635  } \
636  \
637  SACADO_INLINE_FUNCTION \
638  void cache() const {} \
639  \
640  SACADO_INLINE_FUNCTION \
641  value_type val() const { \
642  return VALUE; \
643  } \
644  \
645  SACADO_INLINE_FUNCTION \
646  void computePartials(const value_type& bar, \
647  value_type partials[]) const { \
648  expr2.computePartials(RADJOINT, partials); \
649  } \
650  \
651  SACADO_INLINE_FUNCTION \
652  void getTangents(int i, value_type dots[]) const { \
653  expr2.getTangents(i, dots); \
654  } \
655  \
656  template <int Arg> \
657  SACADO_INLINE_FUNCTION \
658  const value_type& getTangent(int i) const { \
659  return expr2.template getTangent<Arg>(i); \
660  } \
661  \
662  SACADO_INLINE_FUNCTION \
663  bool isLinear() const { \
664  return CONST_LINEAR_1; \
665  } \
666  \
667  SACADO_INLINE_FUNCTION \
668  bool hasFastAccess() const { \
669  return expr2.hasFastAccess(); \
670  } \
671  \
672  SACADO_INLINE_FUNCTION \
673  const value_type dx(int i) const { \
674  return CONST_DX_1; \
675  } \
676  \
677  SACADO_INLINE_FUNCTION \
678  const value_type fastAccessDx(int i) const { \
679  return CONST_FASTACCESSDX_1; \
680  } \
681  \
682  SACADO_INLINE_FUNCTION \
683  const value_type* getDx(int j) const { \
684  return expr2.getDx(j); \
685  } \
686  \
687  SACADO_INLINE_FUNCTION \
688  int numActiveArgs() const { \
689  return expr2.numActiveArgs(); \
690  } \
691  \
692  SACADO_INLINE_FUNCTION \
693  void computeActivePartials(const value_type& bar, \
694  value_type *partials) const { \
695  expr2.computePartials(RADJOINT, partials); \
696  } \
697  protected: \
698  \
699  typename ExprConstRef<ExprT1>::type expr1; \
700  typename ExprConstRef<ExprT2>::type expr2; \
701  \
702  }; \
703  \
704  template <typename T1, typename T2> \
705  SACADO_INLINE_FUNCTION \
706  SACADO_FAD_OP_ENABLE_EXPR_EXPR(OP) \
707  OPNAME (const T1& expr1, const T2& expr2) \
708  { \
709  typedef OP< T1, T2 > expr_t; \
710  \
711  return Expr<expr_t>(expr1, expr2); \
712  } \
713  \
714  template <typename T> \
715  SACADO_INLINE_FUNCTION \
716  Expr< OP< Expr<T>, Expr<T> > > \
717  OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
718  { \
719  typedef OP< Expr<T>, Expr<T> > expr_t; \
720  \
721  return Expr<expr_t>(expr1, expr2); \
722  } \
723  \
724  template <typename T> \
725  SACADO_INLINE_FUNCTION \
726  Expr< OP< ConstExpr<typename Expr<T>::value_type>, \
727  Expr<T> > > \
728  OPNAME (const typename Expr<T>::value_type& c, \
729  const Expr<T>& expr) \
730  { \
731  typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
732  typedef OP< ConstT, Expr<T> > expr_t; \
733  \
734  return Expr<expr_t>(ConstT(c), expr); \
735  } \
736  \
737  template <typename T> \
738  SACADO_INLINE_FUNCTION \
739  Expr< OP< Expr<T>, \
740  ConstExpr<typename Expr<T>::value_type> > > \
741  OPNAME (const Expr<T>& expr, \
742  const typename Expr<T>::value_type& c) \
743  { \
744  typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
745  typedef OP< Expr<T>, ConstT > expr_t; \
746  \
747  return Expr<expr_t>(expr, ConstT(c)); \
748  } \
749  \
750  template <typename T> \
751  SACADO_INLINE_FUNCTION \
752  SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP) \
753  OPNAME (const typename Expr<T>::scalar_type& c, \
754  const Expr<T>& expr) \
755  { \
756  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
757  typedef OP< ConstT, Expr<T> > expr_t; \
758  \
759  return Expr<expr_t>(ConstT(c), expr); \
760  } \
761  \
762  template <typename T> \
763  SACADO_INLINE_FUNCTION \
764  SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP) \
765  OPNAME (const Expr<T>& expr, \
766  const typename Expr<T>::scalar_type& c) \
767  { \
768  typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
769  typedef OP< Expr<T>, ConstT > expr_t; \
770  \
771  return Expr<expr_t>(expr, ConstT(c)); \
772  } \
773  } \
774 }
775 
776 
777 FAD_BINARYOP_MACRO(operator+,
778  AdditionOp,
779  expr1.val() + expr2.val(),
780  bar,
781  bar,
782  expr1.isLinear() && expr2.isLinear(),
783  expr2.isLinear(),
784  expr1.isLinear(),
785  ExprT1::is_linear && ExprT2::is_linear,
786  ExprT2::is_linear,
787  ExprT1::is_linear,
788  expr1.dx(i) + expr2.dx(i),
789  expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
790  expr2.dx(i),
791  expr1.dx(i),
792  expr2.fastAccessDx(i),
793  expr1.fastAccessDx(i))
794 FAD_BINARYOP_MACRO(operator-,
796  expr1.val() - expr2.val(),
797  bar,
798  -bar,
799  expr1.isLinear() && expr2.isLinear(),
800  expr2.isLinear(),
801  expr1.isLinear(),
802  ExprT1::is_linear && ExprT2::is_linear,
803  ExprT2::is_linear,
804  ExprT1::is_linear,
805  expr1.dx(i) - expr2.dx(i),
806  expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
807  -expr2.dx(i),
808  expr1.dx(i),
809  -expr2.fastAccessDx(i),
810  expr1.fastAccessDx(i))
811 FAD_BINARYOP_MACRO(operator*,
813  expr1.val() * expr2.val(),
814  bar*expr2.val(),
815  bar*expr1.val(),
816  false,
817  expr2.isLinear(),
818  expr1.isLinear(),
819  false,
820  ExprT2::is_linear,
821  ExprT1::is_linear,
822  expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
823  expr1.val()*expr2.fastAccessDx(i) +
824  expr1.fastAccessDx(i)*expr2.val(),
825  expr1.val()*expr2.dx(i),
826  expr1.dx(i)*expr2.val(),
827  expr1.val()*expr2.fastAccessDx(i),
828  expr1.fastAccessDx(i)*expr2.val())
829 FAD_BINARYOP_MACRO(operator/,
830  DivisionOp,
831  expr1.val() / expr2.val(),
832  bar/expr2.val(),
833  -bar*expr1.val()/(expr2.val()*expr2.val()),
834  false,
835  false,
836  expr1.isLinear(),
837  false,
838  false,
839  ExprT1::is_linear,
840  (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
841  (expr2.val()*expr2.val()),
842  (expr1.fastAccessDx(i)*expr2.val() -
843  expr2.fastAccessDx(i)*expr1.val()) /
844  (expr2.val()*expr2.val()),
845  -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
846  expr1.dx(i)/expr2.val(),
847  -expr2.fastAccessDx(i)*expr1.val() / (expr2.val()*expr2.val()),
848  expr1.fastAccessDx(i)/expr2.val())
850  Atan2Op,
851  std::atan2(expr1.val(), expr2.val()),
852  bar*expr2.val()/
853  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
854  -bar*expr1.val()/
855  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
856  false,
857  false,
858  false,
859  false,
860  false,
861  false,
862  (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
863  (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
864  (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
865  (-expr1.val()*expr2.dx(i)) / (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
866  (expr2.val()*expr1.dx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
867  (-expr1.val()*expr2.fastAccessDx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
868  (expr2.val()*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()))
870  PowerOp,
871  std::pow(expr1.val(), expr2.val()),
872  expr2.size() == 0 && expr2.val() == value_type(1) ? bar : expr1.val() == value_type(0) ? value_type(0) : value_type(bar*std::pow(expr1.val(),expr2.val())*expr2.val()/expr1.val()),
873  expr1.val() == value_type(0) ? value_type(0) : value_type(bar*std::pow(expr1.val(),expr2.val())*std::log(expr1.val())),
874  false,
875  false,
876  false,
877  false,
878  false,
879  false,
880  expr1.val() == value_type(0) ? value_type(0) : value_type((expr2.dx(i)*std::log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*std::pow(expr1.val(),expr2.val())),
881  expr1.val() == value_type(0) ? value_type(0) : value_type((expr2.fastAccessDx(i)*std::log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*std::pow(expr1.val(),expr2.val())),
882  expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.dx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
883  expr2.val() == value_type(1) ? expr1.dx(i) : expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.val()*expr1.dx(i)/expr1.val()*std::pow(expr1.val(),expr2.val())),
884  expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.fastAccessDx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
885  expr2.val() == value_type(1) ? expr1.fastAccessDx(i) : expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.val()*expr1.fastAccessDx(i)/expr1.val()*std::pow(expr1.val(),expr2.val())))
887  MaxOp,
888  expr1.val() >= expr2.val() ? expr1.val() : expr2.val(),
889  expr1.val() >= expr2.val() ? bar : value_type(0.),
890  expr2.val() > expr1.val() ? bar : value_type(0.),
891  expr1.isLinear() && expr2.isLinear(),
892  expr2.isLinear(),
893  expr1.isLinear(),
894  ExprT1::is_linear && ExprT2::is_linear,
895  ExprT2::is_linear,
896  ExprT1::is_linear,
897  expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
898  expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) :
899  expr2.fastAccessDx(i),
900  expr1.val() >= expr2.val() ? value_type(0) : expr2.dx(i),
901  expr1.val() >= expr2.val() ? expr1.dx(i) : value_type(0),
902  expr1.val() >= expr2.val() ? value_type(0) :
903  expr2.fastAccessDx(i),
904  expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) :
905  value_type(0))
907  MinOp,
908  expr1.val() <= expr2.val() ? expr1.val() : expr2.val(),
909  expr1.val() <= expr2.val() ? bar : value_type(0.),
910  expr2.val() < expr1.val() ? bar : value_type(0.),
911  expr1.isLinear() && expr2.isLinear(),
912  expr2.isLinear(),
913  expr1.isLinear(),
914  ExprT1::is_linear && ExprT2::is_linear,
915  ExprT2::is_linear,
916  ExprT1::is_linear,
917  expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
918  expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) :
919  expr2.fastAccessDx(i),
920  expr1.val() <= expr2.val() ? value_type(0) : expr2.dx(i),
921  expr1.val() <= expr2.val() ? expr1.dx(i) : value_type(0),
922  expr1.val() <= expr2.val() ? value_type(0) :
923  expr2.fastAccessDx(i),
924  expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) :
925  value_type(0))
926 
927 #undef FAD_BINARYOP_MACRO
928 
929 //-------------------------- Relational Operators -----------------------
930 
931 #define FAD_RELOP_MACRO(OP) \
932 namespace Sacado { \
933  namespace ELRFad { \
934  template <typename ExprT1, typename ExprT2> \
935  SACADO_INLINE_FUNCTION \
936  bool \
937  operator OP (const Expr<ExprT1>& expr1, \
938  const Expr<ExprT2>& expr2) \
939  { \
940  return expr1.val() OP expr2.val(); \
941  } \
942  \
943  template <typename ExprT2> \
944  SACADO_INLINE_FUNCTION \
945  bool \
946  operator OP (const typename Expr<ExprT2>::value_type& a, \
947  const Expr<ExprT2>& expr2) \
948  { \
949  return a OP expr2.val(); \
950  } \
951  \
952  template <typename ExprT1> \
953  SACADO_INLINE_FUNCTION \
954  bool \
955  operator OP (const Expr<ExprT1>& expr1, \
956  const typename Expr<ExprT1>::value_type& b) \
957  { \
958  return expr1.val() OP b; \
959  } \
960  } \
961 }
962 
963 FAD_RELOP_MACRO(==)
964 FAD_RELOP_MACRO(!=)
967 FAD_RELOP_MACRO(<=)
968 FAD_RELOP_MACRO(>=)
969 FAD_RELOP_MACRO(<<=)
970 FAD_RELOP_MACRO(>>=)
973 
974 #undef FAD_RELOP_MACRO
975 
976 namespace Sacado {
977 
978  namespace ELRFad {
979 
980  template <typename ExprT>
982  bool operator ! (const Expr<ExprT>& expr)
983  {
984  return ! expr.val();
985  }
986 
987  } // namespace ELRFad
988 
989 } // namespace Sacado
990 
991 //-------------------------- Boolean Operators -----------------------
992 namespace Sacado {
993 
994  namespace ELRFad {
995 
996  template <typename ExprT>
998  bool toBool(const Expr<ExprT>& x) {
999  bool is_zero = (x.val() == 0.0);
1000  for (int i=0; i<x.size(); i++)
1001  is_zero = is_zero && (x.dx(i) == 0.0);
1002  return !is_zero;
1003  }
1004 
1005  } // namespace Fad
1006 
1007 } // namespace Sacado
1008 
1009 #define FAD_BOOL_MACRO(OP) \
1010 namespace Sacado { \
1011  namespace ELRFad { \
1012  template <typename ExprT1, typename ExprT2> \
1013  SACADO_INLINE_FUNCTION \
1014  bool \
1015  operator OP (const Expr<ExprT1>& expr1, \
1016  const Expr<ExprT2>& expr2) \
1017  { \
1018  return toBool(expr1) OP toBool(expr2); \
1019  } \
1020  \
1021  template <typename ExprT2> \
1022  SACADO_INLINE_FUNCTION \
1023  bool \
1024  operator OP (const typename Expr<ExprT2>::value_type& a, \
1025  const Expr<ExprT2>& expr2) \
1026  { \
1027  return a OP toBool(expr2); \
1028  } \
1029  \
1030  template <typename ExprT1> \
1031  SACADO_INLINE_FUNCTION \
1032  bool \
1033  operator OP (const Expr<ExprT1>& expr1, \
1034  const typename Expr<ExprT1>::value_type& b) \
1035  { \
1036  return toBool(expr1) OP b; \
1037  } \
1038  } \
1039 }
1040 
1041 FAD_BOOL_MACRO(&&)
1042 FAD_BOOL_MACRO(||)
1043 
1044 #undef FAD_BOOL_MACRO
1045 
1046 //-------------------------- I/O Operators -----------------------
1047 
1048 namespace Sacado {
1049 
1050  namespace ELRFad {
1051 
1052  template <typename ExprT>
1053  std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
1054  os << x.val() << " [";
1055 
1056  for (int i=0; i< x.size(); i++) {
1057  os << " " << x.dx(i);
1058  }
1059 
1060  os << " ]";
1061  return os;
1062  }
1063 
1064  } // namespace Fad
1065 
1066 } // namespace Sacado
1067 
1068 #endif // SACADO_FAD_OPS_HPP
cbrt(expr.val())
expr expr SinOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
asinh(expr.val())
SACADO_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
asin(expr.val())
cosh(expr.val())
expr expr dx(i)
abs(expr.val())
#define FAD_BOOL_MACRO(OP)
atanh(expr.val())
expr expr CoshOp
expr expr ATanhOp
SACADO_INLINE_FUNCTION bool toBool(const Expr< ExprT > &x)
expr expr TanhOp
expr expr SqrtOp
expr expr ASinhOp
expr bar
atan(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
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()
tanh(expr.val())
expr expr CosOp
#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)
expr expr ATanOp
SimpleFad< ValueT > min(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ACosOp
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
atan2(expr1.val(), expr2.val())
#define FAD_RELOP_MACRO(OP)
sin(expr.val())
Wrapper for a generic expression template.
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
SACADO_INLINE_FUNCTION T safe_sqrt(const T &x)
log(expr.val())
expr expr ACoshOp
expr expr Log10Op
expr expr SinhOp
acosh(expr.val())
acos(expr.val())
SimpleFad< ValueT > max(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ASinOp
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)
expr expr expr bar false
#define SACADO_INLINE_FUNCTION
exp(expr.val())
expr expr expr ExpOp
fabs(expr.val())
expr expr AbsOp
expr expr TanOp
log10(expr.val())
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