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