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)))
332#ifdef HAVE_SACADO_CXX11
334 CbrtOp,
335 std::cbrt(expr.val()),
336 bar/(value_type(3)*std::cbrt(expr.val()*expr.val())),
337 false,
338 expr.dx(i)/(value_type(3)*std::cbrt(expr.val()*expr.val())),
339 expr.fastAccessDx(i)/(value_type(3)*std::cbrt(expr.val()*expr.val())))
340#endif
341
342#undef FAD_UNARYOP_MACRO
343
344#define FAD_BINARYOP_MACRO( \
345 OPNAME,OP,VALUE,LADJOINT,RADJOINT, \
346 LINEAR,CONST_LINEAR_1, CONST_LINEAR_2, \
347 LINEAR_2,CONST_LINEAR_1_2, CONST_LINEAR_2_2, \
348 DX,FASTACCESSDX,CONST_DX_1,CONST_DX_2, \
349 CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
350namespace Sacado { \
351 namespace ELRFad { \
352 \
353 template <typename ExprT1, typename ExprT2> \
354 class OP {}; \
355 \
356 template <typename ExprT1, typename ExprT2> \
357 class Expr< OP<ExprT1,ExprT2> > { \
358 \
359 public: \
360 \
361 typedef typename ExprT1::value_type value_type_1; \
362 typedef typename ExprT2::value_type value_type_2; \
363 typedef typename Sacado::Promote<value_type_1, \
364 value_type_2>::type value_type; \
365 \
366 typedef typename ExprT1::scalar_type scalar_type_1; \
367 typedef typename ExprT2::scalar_type scalar_type_2; \
368 typedef typename Sacado::Promote<scalar_type_1, \
369 scalar_type_2>::type scalar_type; \
370 \
371 typedef typename ExprT1::base_expr_type base_expr_type_1; \
372 typedef typename ExprT2::base_expr_type base_expr_type_2; \
373 typedef typename Sacado::Promote<base_expr_type_1, \
374 base_expr_type_2>::type base_expr_type; \
375 \
376 static const int num_args1 = ExprT1::num_args; \
377 static const int num_args2 = ExprT2::num_args; \
378 static const int num_args = num_args1 + num_args2; \
379 \
380 static const bool is_linear = LINEAR_2; \
381 \
382 SACADO_INLINE_FUNCTION \
383 Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
384 expr1(expr1_), expr2(expr2_) {} \
385 \
386 SACADO_INLINE_FUNCTION \
387 int size() const { \
388 int sz1 = expr1.size(), sz2 = expr2.size(); \
389 return sz1 > sz2 ? sz1 : sz2; \
390 } \
391 \
392 template <int Arg> \
393 SACADO_INLINE_FUNCTION \
394 bool isActive() const { \
395 if (Arg < num_args1) \
396 return expr1.template isActive<Arg>(); \
397 else \
398 return expr2.template isActive<Arg-num_args1>(); \
399 } \
400 \
401 SACADO_INLINE_FUNCTION \
402 bool isActive2(int j) const { \
403 if (j < num_args1) \
404 return expr1.isActive2(j); \
405 else \
406 return expr2.isActive2(j); \
407 } \
408 \
409 SACADO_INLINE_FUNCTION \
410 bool updateValue() const { \
411 return expr1.updateValue() && expr2.updateValue(); \
412 } \
413 \
414 SACADO_INLINE_FUNCTION \
415 void cache() const {} \
416 \
417 SACADO_INLINE_FUNCTION \
418 value_type val() const { \
419 return VALUE; \
420 } \
421 \
422 SACADO_INLINE_FUNCTION \
423 void computePartials(const value_type& bar, \
424 value_type partials[]) const { \
425 if (num_args1 > 0) \
426 expr1.computePartials(LADJOINT, partials); \
427 if (num_args2 > 0) \
428 expr2.computePartials(RADJOINT, partials+num_args1); \
429 } \
430 \
431 SACADO_INLINE_FUNCTION \
432 void getTangents(int i, value_type dots[]) const { \
433 expr1.getTangents(i, dots); \
434 expr2.getTangents(i, dots+num_args1); \
435 } \
436 \
437 template <int Arg> \
438 SACADO_INLINE_FUNCTION \
439 const value_type& getTangent(int i) const { \
440 if (Arg < num_args1) \
441 return expr1.template getTangent<Arg>(i); \
442 else \
443 return expr2.template getTangent<Arg-num_args1>(i); \
444 } \
445 \
446 SACADO_INLINE_FUNCTION \
447 bool isLinear() const { \
448 return LINEAR; \
449 } \
450 \
451 SACADO_INLINE_FUNCTION \
452 bool hasFastAccess() const { \
453 return expr1.hasFastAccess() && expr2.hasFastAccess(); \
454 } \
455 \
456 SACADO_INLINE_FUNCTION \
457 const value_type dx(int i) const { \
458 return DX; \
459 } \
460 \
461 SACADO_INLINE_FUNCTION \
462 const value_type fastAccessDx(int i) const { \
463 return FASTACCESSDX; \
464 } \
465 \
466 SACADO_INLINE_FUNCTION \
467 const value_type* getDx(int j) const { \
468 if (j < num_args1) \
469 return expr1.getDx(j); \
470 else \
471 return expr2.getDx(j-num_args1); \
472 } \
473 \
474 SACADO_INLINE_FUNCTION \
475 int numActiveArgs() const { \
476 return expr1.numActiveArgs() + expr2.numActiveArgs(); \
477 } \
478 \
479 SACADO_INLINE_FUNCTION \
480 void computeActivePartials(const value_type& bar, \
481 value_type *partials) const { \
482 if (expr1.numActiveArgs() > 0) \
483 expr1.computePartials(LADJOINT, partials); \
484 if (expr2.numActiveArgs() > 0) \
485 expr2.computePartials(RADJOINT, partials+expr2.numActiveArgs()); \
486 } \
487 protected: \
488 \
489 typename ExprConstRef<ExprT1>::type expr1; \
490 typename ExprConstRef<ExprT2>::type expr2; \
491 \
492 }; \
493 \
494 template <typename ExprT1, typename T2> \
495 class Expr< OP<ExprT1, ConstExpr<T2> > > { \
496 \
497 public: \
498 \
499 typedef ConstExpr<T2> ExprT2; \
500 typedef typename ExprT1::value_type value_type_1; \
501 typedef typename ExprT2::value_type value_type_2; \
502 typedef typename Sacado::Promote<value_type_1, \
503 value_type_2>::type value_type; \
504 \
505 typedef typename ExprT1::scalar_type scalar_type_1; \
506 typedef typename ExprT2::scalar_type scalar_type_2; \
507 typedef typename Sacado::Promote<scalar_type_1, \
508 scalar_type_2>::type scalar_type; \
509 \
510 typedef typename ExprT1::base_expr_type base_expr_type_1; \
511 typedef typename ExprT2::base_expr_type base_expr_type_2; \
512 typedef typename Sacado::Promote<base_expr_type_1, \
513 base_expr_type_2>::type base_expr_type; \
514 \
515 static const int num_args = ExprT1::num_args; \
516 \
517 static const bool is_linear = CONST_LINEAR_2_2; \
518 \
519 SACADO_INLINE_FUNCTION \
520 Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
521 expr1(expr1_), expr2(expr2_) {} \
522 \
523 SACADO_INLINE_FUNCTION \
524 int size() const { \
525 return expr1.size(); \
526 } \
527 \
528 template <int Arg> \
529 SACADO_INLINE_FUNCTION \
530 bool isActive() const { \
531 return expr1.template isActive<Arg>(); \
532 } \
533 \
534 SACADO_INLINE_FUNCTION \
535 bool isActive2(int j) const { return expr1.isActive2(j); } \
536 \
537 SACADO_INLINE_FUNCTION \
538 bool updateValue() const { \
539 return expr1.updateValue(); \
540 } \
541 \
542 SACADO_INLINE_FUNCTION \
543 void cache() const {} \
544 \
545 SACADO_INLINE_FUNCTION \
546 value_type val() const { \
547 return VALUE; \
548 } \
549 \
550 SACADO_INLINE_FUNCTION \
551 void computePartials(const value_type& bar, \
552 value_type partials[]) const { \
553 expr1.computePartials(LADJOINT, partials); \
554 } \
555 \
556 SACADO_INLINE_FUNCTION \
557 void getTangents(int i, value_type dots[]) const { \
558 expr1.getTangents(i, dots); \
559 } \
560 \
561 template <int Arg> \
562 SACADO_INLINE_FUNCTION \
563 const value_type& getTangent(int i) const { \
564 return expr1.template getTangent<Arg>(i); \
565 } \
566 \
567 SACADO_INLINE_FUNCTION \
568 bool isLinear() const { \
569 return CONST_LINEAR_2; \
570 } \
571 \
572 SACADO_INLINE_FUNCTION \
573 bool hasFastAccess() const { \
574 return expr1.hasFastAccess(); \
575 } \
576 \
577 SACADO_INLINE_FUNCTION \
578 const value_type dx(int i) const { \
579 return CONST_DX_2; \
580 } \
581 \
582 SACADO_INLINE_FUNCTION \
583 const value_type fastAccessDx(int i) const { \
584 return CONST_FASTACCESSDX_2; \
585 } \
586 \
587 SACADO_INLINE_FUNCTION \
588 const value_type* getDx(int j) const { \
589 return expr1.getDx(j); \
590 } \
591 \
592 SACADO_INLINE_FUNCTION \
593 int numActiveArgs() const { \
594 return expr1.numActiveArgs(); \
595 } \
596 \
597 SACADO_INLINE_FUNCTION \
598 void computeActivePartials(const value_type& bar, \
599 value_type *partials) const { \
600 expr1.computePartials(LADJOINT, partials); \
601 } \
602 \
603 protected: \
604 \
605 typename ExprConstRef<ExprT1>::type expr1; \
606 typename ExprConstRef<ExprT2>::type expr2; \
607 \
608 }; \
609 \
610 template <typename T1, typename ExprT2> \
611 class Expr< OP<ConstExpr<T1>,ExprT2> > { \
612 \
613 public: \
614 \
615 typedef ConstExpr<T1> ExprT1; \
616 typedef typename ExprT1::value_type value_type_1; \
617 typedef typename ExprT2::value_type value_type_2; \
618 typedef typename Sacado::Promote<value_type_1, \
619 value_type_2>::type value_type; \
620 \
621 typedef typename ExprT1::scalar_type scalar_type_1; \
622 typedef typename ExprT2::scalar_type scalar_type_2; \
623 typedef typename Sacado::Promote<scalar_type_1, \
624 scalar_type_2>::type scalar_type; \
625 \
626 typedef typename ExprT1::base_expr_type base_expr_type_1; \
627 typedef typename ExprT2::base_expr_type base_expr_type_2; \
628 typedef typename Sacado::Promote<base_expr_type_1, \
629 base_expr_type_2>::type base_expr_type; \
630 \
631 static const int num_args = ExprT2::num_args; \
632 \
633 static const bool is_linear = CONST_LINEAR_1_2; \
634 \
635 SACADO_INLINE_FUNCTION \
636 Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
637 expr1(expr1_), expr2(expr2_) {} \
638 \
639 SACADO_INLINE_FUNCTION \
640 int size() const { \
641 return expr2.size(); \
642 } \
643 \
644 template <int Arg> \
645 SACADO_INLINE_FUNCTION \
646 bool isActive() const { \
647 return expr2.template isActive<Arg>(); \
648 } \
649 \
650 SACADO_INLINE_FUNCTION \
651 bool isActive2(int j) const { return expr2.isActive2(j); } \
652 \
653 SACADO_INLINE_FUNCTION \
654 bool updateValue() const { \
655 return expr2.updateValue(); \
656 } \
657 \
658 SACADO_INLINE_FUNCTION \
659 void cache() const {} \
660 \
661 SACADO_INLINE_FUNCTION \
662 value_type val() const { \
663 return VALUE; \
664 } \
665 \
666 SACADO_INLINE_FUNCTION \
667 void computePartials(const value_type& bar, \
668 value_type partials[]) const { \
669 expr2.computePartials(RADJOINT, partials); \
670 } \
671 \
672 SACADO_INLINE_FUNCTION \
673 void getTangents(int i, value_type dots[]) const { \
674 expr2.getTangents(i, dots); \
675 } \
676 \
677 template <int Arg> \
678 SACADO_INLINE_FUNCTION \
679 const value_type& getTangent(int i) const { \
680 return expr2.template getTangent<Arg>(i); \
681 } \
682 \
683 SACADO_INLINE_FUNCTION \
684 bool isLinear() const { \
685 return CONST_LINEAR_1; \
686 } \
687 \
688 SACADO_INLINE_FUNCTION \
689 bool hasFastAccess() const { \
690 return expr2.hasFastAccess(); \
691 } \
692 \
693 SACADO_INLINE_FUNCTION \
694 const value_type dx(int i) const { \
695 return CONST_DX_1; \
696 } \
697 \
698 SACADO_INLINE_FUNCTION \
699 const value_type fastAccessDx(int i) const { \
700 return CONST_FASTACCESSDX_1; \
701 } \
702 \
703 SACADO_INLINE_FUNCTION \
704 const value_type* getDx(int j) const { \
705 return expr2.getDx(j); \
706 } \
707 \
708 SACADO_INLINE_FUNCTION \
709 int numActiveArgs() const { \
710 return expr2.numActiveArgs(); \
711 } \
712 \
713 SACADO_INLINE_FUNCTION \
714 void computeActivePartials(const value_type& bar, \
715 value_type *partials) const { \
716 expr2.computePartials(RADJOINT, partials); \
717 } \
718 protected: \
719 \
720 typename ExprConstRef<ExprT1>::type expr1; \
721 typename ExprConstRef<ExprT2>::type expr2; \
722 \
723 }; \
724 \
725 template <typename T1, typename T2> \
726 SACADO_INLINE_FUNCTION \
727 SACADO_FAD_OP_ENABLE_EXPR_EXPR(OP) \
728 OPNAME (const T1& expr1, const T2& expr2) \
729 { \
730 typedef OP< T1, T2 > expr_t; \
731 \
732 return Expr<expr_t>(expr1, expr2); \
733 } \
734 \
735 template <typename T> \
736 SACADO_INLINE_FUNCTION \
737 Expr< OP< Expr<T>, Expr<T> > > \
738 OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
739 { \
740 typedef OP< Expr<T>, Expr<T> > expr_t; \
741 \
742 return Expr<expr_t>(expr1, expr2); \
743 } \
744 \
745 template <typename T> \
746 SACADO_INLINE_FUNCTION \
747 Expr< OP< ConstExpr<typename Expr<T>::value_type>, \
748 Expr<T> > > \
749 OPNAME (const typename Expr<T>::value_type& c, \
750 const Expr<T>& expr) \
751 { \
752 typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
753 typedef OP< ConstT, Expr<T> > expr_t; \
754 \
755 return Expr<expr_t>(ConstT(c), expr); \
756 } \
757 \
758 template <typename T> \
759 SACADO_INLINE_FUNCTION \
760 Expr< OP< Expr<T>, \
761 ConstExpr<typename Expr<T>::value_type> > > \
762 OPNAME (const Expr<T>& expr, \
763 const typename Expr<T>::value_type& c) \
764 { \
765 typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
766 typedef OP< Expr<T>, ConstT > expr_t; \
767 \
768 return Expr<expr_t>(expr, ConstT(c)); \
769 } \
770 \
771 template <typename T> \
772 SACADO_INLINE_FUNCTION \
773 SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP) \
774 OPNAME (const typename Expr<T>::scalar_type& c, \
775 const Expr<T>& expr) \
776 { \
777 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
778 typedef OP< ConstT, Expr<T> > expr_t; \
779 \
780 return Expr<expr_t>(ConstT(c), expr); \
781 } \
782 \
783 template <typename T> \
784 SACADO_INLINE_FUNCTION \
785 SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP) \
786 OPNAME (const Expr<T>& expr, \
787 const typename Expr<T>::scalar_type& c) \
788 { \
789 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
790 typedef OP< Expr<T>, ConstT > expr_t; \
791 \
792 return Expr<expr_t>(expr, ConstT(c)); \
793 } \
794 } \
795}
796
797
798FAD_BINARYOP_MACRO(operator+,
800 expr1.val() + expr2.val(),
801 bar,
802 bar,
803 expr1.isLinear() && expr2.isLinear(),
804 expr2.isLinear(),
805 expr1.isLinear(),
806 ExprT1::is_linear && ExprT2::is_linear,
807 ExprT2::is_linear,
808 ExprT1::is_linear,
809 expr1.dx(i) + expr2.dx(i),
810 expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
811 expr2.dx(i),
812 expr1.dx(i),
813 expr2.fastAccessDx(i),
814 expr1.fastAccessDx(i))
815FAD_BINARYOP_MACRO(operator-,
817 expr1.val() - expr2.val(),
818 bar,
819 -bar,
820 expr1.isLinear() && expr2.isLinear(),
821 expr2.isLinear(),
822 expr1.isLinear(),
823 ExprT1::is_linear && ExprT2::is_linear,
824 ExprT2::is_linear,
825 ExprT1::is_linear,
826 expr1.dx(i) - expr2.dx(i),
827 expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
828 -expr2.dx(i),
829 expr1.dx(i),
830 -expr2.fastAccessDx(i),
831 expr1.fastAccessDx(i))
832FAD_BINARYOP_MACRO(operator*,
834 expr1.val() * expr2.val(),
835 bar*expr2.val(),
836 bar*expr1.val(),
837 false,
838 expr2.isLinear(),
839 expr1.isLinear(),
840 false,
841 ExprT2::is_linear,
842 ExprT1::is_linear,
843 expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
844 expr1.val()*expr2.fastAccessDx(i) +
845 expr1.fastAccessDx(i)*expr2.val(),
846 expr1.val()*expr2.dx(i),
847 expr1.dx(i)*expr2.val(),
848 expr1.val()*expr2.fastAccessDx(i),
849 expr1.fastAccessDx(i)*expr2.val())
850FAD_BINARYOP_MACRO(operator/,
852 expr1.val() / expr2.val(),
853 bar/expr2.val(),
854 -bar*expr1.val()/(expr2.val()*expr2.val()),
855 false,
856 false,
857 expr1.isLinear(),
858 false,
859 false,
860 ExprT1::is_linear,
861 (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
862 (expr2.val()*expr2.val()),
863 (expr1.fastAccessDx(i)*expr2.val() -
864 expr2.fastAccessDx(i)*expr1.val()) /
865 (expr2.val()*expr2.val()),
866 -expr2.dx(i)*expr1.val() / (expr2.val()*expr2.val()),
867 expr1.dx(i)/expr2.val(),
868 -expr2.fastAccessDx(i)*expr1.val() / (expr2.val()*expr2.val()),
869 expr1.fastAccessDx(i)/expr2.val())
871 Atan2Op,
872 std::atan2(expr1.val(), expr2.val()),
873 bar*expr2.val()/
874 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
875 -bar*expr1.val()/
876 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
877 false,
878 false,
879 false,
880 false,
881 false,
882 false,
883 (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
884 (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
885 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
886 (-expr1.val()*expr2.dx(i)) / (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
887 (expr2.val()*expr1.dx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
888 (-expr1.val()*expr2.fastAccessDx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
889 (expr2.val()*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + expr2.val()*expr2.val()))
891 PowerOp,
892 std::pow(expr1.val(), expr2.val()),
893 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()),
894 expr1.val() == value_type(0) ? value_type(0) : value_type(bar*std::pow(expr1.val(),expr2.val())*std::log(expr1.val())),
895 false,
896 false,
897 false,
898 false,
899 false,
900 false,
901 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())),
902 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())),
903 expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.dx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
904 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())),
905 expr1.val() == value_type(0) ? value_type(0) : value_type(expr2.fastAccessDx(i)*std::log(expr1.val())*std::pow(expr1.val(),expr2.val())),
906 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())))
908 MaxOp,
909 expr1.val() >= expr2.val() ? expr1.val() : expr2.val(),
910 expr1.val() >= expr2.val() ? bar : value_type(0.),
911 expr2.val() > expr1.val() ? bar : value_type(0.),
912 expr1.isLinear() && expr2.isLinear(),
913 expr2.isLinear(),
914 expr1.isLinear(),
915 ExprT1::is_linear && ExprT2::is_linear,
916 ExprT2::is_linear,
917 ExprT1::is_linear,
918 expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
919 expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) :
920 expr2.fastAccessDx(i),
921 expr1.val() >= expr2.val() ? value_type(0) : expr2.dx(i),
922 expr1.val() >= expr2.val() ? expr1.dx(i) : value_type(0),
923 expr1.val() >= expr2.val() ? value_type(0) :
924 expr2.fastAccessDx(i),
925 expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) :
926 value_type(0))
928 MinOp,
929 expr1.val() <= expr2.val() ? expr1.val() : expr2.val(),
930 expr1.val() <= expr2.val() ? bar : value_type(0.),
931 expr2.val() < expr1.val() ? bar : value_type(0.),
932 expr1.isLinear() && expr2.isLinear(),
933 expr2.isLinear(),
934 expr1.isLinear(),
935 ExprT1::is_linear && ExprT2::is_linear,
936 ExprT2::is_linear,
937 ExprT1::is_linear,
938 expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
939 expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) :
940 expr2.fastAccessDx(i),
941 expr1.val() <= expr2.val() ? value_type(0) : expr2.dx(i),
942 expr1.val() <= expr2.val() ? expr1.dx(i) : value_type(0),
943 expr1.val() <= expr2.val() ? value_type(0) :
944 expr2.fastAccessDx(i),
945 expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) :
946 value_type(0))
947
948#undef FAD_BINARYOP_MACRO
949
950//-------------------------- Relational Operators -----------------------
951
952#define FAD_RELOP_MACRO(OP) \
953namespace Sacado { \
954 namespace ELRFad { \
955 template <typename ExprT1, typename ExprT2> \
956 SACADO_INLINE_FUNCTION \
957 bool \
958 operator OP (const Expr<ExprT1>& expr1, \
959 const Expr<ExprT2>& expr2) \
960 { \
961 return expr1.val() OP expr2.val(); \
962 } \
963 \
964 template <typename ExprT2> \
965 SACADO_INLINE_FUNCTION \
966 bool \
967 operator OP (const typename Expr<ExprT2>::value_type& a, \
968 const Expr<ExprT2>& expr2) \
969 { \
970 return a OP expr2.val(); \
971 } \
972 \
973 template <typename ExprT1> \
974 SACADO_INLINE_FUNCTION \
975 bool \
976 operator OP (const Expr<ExprT1>& expr1, \
977 const typename Expr<ExprT1>::value_type& b) \
978 { \
979 return expr1.val() OP b; \
980 } \
981 } \
982}
983
994
995#undef FAD_RELOP_MACRO
996
997namespace Sacado {
998
999 namespace ELRFad {
1000
1001 template <typename ExprT>
1003 bool operator ! (const Expr<ExprT>& expr)
1004 {
1005 return ! expr.val();
1006 }
1007
1008 } // namespace ELRFad
1009
1010} // namespace Sacado
1011
1012//-------------------------- Boolean Operators -----------------------
1013namespace Sacado {
1014
1015 namespace ELRFad {
1016
1017 template <typename ExprT>
1019 bool toBool(const Expr<ExprT>& x) {
1020 bool is_zero = (x.val() == 0.0);
1021 for (int i=0; i<x.size(); i++)
1022 is_zero = is_zero && (x.dx(i) == 0.0);
1023 return !is_zero;
1024 }
1025
1026 } // namespace Fad
1027
1028} // namespace Sacado
1029
1030#define FAD_BOOL_MACRO(OP) \
1031namespace Sacado { \
1032 namespace ELRFad { \
1033 template <typename ExprT1, typename ExprT2> \
1034 SACADO_INLINE_FUNCTION \
1035 bool \
1036 operator OP (const Expr<ExprT1>& expr1, \
1037 const Expr<ExprT2>& expr2) \
1038 { \
1039 return toBool(expr1) OP toBool(expr2); \
1040 } \
1041 \
1042 template <typename ExprT2> \
1043 SACADO_INLINE_FUNCTION \
1044 bool \
1045 operator OP (const typename Expr<ExprT2>::value_type& a, \
1046 const Expr<ExprT2>& expr2) \
1047 { \
1048 return a OP toBool(expr2); \
1049 } \
1050 \
1051 template <typename ExprT1> \
1052 SACADO_INLINE_FUNCTION \
1053 bool \
1054 operator OP (const Expr<ExprT1>& expr1, \
1055 const typename Expr<ExprT1>::value_type& b) \
1056 { \
1057 return toBool(expr1) OP b; \
1058 } \
1059 } \
1060}
1061
1064
1065#undef FAD_BOOL_MACRO
1066
1067//-------------------------- I/O Operators -----------------------
1068
1069namespace Sacado {
1070
1071 namespace ELRFad {
1072
1073 template <typename ExprT>
1074 std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
1075 os << x.val() << " [";
1076
1077 for (int i=0; i< x.size(); i++) {
1078 os << " " << x.dx(i);
1079 }
1080
1081 os << " ]";
1082 return os;
1083 }
1084
1085 } // namespace Fad
1086
1087} // namespace Sacado
1088
1089#endif // SACADO_FAD_OPS_HPP
#define SACADO_INLINE_FUNCTION
expr expr expr bar bar expr expr expr Log10Op
expr expr expr bar LogOp
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
#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)
expr expr expr ExpOp
expr expr dx(i)
#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)
expr expr SinOp
expr expr SinhOp
expr expr SqrtOp
expr expr ACosOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
expr expr ATanOp
expr expr ACoshOp
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
expr expr TanhOp
expr expr TanOp
expr expr CoshOp
expr expr ASinhOp
expr expr AbsOp
expr expr ATanhOp
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
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
std::ostream & operator<<(std::ostream &os, const GeneralFad< T, Storage > &x)
SACADO_INLINE_FUNCTION bool toBool(const Expr< T > &xx)