Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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) \
61namespace 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
173 UnaryPlusOp,
174 expr.val(),
175 bar,
176 true,
177 expr.dx(i),
178 expr.fastAccessDx(i))
179FAD_UNARYOP_MACRO(operator-,
181 -expr.val(),
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())))
214FAD_UNARYOP_MACRO(safe_sqrt,
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) \
348namespace 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
796FAD_BINARYOP_MACRO(operator+,
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))
813FAD_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))
830FAD_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())
848FAD_BINARYOP_MACRO(operator/,
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) \
951namespace 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
992
993#undef FAD_RELOP_MACRO
994
995namespace 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 -----------------------
1011namespace 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) \
1029namespace 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
1062
1063#undef FAD_BOOL_MACRO
1064
1065//-------------------------- I/O Operators -----------------------
1066
1067namespace 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
#define SACADO_INLINE_FUNCTION
#define FAD_RELOP_MACRO(OP)
#define FAD_UNARYOP_MACRO(OPNAME, OP, VALUE, ADJOINT, LINEAR, DX, FASTACCESSDX)
expr expr expr bar false
#define FAD_BOOL_MACRO(OP)
#define FAD_BINARYOP_MACRO( OPNAME, OP, VALUE, LADJOINT, RADJOINT, LINEAR, CONST_LINEAR_1, CONST_LINEAR_2, LINEAR_2, CONST_LINEAR_1_2, CONST_LINEAR_2_2, DX, FASTACCESSDX, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
fabs(expr.val())
asinh(expr.val())
log(expr.val())
tan(expr.val())
abs(expr.val())
expr expr SinOp
expr expr SinhOp
cos(expr.val())
expr expr SqrtOp
expr expr dx(i)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
expr expr ACosOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
expr expr ATanOp
cosh(expr.val())
expr expr ACoshOp
acos(expr.val())
sin(expr.val())
sinh(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
expr expr ASinOp
cbrt(expr.val())
log10(expr.val())
expr expr TanhOp
expr expr TanOp
expr expr CoshOp
#define FAD_BINARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, CDX1, CDX2, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
#define FAD_RELOP_MACRO(OP)
expr expr ASinhOp
exp(expr.val())
atan(expr.val())
expr expr Log10Op
#define FAD_BOOL_MACRO(OP)
acosh(expr.val())
expr expr expr ExpOp
expr expr AbsOp
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
atan2(expr1.val(), expr2.val())
sqrt(expr.val())
expr expr ATanhOp
tanh(expr.val())
atanh(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 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 Atan2Op
expr expr CosOp
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
asin(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
adouble max(const adouble &a, const adouble &b)
adouble min(const adouble &a, const adouble &b)
Wrapper for a generic expression template.
Namespace for expression-level reverse forward-mode AD classes.
std::ostream & operator<<(std::ostream &os, const GeneralFad< T, Storage > &x)
SACADO_INLINE_FUNCTION bool toBool(const Expr< T > &xx)