Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_Fad_Ops.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Sacado Package
5// Copyright (2006) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25// (etphipp@sandia.gov).
26//
27// ***********************************************************************
28//
29// The forward-mode AD classes in Sacado are a derivative work of the
30// expression template classes in the Fad package by Nicolas Di Cesare.
31// The following banner is included in the original Fad source code:
32//
33// ************ DO NOT REMOVE THIS BANNER ****************
34//
35// Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
36// http://www.ann.jussieu.fr/~dicesare
37//
38// CEMRACS 98 : C++ courses,
39// templates : new C++ techniques
40// for scientific computing
41//
42//********************************************************
43//
44// A short implementation ( not all operators and
45// functions are overloaded ) of 1st order Automatic
46// Differentiation in forward mode (FAD) using
47// EXPRESSION TEMPLATES.
48//
49//********************************************************
50// @HEADER
51
52#ifndef SACADO_FAD_OPS_HPP
53#define SACADO_FAD_OPS_HPP
54
57#include "Sacado_cmath.hpp"
58#include <ostream> // for std::ostream
59
60#define FAD_UNARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX) \
61namespace Sacado { \
62 namespace Fad { \
63 \
64 template <typename ExprT> \
65 class OP {}; \
66 \
67 template <typename ExprT> \
68 struct ExprSpec< OP<ExprT> > { \
69 typedef typename ExprSpec<ExprT>::type type; \
70 }; \
71 \
72 template <typename ExprT> \
73 class Expr< OP<ExprT>,ExprSpecDefault > { \
74 public: \
75 \
76 typedef typename ExprT::value_type value_type; \
77 typedef typename ExprT::scalar_type scalar_type; \
78 typedef typename ExprT::base_expr_type base_expr_type; \
79 \
80 SACADO_INLINE_FUNCTION \
81 explicit Expr(const ExprT& expr_) : expr(expr_) {} \
82 \
83 SACADO_INLINE_FUNCTION \
84 int size() const { return expr.size(); } \
85 \
86 SACADO_INLINE_FUNCTION \
87 bool hasFastAccess() const { return expr.hasFastAccess(); } \
88 \
89 SACADO_INLINE_FUNCTION \
90 bool isPassive() const { return expr.isPassive();} \
91 \
92 SACADO_INLINE_FUNCTION \
93 bool updateValue() const { return expr.updateValue(); } \
94 \
95 SACADO_INLINE_FUNCTION \
96 void cache() const {} \
97 \
98 SACADO_INLINE_FUNCTION \
99 value_type val() const { \
100 USING \
101 return VALUE; \
102 } \
103 \
104 SACADO_INLINE_FUNCTION \
105 value_type dx(int i) const { \
106 USING \
107 return DX; \
108 } \
109 \
110 SACADO_INLINE_FUNCTION \
111 value_type fastAccessDx(int i) const { \
112 USING \
113 return FASTACCESSDX; \
114 } \
115 \
116 protected: \
117 \
118 const ExprT& expr; \
119 }; \
120 \
121 template <typename T> \
122 SACADO_INLINE_FUNCTION \
123 Expr< OP< Expr<T> > > \
124 OPNAME (const Expr<T>& expr) \
125 { \
126 typedef OP< Expr<T> > expr_t; \
127 \
128 return Expr<expr_t>(expr); \
129 } \
130 } \
131 \
132}
133
135 UnaryPlusOp,
136 ;,
137 expr.val(),
138 expr.dx(i),
139 expr.fastAccessDx(i))
140FAD_UNARYOP_MACRO(operator-,
142 ;,
143 -expr.val(),
144 -expr.dx(i),
145 -expr.fastAccessDx(i))
148 using std::exp;,
149 exp(expr.val()),
150 exp(expr.val())*expr.dx(i),
151 exp(expr.val())*expr.fastAccessDx(i))
154 using std::log;,
155 log(expr.val()),
156 expr.dx(i)/expr.val(),
157 expr.fastAccessDx(i)/expr.val())
160 using std::log10; using std::log;,
161 log10(expr.val()),
162 expr.dx(i)/( log(value_type(10))*expr.val()),
163 expr.fastAccessDx(i) / ( log(value_type(10))*expr.val()))
166 using std::sqrt;,
167 sqrt(expr.val()),
168 expr.dx(i)/(value_type(2)* sqrt(expr.val())),
169 expr.fastAccessDx(i)/(value_type(2)* sqrt(expr.val())))
172 using std::cos; using std::sin;,
173 cos(expr.val()),
174 -expr.dx(i)* sin(expr.val()),
175 -expr.fastAccessDx(i)* sin(expr.val()))
178 using std::cos; using std::sin;,
179 sin(expr.val()),
180 expr.dx(i)* cos(expr.val()),
181 expr.fastAccessDx(i)* cos(expr.val()))
184 using std::tan;,
185 std::tan(expr.val()),
186 expr.dx(i)*
187 (value_type(1)+ tan(expr.val())* tan(expr.val())),
188 expr.fastAccessDx(i)*
189 (value_type(1)+ tan(expr.val())* tan(expr.val())))
192 using std::acos; using std::sqrt;,
193 acos(expr.val()),
194 -expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
195 -expr.fastAccessDx(i) /
196 sqrt(value_type(1)-expr.val()*expr.val()))
199 using std::asin; using std::sqrt;,
200 asin(expr.val()),
201 expr.dx(i)/ sqrt(value_type(1)-expr.val()*expr.val()),
202 expr.fastAccessDx(i) /
203 sqrt(value_type(1)-expr.val()*expr.val()))
206 using std::atan;,
207 atan(expr.val()),
208 expr.dx(i)/(value_type(1)+expr.val()*expr.val()),
209 expr.fastAccessDx(i)/(value_type(1)+expr.val()*expr.val()))
212 using std::cosh; using std::sinh;,
213 cosh(expr.val()),
214 expr.dx(i)* sinh(expr.val()),
215 expr.fastAccessDx(i)* sinh(expr.val()))
218 using std::cosh; using std::sinh;,
219 sinh(expr.val()),
220 expr.dx(i)* cosh(expr.val()),
221 expr.fastAccessDx(i)* cosh(expr.val()))
224 using std::tanh;,
225 tanh(expr.val()),
226 expr.dx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())),
227 expr.fastAccessDx(i)*(value_type(1)-tanh(expr.val())*tanh(expr.val())))
230 using std::acosh; using std::sqrt;,
231 acosh(expr.val()),
232 expr.dx(i)/ sqrt((expr.val()-value_type(1)) *
233 (expr.val()+value_type(1))),
234 expr.fastAccessDx(i)/ sqrt((expr.val()-value_type(1)) *
235 (expr.val()+value_type(1))))
238 using std::asinh; using std::sqrt;,
239 asinh(expr.val()),
240 expr.dx(i)/ sqrt(value_type(1)+expr.val()*expr.val()),
241 expr.fastAccessDx(i)/ sqrt(value_type(1)+
242 expr.val()*expr.val()))
245 using std::atanh;,
246 atanh(expr.val()),
247 expr.dx(i)/(value_type(1)-expr.val()*expr.val()),
248 expr.fastAccessDx(i)/(value_type(1)-
249 expr.val()*expr.val()))
252 using std::abs; using Sacado::if_then_else;,
253 abs(expr.val()),
254 if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
255 if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
258 using std::fabs; using Sacado::if_then_else;,
259 fabs(expr.val()),
260 if_then_else( expr.val() >= 0, expr.dx(i), value_type(-expr.dx(i)) ),
261 if_then_else( expr.val() >= 0, expr.fastAccessDx(i), value_type(-expr.fastAccessDx(i)) ) )
262#ifdef HAVE_SACADO_CXX11
264 CbrtOp,
265 using std::cbrt;,
266 cbrt(expr.val()),
267 expr.dx(i)/(value_type(3)*cbrt(expr.val()*expr.val())),
268 expr.fastAccessDx(i)/(value_type(3)*cbrt(expr.val()*expr.val())))
269#endif
270
271#undef FAD_UNARYOP_MACRO
272
273// Special handling for safe_sqrt() to provide specializations of SafeSqrtOp for
274// "simd" value types that use if_then_else(). The only reason for not using
275// if_then_else() always is to avoid evaluating the derivative if the value is
276// zero to avoid throwing FPEs.
277namespace Sacado {
278 namespace Fad {
279
280 template <typename ExprT, bool is_simd>
281 class SafeSqrtOp {};
282
283 template <typename ExprT>
284 struct ExprSpec< SafeSqrtOp<ExprT> > {
285 typedef typename ExprSpec<ExprT>::type type;
286 };
287
288 //
289 // Implementation for simd type using if_then_else()
290 //
291 template <typename ExprT>
292 class Expr< SafeSqrtOp<ExprT,true>,ExprSpecDefault > {
293 public:
294
295 typedef typename ExprT::value_type value_type;
296 typedef typename ExprT::scalar_type scalar_type;
297 typedef typename ExprT::base_expr_type base_expr_type;
298
300 explicit Expr(const ExprT& expr_) : expr(expr_) {}
301
303 int size() const { return expr.size(); }
304
306 bool hasFastAccess() const { return expr.hasFastAccess(); }
307
309 bool isPassive() const { return expr.isPassive();}
310
312 bool updateValue() const { return expr.updateValue(); }
313
315 void cache() const {}
316
318 value_type val() const {
319 using std::sqrt;
320 return sqrt(expr.val());
321 }
322
324 value_type dx(int i) const {
325 using std::sqrt; using Sacado::if_then_else;
326 return if_then_else(
327 expr.val() == value_type(0.0), value_type(0.0),
328 value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val()))));
329 }
330
332 value_type fastAccessDx(int i) const {
333 using std::sqrt; using Sacado::if_then_else;
334 return if_then_else(
335 expr.val() == value_type(0.0), value_type(0.0),
336 value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val()))));
337 }
338
339 protected:
340
341 const ExprT& expr;
342 };
343
344 //
345 // Specialization for scalar types using ternary operator
346 //
347 template <typename ExprT>
348 class Expr< SafeSqrtOp<ExprT,false>,ExprSpecDefault > {
349 public:
350
351 typedef typename ExprT::value_type value_type;
352 typedef typename ExprT::scalar_type scalar_type;
353 typedef typename ExprT::base_expr_type base_expr_type;
354
356 explicit Expr(const ExprT& expr_) : expr(expr_) {}
357
359 int size() const { return expr.size(); }
360
362 bool hasFastAccess() const { return expr.hasFastAccess(); }
363
365 bool isPassive() const { return expr.isPassive();}
366
368 bool updateValue() const { return expr.updateValue(); }
369
371 void cache() const {}
372
374 value_type val() const {
375 using std::sqrt;
376 return sqrt(expr.val());
377 }
378
380 value_type dx(int i) const {
381 using std::sqrt;
382 return expr.val() == value_type(0.0) ? value_type(0.0) :
383 value_type(expr.dx(i)/(value_type(2)*sqrt(expr.val())));
384 }
385
387 value_type fastAccessDx(int i) const {
388 using std::sqrt;
389 return expr.val() == value_type(0.0) ? value_type(0.0) :
390 value_type(expr.fastAccessDx(i)/(value_type(2)*sqrt(expr.val())));
391 }
392
393 protected:
394
395 const ExprT& expr;
396 };
397
398 template <typename T>
400 Expr< SafeSqrtOp< Expr<T> > >
401 safe_sqrt (const Expr<T>& expr)
402 {
403 typedef SafeSqrtOp< Expr<T> > expr_t;
404
405 return Expr<expr_t>(expr);
406 }
407 }
408
409}
410
411#define FAD_BINARYOP_MACRO(OPNAME,OP,USING,VALUE,DX,FASTACCESSDX,VAL_CONST_DX_1,VAL_CONST_DX_2,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
412namespace Sacado { \
413 namespace Fad { \
414 \
415 template <typename ExprT1, typename ExprT2> \
416 class OP {}; \
417 \
418 template <typename ExprT1, typename ExprT2> \
419 struct ExprSpec< OP< ExprT1, ExprT2 > > { \
420 typedef typename ExprSpec<ExprT1>::type type; \
421 }; \
422 \
423 template <typename ExprT1, typename ExprT2> \
424 class Expr< OP< ExprT1, ExprT2 >,ExprSpecDefault > { \
425 \
426 public: \
427 \
428 typedef typename ExprT1::value_type value_type_1; \
429 typedef typename ExprT2::value_type value_type_2; \
430 typedef typename Sacado::Promote<value_type_1, \
431 value_type_2>::type value_type; \
432 \
433 typedef typename ExprT1::scalar_type scalar_type_1; \
434 typedef typename ExprT2::scalar_type scalar_type_2; \
435 typedef typename Sacado::Promote<scalar_type_1, \
436 scalar_type_2>::type scalar_type; \
437 \
438 typedef typename ExprT1::base_expr_type base_expr_type_1; \
439 typedef typename ExprT2::base_expr_type base_expr_type_2; \
440 typedef typename Sacado::Promote<base_expr_type_1, \
441 base_expr_type_2>::type base_expr_type; \
442 \
443 SACADO_INLINE_FUNCTION \
444 Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
445 expr1(expr1_), expr2(expr2_) {} \
446 \
447 SACADO_INLINE_FUNCTION \
448 int size() const { \
449 int sz1 = expr1.size(), sz2 = expr2.size(); \
450 return sz1 > sz2 ? sz1 : sz2; \
451 } \
452 \
453 SACADO_INLINE_FUNCTION \
454 bool hasFastAccess() const { \
455 return expr1.hasFastAccess() && expr2.hasFastAccess(); \
456 } \
457 \
458 SACADO_INLINE_FUNCTION \
459 bool isPassive() const { \
460 return expr1.isPassive() && expr2.isPassive(); \
461 } \
462 \
463 SACADO_INLINE_FUNCTION \
464 bool updateValue() const { \
465 return expr1.updateValue() && expr2.updateValue(); \
466 } \
467 \
468 SACADO_INLINE_FUNCTION \
469 void cache() const {} \
470 \
471 SACADO_INLINE_FUNCTION \
472 const value_type val() const { \
473 USING \
474 return VALUE; \
475 } \
476 \
477 SACADO_INLINE_FUNCTION \
478 const value_type dx(int i) const { \
479 USING \
480 return DX; \
481 } \
482 \
483 SACADO_INLINE_FUNCTION \
484 const value_type fastAccessDx(int i) const { \
485 USING \
486 return FASTACCESSDX; \
487 } \
488 \
489 protected: \
490 \
491 const ExprT1& expr1; \
492 const ExprT2& expr2; \
493 \
494 }; \
495 \
496 template <typename ExprT1, typename T2> \
497 struct ExprSpec< OP< ExprT1, ConstExpr<T2> > > { \
498 typedef typename ExprSpec<ExprT1>::type type; \
499 }; \
500 \
501 template <typename ExprT1, typename T2> \
502 class Expr< OP< ExprT1, ConstExpr<T2> >,ExprSpecDefault > { \
503 \
504 public: \
505 \
506 typedef ConstExpr<T2> ConstT; \
507 typedef ConstExpr<T2> ExprT2; \
508 typedef typename ExprT1::value_type value_type_1; \
509 typedef typename ExprT2::value_type value_type_2; \
510 typedef typename Sacado::Promote<value_type_1, \
511 value_type_2>::type value_type; \
512 \
513 typedef typename ExprT1::scalar_type scalar_type_1; \
514 typedef typename ExprT2::scalar_type scalar_type_2; \
515 typedef typename Sacado::Promote<scalar_type_1, \
516 scalar_type_2>::type scalar_type; \
517 \
518 typedef typename ExprT1::base_expr_type base_expr_type_1; \
519 typedef typename ExprT2::base_expr_type base_expr_type_2; \
520 typedef typename Sacado::Promote<base_expr_type_1, \
521 base_expr_type_2>::type base_expr_type; \
522 \
523 SACADO_INLINE_FUNCTION \
524 Expr(const ExprT1& expr1_, const ConstT& c_) : \
525 expr1(expr1_), c(c_) {} \
526 \
527 SACADO_INLINE_FUNCTION \
528 int size() const { \
529 return expr1.size(); \
530 } \
531 \
532 SACADO_INLINE_FUNCTION \
533 bool hasFastAccess() const { \
534 return expr1.hasFastAccess(); \
535 } \
536 \
537 SACADO_INLINE_FUNCTION \
538 bool isPassive() const { \
539 return expr1.isPassive(); \
540 } \
541 \
542 SACADO_INLINE_FUNCTION \
543 bool updateValue() const { return expr1.updateValue(); } \
544 \
545 SACADO_INLINE_FUNCTION \
546 void cache() const {} \
547 \
548 SACADO_INLINE_FUNCTION \
549 const value_type val() const { \
550 USING \
551 return VAL_CONST_DX_2; \
552 } \
553 \
554 SACADO_INLINE_FUNCTION \
555 const value_type dx(int i) const { \
556 USING \
557 return CONST_DX_2; \
558 } \
559 \
560 SACADO_INLINE_FUNCTION \
561 const value_type fastAccessDx(int i) const { \
562 USING \
563 return CONST_FASTACCESSDX_2; \
564 } \
565 \
566 protected: \
567 \
568 const ExprT1& expr1; \
569 ConstT c; \
570 }; \
571 \
572 template <typename T1, typename ExprT2> \
573 struct ExprSpec< OP< ConstExpr<T1>, ExprT2 > > { \
574 typedef typename ExprSpec<ExprT2>::type type; \
575 }; \
576 \
577 template <typename T1, typename ExprT2> \
578 class Expr< OP< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > { \
579 \
580 public: \
581 \
582 typedef ConstExpr<T1> ConstT; \
583 typedef ConstExpr<T1> ExprT1; \
584 typedef typename ExprT1::value_type value_type_1; \
585 typedef typename ExprT2::value_type value_type_2; \
586 typedef typename Sacado::Promote<value_type_1, \
587 value_type_2>::type value_type; \
588 \
589 typedef typename ExprT1::scalar_type scalar_type_1; \
590 typedef typename ExprT2::scalar_type scalar_type_2; \
591 typedef typename Sacado::Promote<scalar_type_1, \
592 scalar_type_2>::type scalar_type; \
593 \
594 typedef typename ExprT1::base_expr_type base_expr_type_1; \
595 typedef typename ExprT2::base_expr_type base_expr_type_2; \
596 typedef typename Sacado::Promote<base_expr_type_1, \
597 base_expr_type_2>::type base_expr_type; \
598 \
599 \
600 SACADO_INLINE_FUNCTION \
601 Expr(const ConstT& c_, const ExprT2& expr2_) : \
602 c(c_), expr2(expr2_) {} \
603 \
604 SACADO_INLINE_FUNCTION \
605 int size() const { \
606 return expr2.size(); \
607 } \
608 \
609 SACADO_INLINE_FUNCTION \
610 bool hasFastAccess() const { \
611 return expr2.hasFastAccess(); \
612 } \
613 \
614 SACADO_INLINE_FUNCTION \
615 bool isPassive() const { \
616 return expr2.isPassive(); \
617 } \
618 \
619 SACADO_INLINE_FUNCTION \
620 bool updateValue() const { return expr2.updateValue(); } \
621 \
622 SACADO_INLINE_FUNCTION \
623 void cache() const {} \
624 \
625 SACADO_INLINE_FUNCTION \
626 const value_type val() const { \
627 USING \
628 return VAL_CONST_DX_1; \
629 } \
630 \
631 SACADO_INLINE_FUNCTION \
632 const value_type dx(int i) const { \
633 USING \
634 return CONST_DX_1; \
635 } \
636 \
637 SACADO_INLINE_FUNCTION \
638 const value_type fastAccessDx(int i) const { \
639 USING \
640 return CONST_FASTACCESSDX_1; \
641 } \
642 \
643 protected: \
644 \
645 ConstT c; \
646 const ExprT2& expr2; \
647 }; \
648 \
649 template <typename T1, typename T2> \
650 SACADO_INLINE_FUNCTION \
651 typename mpl::enable_if_c< \
652 ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value, \
653 Expr< OP< Expr<T1>, Expr<T2> > > \
654 >::type \
655 /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(OP)*/ \
656 OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2) \
657 { \
658 typedef OP< Expr<T1>, Expr<T2> > expr_t; \
659 \
660 return Expr<expr_t>(expr1, expr2); \
661 } \
662 \
663 template <typename T> \
664 SACADO_INLINE_FUNCTION \
665 Expr< OP< Expr<T>, Expr<T> > > \
666 OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
667 { \
668 typedef OP< Expr<T>, Expr<T> > expr_t; \
669 \
670 return Expr<expr_t>(expr1, expr2); \
671 } \
672 \
673 template <typename T> \
674 SACADO_INLINE_FUNCTION \
675 Expr< OP< ConstExpr<typename Expr<T>::value_type>, \
676 Expr<T> > > \
677 OPNAME (const typename Expr<T>::value_type& c, \
678 const Expr<T>& expr) \
679 { \
680 typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
681 typedef OP< ConstT, Expr<T> > expr_t; \
682 \
683 return Expr<expr_t>(ConstT(c), expr); \
684 } \
685 \
686 template <typename T> \
687 SACADO_INLINE_FUNCTION \
688 Expr< OP< Expr<T>, \
689 ConstExpr<typename Expr<T>::value_type> > > \
690 OPNAME (const Expr<T>& expr, \
691 const typename Expr<T>::value_type& c) \
692 { \
693 typedef ConstExpr<typename Expr<T>::value_type> ConstT; \
694 typedef OP< Expr<T>, ConstT > expr_t; \
695 \
696 return Expr<expr_t>(expr, ConstT(c)); \
697 } \
698 \
699 template <typename T> \
700 SACADO_INLINE_FUNCTION \
701 SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP) \
702 OPNAME (const typename Expr<T>::scalar_type& c, \
703 const Expr<T>& expr) \
704 { \
705 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
706 typedef OP< ConstT, Expr<T> > expr_t; \
707 \
708 return Expr<expr_t>(ConstT(c), expr); \
709 } \
710 \
711 template <typename T> \
712 SACADO_INLINE_FUNCTION \
713 SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP) \
714 OPNAME (const Expr<T>& expr, \
715 const typename Expr<T>::scalar_type& c) \
716 { \
717 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT; \
718 typedef OP< Expr<T>, ConstT > expr_t; \
719 \
720 return Expr<expr_t>(expr, ConstT(c)); \
721 } \
722 } \
723 \
724}
725
726
729 ;,
730 expr1.val() + expr2.val(),
731 expr1.dx(i) + expr2.dx(i),
732 expr1.fastAccessDx(i) + expr2.fastAccessDx(i),
733 c.val() + expr2.val(),
734 expr1.val() + c.val(),
735 expr2.dx(i),
736 expr1.dx(i),
737 expr2.fastAccessDx(i),
738 expr1.fastAccessDx(i))
739FAD_BINARYOP_MACRO(operator-,
741 ;,
742 expr1.val() - expr2.val(),
743 expr1.dx(i) - expr2.dx(i),
744 expr1.fastAccessDx(i) - expr2.fastAccessDx(i),
745 c.val() - expr2.val(),
746 expr1.val() - c.val(),
747 -expr2.dx(i),
748 expr1.dx(i),
749 -expr2.fastAccessDx(i),
750 expr1.fastAccessDx(i))
751// FAD_BINARYOP_MACRO(operator*,
752// MultiplicationOp,
753// ;,
754// expr1.val() * expr2.val(),
755// expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val(),
756// expr1.val()*expr2.fastAccessDx(i) +
757// expr1.fastAccessDx(i)*expr2.val(),
758// c.val() * expr2.val(),
759// expr1.val() * c.val(),
760// c.val()*expr2.dx(i),
761// expr1.dx(i)*c.val(),
762// c.val()*expr2.fastAccessDx(i),
763// expr1.fastAccessDx(i)*c.val())
764FAD_BINARYOP_MACRO(operator/,
766 ;,
767 expr1.val() / expr2.val(),
768 (expr1.dx(i)*expr2.val() - expr2.dx(i)*expr1.val()) /
769 (expr2.val()*expr2.val()),
770 (expr1.fastAccessDx(i)*expr2.val() -
771 expr2.fastAccessDx(i)*expr1.val()) /
772 (expr2.val()*expr2.val()),
773 c.val() / expr2.val(),
774 expr1.val() / c.val(),
775 -expr2.dx(i)*c.val() / (expr2.val()*expr2.val()),
776 expr1.dx(i)/c.val(),
777 -expr2.fastAccessDx(i)*c.val() / (expr2.val()*expr2.val()),
778 expr1.fastAccessDx(i)/c.val())
781 using std::atan2;,
782 atan2(expr1.val(), expr2.val()),
783 (expr2.val()*expr1.dx(i) - expr1.val()*expr2.dx(i))/
784 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
785 (expr2.val()*expr1.fastAccessDx(i) - expr1.val()*expr2.fastAccessDx(i))/
786 (expr1.val()*expr1.val() + expr2.val()*expr2.val()),
787 atan2(c.val(), expr2.val()),
788 atan2(expr1.val(), c.val()),
789 (-c.val()*expr2.dx(i)) / (c.val()*c.val() + expr2.val()*expr2.val()),
790 (c.val()*expr1.dx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()),
791 (-c.val()*expr2.fastAccessDx(i))/ (c.val()*c.val() + expr2.val()*expr2.val()),
792 (c.val()*expr1.fastAccessDx(i))/ (expr1.val()*expr1.val() + c.val()*c.val()))
793// FAD_BINARYOP_MACRO(pow,
794// PowerOp,
795// using std::pow; using std::log; using Sacado::if_then_else;,
796// pow(expr1.val(), expr2.val()),
797// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val())) ),
798// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type((expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val())) ),
799// pow(c.val(), expr2.val()),
800// pow(expr1.val(), c.val()),
801// if_then_else( c.val() == value_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val())) ),
802// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val())) ),
803// if_then_else( c.val() == value_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val())) ),
804// if_then_else( expr1.val() == value_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()))) )
807 using Sacado::if_then_else;,
808 if_then_else( expr1.val() >= expr2.val(), expr1.val(), expr2.val() ),
809 if_then_else( expr1.val() >= expr2.val(), expr1.dx(i), expr2.dx(i) ),
810 if_then_else( expr1.val() >= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
811 if_then_else( c.val() >= expr2.val(), value_type(c.val()), expr2.val() ),
812 if_then_else( expr1.val() >= c.val(), expr1.val(), value_type(c.val()) ),
813 if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.dx(i) ),
814 if_then_else( expr1.val() >= c.val(), expr1.dx(i), value_type(0.0) ),
815 if_then_else( c.val() >= expr2.val(), value_type(0.0), expr2.fastAccessDx(i) ),
816 if_then_else( expr1.val() >= c.val(), expr1.fastAccessDx(i), value_type(0.0) ) )
819 using Sacado::if_then_else;,
820 if_then_else( expr1.val() <= expr2.val(), expr1.val(), expr2.val() ),
821 if_then_else( expr1.val() <= expr2.val(), expr1.dx(i), expr2.dx(i) ),
822 if_then_else( expr1.val() <= expr2.val(), expr1.fastAccessDx(i), expr2.fastAccessDx(i) ),
823 if_then_else( c.val() <= expr2.val(), value_type(c.val()), expr2.val() ),
824 if_then_else( expr1.val() <= c.val(), expr1.val(), value_type(c.val()) ),
825 if_then_else( c.val() <= expr2.val(), value_type(0), expr2.dx(i) ),
826 if_then_else( expr1.val() <= c.val(), expr1.dx(i), value_type(0) ),
827 if_then_else( c.val() <= expr2.val(), value_type(0), expr2.fastAccessDx(i) ),
828 if_then_else( expr1.val() <= c.val(), expr1.fastAccessDx(i), value_type(0) ) )
829
830
831#undef FAD_BINARYOP_MACRO
832
833namespace Sacado {
834 namespace Fad {
835
836 template <typename ExprT1, typename ExprT2>
837 class MultiplicationOp {};
838
839 template <typename ExprT1, typename ExprT2>
840 struct ExprSpec< MultiplicationOp< ExprT1, ExprT2 > > {
841 typedef typename ExprSpec<ExprT1>::type type;
842 };
843
844 template <typename ExprT1, typename ExprT2>
845 class Expr< MultiplicationOp< ExprT1, ExprT2 >,ExprSpecDefault > {
846
847 public:
848
849 typedef typename ExprT1::value_type value_type_1;
850 typedef typename ExprT2::value_type value_type_2;
851 typedef typename Sacado::Promote<value_type_1,
852 value_type_2>::type value_type;
853
854 typedef typename ExprT1::scalar_type scalar_type_1;
855 typedef typename ExprT2::scalar_type scalar_type_2;
856 typedef typename Sacado::Promote<scalar_type_1,
857 scalar_type_2>::type scalar_type;
858
859 typedef typename ExprT1::base_expr_type base_expr_type_1;
860 typedef typename ExprT2::base_expr_type base_expr_type_2;
861 typedef typename Sacado::Promote<base_expr_type_1,
862 base_expr_type_2>::type base_expr_type;
863
865 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
866 expr1(expr1_), expr2(expr2_) {}
867
869 int size() const {
870 int sz1 = expr1.size(), sz2 = expr2.size();
871 return sz1 > sz2 ? sz1 : sz2;
872 }
873
875 bool hasFastAccess() const {
876 return expr1.hasFastAccess() && expr2.hasFastAccess();
877 }
878
880 bool isPassive() const {
881 return expr1.isPassive() && expr2.isPassive();
882 }
883
885 bool updateValue() const {
886 return expr1.updateValue() && expr2.updateValue();
887 }
888
890 void cache() const {}
891
893 const value_type val() const {
894 return expr1.val()*expr2.val();
895 }
896
898 const value_type dx(int i) const {
899 if (expr1.size() > 0 && expr2.size() > 0)
900 return expr1.val()*expr2.dx(i) + expr1.dx(i)*expr2.val();
901 else if (expr1.size() > 0)
902 return expr1.dx(i)*expr2.val();
903 else
904 return expr1.val()*expr2.dx(i);
905 }
906
908 const value_type fastAccessDx(int i) const {
909 return expr1.val()*expr2.fastAccessDx(i) +
910 expr1.fastAccessDx(i)*expr2.val();
911 }
912
913 protected:
914
915 const ExprT1& expr1;
916 const ExprT2& expr2;
917
918 };
919
920 template <typename ExprT1, typename T2>
921 struct ExprSpec< MultiplicationOp< ExprT1, ConstExpr<T2> > > {
922 typedef typename ExprSpec<ExprT1>::type type;
923 };
924
925 template <typename ExprT1, typename T2>
926 class Expr< MultiplicationOp< ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
927
928 public:
929
930 typedef ConstExpr<T2> ConstT;
931 typedef ConstExpr<T2> ExprT2;
932 typedef typename ExprT1::value_type value_type_1;
933 typedef typename ExprT2::value_type value_type_2;
934 typedef typename Sacado::Promote<value_type_1,
935 value_type_2>::type value_type;
936
937 typedef typename ExprT1::scalar_type scalar_type_1;
938 typedef typename ExprT2::scalar_type scalar_type_2;
939 typedef typename Sacado::Promote<scalar_type_1,
940 scalar_type_2>::type scalar_type;
941
942 typedef typename ExprT1::base_expr_type base_expr_type_1;
943 typedef typename ExprT2::base_expr_type base_expr_type_2;
944 typedef typename Sacado::Promote<base_expr_type_1,
945 base_expr_type_2>::type base_expr_type;
946
948 Expr(const ExprT1& expr1_, const ConstT& c_) :
949 expr1(expr1_), c(c_) {}
950
952 int size() const {
953 return expr1.size();
954 }
955
957 bool hasFastAccess() const {
958 return expr1.hasFastAccess();
959 }
960
962 bool isPassive() const {
963 return expr1.isPassive();
964 }
965
967 bool updateValue() const { return expr1.updateValue(); }
968
970 void cache() const {}
971
973 const value_type val() const {
974 return expr1.val()*c.val();
975 }
976
978 const value_type dx(int i) const {
979 return expr1.dx(i)*c.val();
980 }
981
983 const value_type fastAccessDx(int i) const {
984 return expr1.fastAccessDx(i)*c.val();
985 }
986
987 protected:
988
989 const ExprT1& expr1;
990 ConstT c;
991 };
992
993 template <typename T1, typename ExprT2>
994 struct ExprSpec< MultiplicationOp< ConstExpr<T1>, ExprT2 > > {
995 typedef typename ExprSpec<ExprT2>::type type;
996 };
997
998 template <typename T1, typename ExprT2>
999 class Expr< MultiplicationOp< ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
1000
1001 public:
1002
1003 typedef ConstExpr<T1> ConstT;
1004 typedef ConstExpr<T1> ExprT1;
1005 typedef typename ExprT1::value_type value_type_1;
1006 typedef typename ExprT2::value_type value_type_2;
1007 typedef typename Sacado::Promote<value_type_1,
1008 value_type_2>::type value_type;
1009
1010 typedef typename ExprT1::scalar_type scalar_type_1;
1011 typedef typename ExprT2::scalar_type scalar_type_2;
1012 typedef typename Sacado::Promote<scalar_type_1,
1013 scalar_type_2>::type scalar_type;
1014
1015 typedef typename ExprT1::base_expr_type base_expr_type_1;
1016 typedef typename ExprT2::base_expr_type base_expr_type_2;
1017 typedef typename Sacado::Promote<base_expr_type_1,
1018 base_expr_type_2>::type base_expr_type;
1019
1021 Expr(const ConstT& c_, const ExprT2& expr2_) :
1022 c(c_), expr2(expr2_) {}
1023
1025 int size() const {
1026 return expr2.size();
1027 }
1028
1030 bool hasFastAccess() const {
1031 return expr2.hasFastAccess();
1032 }
1033
1035 bool isPassive() const {
1036 return expr2.isPassive();
1037 }
1038
1040 bool updateValue() const { return expr2.updateValue(); }
1041
1043 void cache() const {}
1044
1046 const value_type val() const {
1047 return c.val()*expr2.val();
1048 }
1049
1051 const value_type dx(int i) const {
1052 return c.val()*expr2.dx(i);
1053 }
1054
1056 const value_type fastAccessDx(int i) const {
1057 return c.val()*expr2.fastAccessDx(i);
1058 }
1059
1060 protected:
1061
1062 ConstT c;
1063 const ExprT2& expr2;
1064 };
1065
1066 template <typename T1, typename T2>
1068 typename mpl::enable_if_c<
1069 ExprLevel< Expr<T1> >::value == ExprLevel< Expr<T2> >::value,
1070 Expr< MultiplicationOp< Expr<T1>, Expr<T2> > >
1071 >::type
1072 /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(MultiplicationOp)*/
1073 operator* (const Expr<T1>& expr1, const Expr<T2>& expr2)
1074 {
1075 typedef MultiplicationOp< Expr<T1>, Expr<T2> > expr_t;
1076
1077 return Expr<expr_t>(expr1, expr2);
1078 }
1079
1080 template <typename T>
1082 Expr< MultiplicationOp< Expr<T>, Expr<T> > >
1083 operator* (const Expr<T>& expr1, const Expr<T>& expr2)
1084 {
1085 typedef MultiplicationOp< Expr<T>, Expr<T> > expr_t;
1086
1087 return Expr<expr_t>(expr1, expr2);
1088 }
1089
1090 template <typename T>
1092 Expr< MultiplicationOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
1093 operator* (const typename Expr<T>::value_type& c,
1094 const Expr<T>& expr)
1095 {
1096 typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1097 typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1098
1099 return Expr<expr_t>(ConstT(c), expr);
1100 }
1101
1102 template <typename T>
1104 Expr< MultiplicationOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
1105 operator* (const Expr<T>& expr,
1106 const typename Expr<T>::value_type& c)
1107 {
1108 typedef ConstExpr<typename Expr<T>::value_type> ConstT;
1109 typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1110
1111 return Expr<expr_t>(expr, ConstT(c));
1112 }
1113
1114 template <typename T>
1117 operator* (const typename Expr<T>::scalar_type& c,
1118 const Expr<T>& expr)
1119 {
1120 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1121 typedef MultiplicationOp< ConstT, Expr<T> > expr_t;
1122
1123 return Expr<expr_t>(ConstT(c), expr);
1124 }
1125
1126 template <typename T>
1129 operator* (const Expr<T>& expr,
1130 const typename Expr<T>::scalar_type& c)
1131 {
1132 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
1133 typedef MultiplicationOp< Expr<T>, ConstT > expr_t;
1134
1135 return Expr<expr_t>(expr, ConstT(c));
1136 }
1137 }
1138}
1139
1140// Special handling for std::pow() to provide specializations of PowerOp for
1141// "simd" value types that use if_then_else(). The only reason for not using
1142// if_then_else() always is to avoid evaluating the derivative if the value is
1143// zero to avoid throwing FPEs.
1144namespace Sacado {
1145 namespace Fad {
1146
1147 template <typename ExprT1, typename ExprT2, typename Impl>
1148 class PowerOp {};
1149
1150 template <typename ExprT1, typename ExprT2>
1151 struct ExprSpec< PowerOp< ExprT1, ExprT2 > > {
1153 };
1154
1155 template <typename ExprT1, typename T2>
1156 struct ExprSpec<PowerOp< ExprT1, ConstExpr<T2> > > {
1158 };
1159
1160 template <typename T1, typename ExprT2>
1161 struct ExprSpec< PowerOp< ConstExpr<T1>, ExprT2 > > {
1163 };
1164
1165 //
1166 // Implementation for simd type using if_then_else()
1167 //
1168 template <typename ExprT1, typename ExprT2>
1169 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Simd >, ExprSpecDefault > {
1170
1171 public:
1172
1173 typedef typename ExprT1::value_type value_type_1;
1174 typedef typename ExprT2::value_type value_type_2;
1175 typedef typename Sacado::Promote<value_type_1,
1177
1178 typedef typename ExprT1::scalar_type scalar_type_1;
1179 typedef typename ExprT2::scalar_type scalar_type_2;
1180 typedef typename Sacado::Promote<scalar_type_1,
1182
1183 typedef typename ExprT1::base_expr_type base_expr_type_1;
1184 typedef typename ExprT2::base_expr_type base_expr_type_2;
1185 typedef typename Sacado::Promote<base_expr_type_1,
1187
1189 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1190 expr1(expr1_), expr2(expr2_) {}
1191
1193 int size() const {
1194 int sz1 = expr1.size(), sz2 = expr2.size();
1195 return sz1 > sz2 ? sz1 : sz2;
1196 }
1197
1199 bool hasFastAccess() const {
1200 return expr1.hasFastAccess() && expr2.hasFastAccess();
1201 }
1202
1204 bool isPassive() const {
1205 return expr1.isPassive() && expr2.isPassive();
1206 }
1207
1209 bool updateValue() const {
1210 return expr1.updateValue() && expr2.updateValue();
1211 }
1212
1214 void cache() const {}
1215
1217 const value_type val() const {
1218 using std::pow;
1219 return pow(expr1.val(), expr2.val());
1220 }
1221
1223 const value_type dx(int i) const {
1224 using std::pow; using std::log; using Sacado::if_then_else;
1225 const int sz1 = expr1.size(), sz2 = expr2.size();
1226 if (sz1 > 0 && sz2 > 0)
1227 return if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val())) );
1228 else if (sz1 > 0)
1229 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1230 // It seems less accurate and caused convergence problems in some codes
1231 return if_then_else( expr2.val() == scalar_type(1.0), expr1.dx(i), if_then_else(expr1.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val())) ));
1232 else
1233 return if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val())) );
1234 }
1235
1237 const value_type fastAccessDx(int i) const {
1238 using std::pow; using std::log; using Sacado::if_then_else;
1239 return if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type((expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val())) );
1240 }
1241
1242 protected:
1243
1244 const ExprT1& expr1;
1245 const ExprT2& expr2;
1246
1247 };
1248
1249 template <typename ExprT1, typename T2>
1250 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Simd >,
1251 ExprSpecDefault > {
1252
1253 public:
1254
1257 typedef typename ExprT1::value_type value_type_1;
1259 typedef typename Sacado::Promote<value_type_1,
1261
1262 typedef typename ExprT1::scalar_type scalar_type_1;
1264 typedef typename Sacado::Promote<scalar_type_1,
1266
1267 typedef typename ExprT1::base_expr_type base_expr_type_1;
1269 typedef typename Sacado::Promote<base_expr_type_1,
1271
1273 Expr(const ExprT1& expr1_, const ConstT& c_) :
1274 expr1(expr1_), c(c_) {}
1275
1277 int size() const {
1278 return expr1.size();
1279 }
1280
1282 bool hasFastAccess() const {
1283 return expr1.hasFastAccess();
1284 }
1285
1287 bool isPassive() const {
1288 return expr1.isPassive();
1289 }
1290
1292 bool updateValue() const { return expr1.updateValue(); }
1293
1295 void cache() const {}
1296
1298 const value_type val() const {
1299 using std::pow;
1300 return pow(expr1.val(), c.val());
1301 }
1302
1304 const value_type dx(int i) const {
1305 using std::pow; using Sacado::if_then_else;
1306 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1307 // It seems less accurate and caused convergence problems in some codes
1308 return if_then_else( c.val() == scalar_type(1.0), expr1.dx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val())) ));
1309 }
1310
1312 const value_type fastAccessDx(int i) const {
1313 using std::pow; using Sacado::if_then_else;
1314 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1315 // It seems less accurate and caused convergence problems in some codes
1316 return if_then_else( c.val() == scalar_type(1.0), expr1.fastAccessDx(i), if_then_else( expr1.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()))));
1317 }
1318
1319 protected:
1320
1321 const ExprT1& expr1;
1323 };
1324
1325 template <typename T1, typename ExprT2>
1326 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Simd >,
1327 ExprSpecDefault > {
1328
1329 public:
1330
1334 typedef typename ExprT2::value_type value_type_2;
1335 typedef typename Sacado::Promote<value_type_1,
1337
1339 typedef typename ExprT2::scalar_type scalar_type_2;
1340 typedef typename Sacado::Promote<scalar_type_1,
1342
1344 typedef typename ExprT2::base_expr_type base_expr_type_2;
1345 typedef typename Sacado::Promote<base_expr_type_1,
1347
1348
1350 Expr(const ConstT& c_, const ExprT2& expr2_) :
1351 c(c_), expr2(expr2_) {}
1352
1354 int size() const {
1355 return expr2.size();
1356 }
1357
1359 bool hasFastAccess() const {
1360 return expr2.hasFastAccess();
1361 }
1362
1364 bool isPassive() const {
1365 return expr2.isPassive();
1366 }
1367
1369 bool updateValue() const { return expr2.updateValue(); }
1370
1372 void cache() const {}
1373
1375 const value_type val() const {
1376 using std::pow;
1377 return pow(c.val(), expr2.val());
1378 }
1379
1381 const value_type dx(int i) const {
1382 using std::pow; using std::log; using Sacado::if_then_else;
1383 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val())) );
1384 }
1385
1387 const value_type fastAccessDx(int i) const {
1388 using std::pow; using std::log; using Sacado::if_then_else;
1389 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val())) );
1390 }
1391
1392 protected:
1393
1395 const ExprT2& expr2;
1396 };
1397
1398 //
1399 // Specialization for scalar types using ternary operator
1400 //
1401
1402 template <typename ExprT1, typename ExprT2>
1403 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Scalar >,
1404 ExprSpecDefault > {
1405
1406 public:
1407
1408 typedef typename ExprT1::value_type value_type_1;
1409 typedef typename ExprT2::value_type value_type_2;
1410 typedef typename Sacado::Promote<value_type_1,
1412
1413 typedef typename ExprT1::scalar_type scalar_type_1;
1414 typedef typename ExprT2::scalar_type scalar_type_2;
1415 typedef typename Sacado::Promote<scalar_type_1,
1417
1418 typedef typename ExprT1::base_expr_type base_expr_type_1;
1419 typedef typename ExprT2::base_expr_type base_expr_type_2;
1420 typedef typename Sacado::Promote<base_expr_type_1,
1422
1424 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1425 expr1(expr1_), expr2(expr2_) {}
1426
1428 int size() const {
1429 int sz1 = expr1.size(), sz2 = expr2.size();
1430 return sz1 > sz2 ? sz1 : sz2;
1431 }
1432
1434 bool hasFastAccess() const {
1435 return expr1.hasFastAccess() && expr2.hasFastAccess();
1436 }
1437
1439 bool isPassive() const {
1440 return expr1.isPassive() && expr2.isPassive();
1441 }
1442
1444 bool updateValue() const {
1445 return expr1.updateValue() && expr2.updateValue();
1446 }
1447
1449 void cache() const {}
1450
1452 const value_type val() const {
1453 using std::pow;
1454 return pow(expr1.val(), expr2.val());
1455 }
1456
1458 const value_type dx(int i) const {
1459 using std::pow; using std::log;
1460 const int sz1 = expr1.size(), sz2 = expr2.size();
1461 if (sz1 > 0 && sz2 > 0)
1462 return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type((expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val()));
1463 else if (sz1 > 0)
1464 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1465 // It seems less accurate and caused convergence problems in some codes
1466 return expr2.val() == scalar_type(1.0) ? expr1.dx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),expr2.val()));
1467 else
1468 return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val()));
1469 }
1470
1472 const value_type fastAccessDx(int i) const {
1473 using std::pow; using std::log;
1474 return expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type((expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val()));
1475 }
1476
1477 protected:
1478
1479 const ExprT1& expr1;
1480 const ExprT2& expr2;
1481
1482 };
1483
1484 template <typename ExprT1, typename T2>
1485 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Scalar >,
1486 ExprSpecDefault > {
1487
1488 public:
1489
1492 typedef typename ExprT1::value_type value_type_1;
1494 typedef typename Sacado::Promote<value_type_1,
1496
1497 typedef typename ExprT1::scalar_type scalar_type_1;
1499 typedef typename Sacado::Promote<scalar_type_1,
1501
1502 typedef typename ExprT1::base_expr_type base_expr_type_1;
1504 typedef typename Sacado::Promote<base_expr_type_1,
1506
1508 Expr(const ExprT1& expr1_, const ConstT& c_) :
1509 expr1(expr1_), c(c_) {}
1510
1512 int size() const {
1513 return expr1.size();
1514 }
1515
1517 bool hasFastAccess() const {
1518 return expr1.hasFastAccess();
1519 }
1520
1522 bool isPassive() const {
1523 return expr1.isPassive();
1524 }
1525
1527 bool updateValue() const { return expr1.updateValue(); }
1528
1530 void cache() const {}
1531
1533 const value_type val() const {
1534 using std::pow;
1535 return pow(expr1.val(), c.val());
1536 }
1537
1539 const value_type dx(int i) const {
1540 using std::pow;
1541 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1542 // It seems less accurate and caused convergence problems in some codes
1543 return c.val() == scalar_type(1.0) ? expr1.dx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)/expr1.val()*pow(expr1.val(),c.val()));
1544 }
1545
1547 const value_type fastAccessDx(int i) const {
1548 using std::pow;
1549 // Don't use formula (a(x)^b)' = b*a(x)^{b-1}*a'(x)
1550 // It seems less accurate and caused convergence problems in some codes
1551 return c.val() == scalar_type(1.0) ? expr1.fastAccessDx(i) : expr1.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)/expr1.val()*pow(expr1.val(),c.val()));
1552 }
1553
1554 protected:
1555
1556 const ExprT1& expr1;
1558 };
1559
1560 template <typename T1, typename ExprT2>
1561 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Scalar >,
1562 ExprSpecDefault > {
1563
1564 public:
1565
1569 typedef typename ExprT2::value_type value_type_2;
1570 typedef typename Sacado::Promote<value_type_1,
1572
1574 typedef typename ExprT2::scalar_type scalar_type_2;
1575 typedef typename Sacado::Promote<scalar_type_1,
1577
1579 typedef typename ExprT2::base_expr_type base_expr_type_2;
1580 typedef typename Sacado::Promote<base_expr_type_1,
1582
1583
1585 Expr(const ConstT& c_, const ExprT2& expr2_) :
1586 c(c_), expr2(expr2_) {}
1587
1589 int size() const {
1590 return expr2.size();
1591 }
1592
1594 bool hasFastAccess() const {
1595 return expr2.hasFastAccess();
1596 }
1597
1599 bool isPassive() const {
1600 return expr2.isPassive();
1601 }
1602
1604 bool updateValue() const { return expr2.updateValue(); }
1605
1607 void cache() const {}
1608
1610 const value_type val() const {
1611 using std::pow;
1612 return pow(c.val(), expr2.val());
1613 }
1614
1616 const value_type dx(int i) const {
1617 using std::pow; using std::log;
1618 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val()));
1619 }
1620
1622 const value_type fastAccessDx(int i) const {
1623 using std::pow; using std::log;
1624 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val()));
1625 }
1626
1627 protected:
1628
1630 const ExprT2& expr2;
1631 };
1632
1633 //
1634 // Specialization for nested derivatives. This version does not use
1635 // if_then_else/ternary-operator on the base so that nested derivatives
1636 // are correct.
1637 //
1638 template <typename ExprT1, typename ExprT2>
1639 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::Nested >,
1640 ExprSpecDefault > {
1641
1642 public:
1643
1644 typedef typename ExprT1::value_type value_type_1;
1645 typedef typename ExprT2::value_type value_type_2;
1646 typedef typename Sacado::Promote<value_type_1,
1648
1649 typedef typename ExprT1::scalar_type scalar_type_1;
1650 typedef typename ExprT2::scalar_type scalar_type_2;
1651 typedef typename Sacado::Promote<scalar_type_1,
1653
1654 typedef typename ExprT1::base_expr_type base_expr_type_1;
1655 typedef typename ExprT2::base_expr_type base_expr_type_2;
1656 typedef typename Sacado::Promote<base_expr_type_1,
1658
1660 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1661 expr1(expr1_), expr2(expr2_) {}
1662
1664 int size() const {
1665 int sz1 = expr1.size(), sz2 = expr2.size();
1666 return sz1 > sz2 ? sz1 : sz2;
1667 }
1668
1670 bool hasFastAccess() const {
1671 return expr1.hasFastAccess() && expr2.hasFastAccess();
1672 }
1673
1675 bool isPassive() const {
1676 return expr1.isPassive() && expr2.isPassive();
1677 }
1678
1680 bool updateValue() const {
1681 return expr1.updateValue() && expr2.updateValue();
1682 }
1683
1685 void cache() const {}
1686
1688 const value_type val() const {
1689 using std::pow;
1690 return pow(expr1.val(), expr2.val());
1691 }
1692
1694 const value_type dx(int i) const {
1695 using std::pow; using std::log;
1696 const int sz1 = expr1.size(), sz2 = expr2.size();
1697 if (sz1 > 0 && sz2 > 0)
1698 return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1699 else if (sz1 > 0)
1700 return expr2.val() == scalar_type(0.0) ? value_type(0.0) : value_type((expr2.val()*expr1.dx(i))*pow(expr1.val(),expr2.val()-scalar_type(1.0)));
1701 else
1702 return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1703 }
1704
1706 const value_type fastAccessDx(int i) const {
1707 using std::pow; using std::log;
1708 return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1709 }
1710
1711 protected:
1712
1713 const ExprT1& expr1;
1714 const ExprT2& expr2;
1715
1716 };
1717
1718 template <typename ExprT1, typename T2>
1719 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::Nested >,
1720 ExprSpecDefault > {
1721
1722 public:
1723
1726 typedef typename ExprT1::value_type value_type_1;
1728 typedef typename Sacado::Promote<value_type_1,
1730
1731 typedef typename ExprT1::scalar_type scalar_type_1;
1733 typedef typename Sacado::Promote<scalar_type_1,
1735
1736 typedef typename ExprT1::base_expr_type base_expr_type_1;
1738 typedef typename Sacado::Promote<base_expr_type_1,
1740
1742 Expr(const ExprT1& expr1_, const ConstT& c_) :
1743 expr1(expr1_), c(c_) {}
1744
1746 int size() const {
1747 return expr1.size();
1748 }
1749
1751 bool hasFastAccess() const {
1752 return expr1.hasFastAccess();
1753 }
1754
1756 bool isPassive() const {
1757 return expr1.isPassive();
1758 }
1759
1761 bool updateValue() const { return expr1.updateValue(); }
1762
1764 void cache() const {}
1765
1767 const value_type val() const {
1768 using std::pow;
1769 return pow(expr1.val(), c.val());
1770 }
1771
1773 const value_type dx(int i) const {
1774 using std::pow;
1775 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
1776 }
1777
1779 const value_type fastAccessDx(int i) const {
1780 using std::pow;
1781 return c.val() == scalar_type(0.0) ? value_type(0.0) : value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0)));
1782 }
1783
1784 protected:
1785
1786 const ExprT1& expr1;
1788 };
1789
1790 template <typename T1, typename ExprT2>
1791 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::Nested >,
1792 ExprSpecDefault > {
1793
1794 public:
1795
1799 typedef typename ExprT2::value_type value_type_2;
1800 typedef typename Sacado::Promote<value_type_1,
1802
1804 typedef typename ExprT2::scalar_type scalar_type_2;
1805 typedef typename Sacado::Promote<scalar_type_1,
1807
1809 typedef typename ExprT2::base_expr_type base_expr_type_2;
1810 typedef typename Sacado::Promote<base_expr_type_1,
1812
1813
1815 Expr(const ConstT& c_, const ExprT2& expr2_) :
1816 c(c_), expr2(expr2_) {}
1817
1819 int size() const {
1820 return expr2.size();
1821 }
1822
1824 bool hasFastAccess() const {
1825 return expr2.hasFastAccess();
1826 }
1827
1829 bool isPassive() const {
1830 return expr2.isPassive();
1831 }
1832
1834 bool updateValue() const { return expr2.updateValue(); }
1835
1837 void cache() const {}
1838
1840 const value_type val() const {
1841 using std::pow;
1842 return pow(c.val(), expr2.val());
1843 }
1844
1846 const value_type dx(int i) const {
1847 using std::pow; using std::log;
1848 return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
1849 }
1850
1852 const value_type fastAccessDx(int i) const {
1853 using std::pow; using std::log;
1854 return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
1855 }
1856
1857 protected:
1858
1860 const ExprT2& expr2;
1861 };
1862
1863 //
1864 // Specialization for nested derivatives. This version does not use
1865 // if_then_else/ternary-operator on the base so that nested derivatives
1866 // are correct.
1867 //
1868 template <typename ExprT1, typename ExprT2>
1869 class Expr< PowerOp< ExprT1, ExprT2, PowerImpl::NestedSimd >,
1870 ExprSpecDefault > {
1871
1872 public:
1873
1874 typedef typename ExprT1::value_type value_type_1;
1875 typedef typename ExprT2::value_type value_type_2;
1876 typedef typename Sacado::Promote<value_type_1,
1878
1879 typedef typename ExprT1::scalar_type scalar_type_1;
1880 typedef typename ExprT2::scalar_type scalar_type_2;
1881 typedef typename Sacado::Promote<scalar_type_1,
1883
1884 typedef typename ExprT1::base_expr_type base_expr_type_1;
1885 typedef typename ExprT2::base_expr_type base_expr_type_2;
1886 typedef typename Sacado::Promote<base_expr_type_1,
1888
1890 Expr(const ExprT1& expr1_, const ExprT2& expr2_) :
1891 expr1(expr1_), expr2(expr2_) {}
1892
1894 int size() const {
1895 int sz1 = expr1.size(), sz2 = expr2.size();
1896 return sz1 > sz2 ? sz1 : sz2;
1897 }
1898
1900 bool hasFastAccess() const {
1901 return expr1.hasFastAccess() && expr2.hasFastAccess();
1902 }
1903
1905 bool isPassive() const {
1906 return expr1.isPassive() && expr2.isPassive();
1907 }
1908
1910 bool updateValue() const {
1911 return expr1.updateValue() && expr2.updateValue();
1912 }
1913
1915 void cache() const {}
1916
1918 const value_type val() const {
1919 using std::pow;
1920 return pow(expr1.val(), expr2.val());
1921 }
1922
1924 const value_type dx(int i) const {
1925 using std::pow; using std::log; using Sacado::if_then_else;
1926 const int sz1 = expr1.size(), sz2 = expr2.size();
1927 if (sz1 > 0 && sz2 > 0)
1928 return (expr2.dx(i)*log(expr1.val())+expr2.val()*expr1.dx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1929 else if (sz1 > 0)
1930 return if_then_else( expr2.val() == scalar_type(0.0), value_type(0.0), value_type((expr2.val()*expr1.dx(i))*pow(expr1.val(),expr2.val()-scalar_type(1.0))));
1931 else
1932 return expr2.dx(i)*log(expr1.val())*pow(expr1.val(),expr2.val());
1933 }
1934
1936 const value_type fastAccessDx(int i) const {
1937 using std::pow; using std::log;
1938 return (expr2.fastAccessDx(i)*log(expr1.val())+expr2.val()*expr1.fastAccessDx(i)/expr1.val())*pow(expr1.val(),expr2.val());
1939 }
1940
1941 protected:
1942
1943 const ExprT1& expr1;
1944 const ExprT2& expr2;
1945
1946 };
1947
1948 template <typename ExprT1, typename T2>
1949 class Expr< PowerOp< ExprT1, ConstExpr<T2>, PowerImpl::NestedSimd >,
1950 ExprSpecDefault > {
1951
1952 public:
1953
1956 typedef typename ExprT1::value_type value_type_1;
1958 typedef typename Sacado::Promote<value_type_1,
1960
1961 typedef typename ExprT1::scalar_type scalar_type_1;
1963 typedef typename Sacado::Promote<scalar_type_1,
1965
1966 typedef typename ExprT1::base_expr_type base_expr_type_1;
1968 typedef typename Sacado::Promote<base_expr_type_1,
1970
1972 Expr(const ExprT1& expr1_, const ConstT& c_) :
1973 expr1(expr1_), c(c_) {}
1974
1976 int size() const {
1977 return expr1.size();
1978 }
1979
1981 bool hasFastAccess() const {
1982 return expr1.hasFastAccess();
1983 }
1984
1986 bool isPassive() const {
1987 return expr1.isPassive();
1988 }
1989
1991 bool updateValue() const { return expr1.updateValue(); }
1992
1994 void cache() const {}
1995
1997 const value_type val() const {
1998 using std::pow;
1999 return pow(expr1.val(), c.val());
2000 }
2001
2003 const value_type dx(int i) const {
2004 using std::pow; using Sacado::if_then_else;
2005 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.dx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))));
2006 }
2007
2009 const value_type fastAccessDx(int i) const {
2010 using std::pow; using Sacado::if_then_else;
2011 return if_then_else( c.val() == scalar_type(0.0), value_type(0.0), value_type(c.val()*expr1.fastAccessDx(i)*pow(expr1.val(),c.val()-scalar_type(1.0))));
2012 }
2013
2014 protected:
2015
2016 const ExprT1& expr1;
2018 };
2019
2020 template <typename T1, typename ExprT2>
2021 class Expr< PowerOp< ConstExpr<T1>, ExprT2, PowerImpl::NestedSimd >,
2022 ExprSpecDefault > {
2023
2024 public:
2025
2029 typedef typename ExprT2::value_type value_type_2;
2030 typedef typename Sacado::Promote<value_type_1,
2032
2034 typedef typename ExprT2::scalar_type scalar_type_2;
2035 typedef typename Sacado::Promote<scalar_type_1,
2037
2039 typedef typename ExprT2::base_expr_type base_expr_type_2;
2040 typedef typename Sacado::Promote<base_expr_type_1,
2042
2043
2045 Expr(const ConstT& c_, const ExprT2& expr2_) :
2046 c(c_), expr2(expr2_) {}
2047
2049 int size() const {
2050 return expr2.size();
2051 }
2052
2054 bool hasFastAccess() const {
2055 return expr2.hasFastAccess();
2056 }
2057
2059 bool isPassive() const {
2060 return expr2.isPassive();
2061 }
2062
2064 bool updateValue() const { return expr2.updateValue(); }
2065
2067 void cache() const {}
2068
2070 const value_type val() const {
2071 using std::pow;
2072 return pow(c.val(), expr2.val());
2073 }
2074
2076 const value_type dx(int i) const {
2077 using std::pow; using std::log;
2078 return expr2.dx(i)*log(c.val())*pow(c.val(),expr2.val());
2079 }
2080
2082 const value_type fastAccessDx(int i) const {
2083 using std::pow; using std::log;
2084 return expr2.fastAccessDx(i)*log(c.val())*pow(c.val(),expr2.val());
2085 }
2086
2087 protected:
2088
2090 const ExprT2& expr2;
2091 };
2092
2093 template <typename T1, typename T2>
2095 typename mpl::enable_if_c<
2098 >::type
2099 /*SACADO_FAD_OP_ENABLE_EXPR_EXPR(PowerOp)*/
2100 pow (const Expr<T1>& expr1, const Expr<T2>& expr2)
2101 {
2103
2104 return Expr<expr_t>(expr1, expr2);
2105 }
2106
2107 template <typename T>
2109 Expr< PowerOp< Expr<T>, Expr<T> > >
2110 pow (const Expr<T>& expr1, const Expr<T>& expr2)
2111 {
2112 typedef PowerOp< Expr<T>, Expr<T> > expr_t;
2113
2114 return Expr<expr_t>(expr1, expr2);
2115 }
2116
2117 template <typename T>
2119 Expr< PowerOp< ConstExpr<typename Expr<T>::value_type>, Expr<T> > >
2120 pow (const typename Expr<T>::value_type& c,
2121 const Expr<T>& expr)
2122 {
2125
2126 return Expr<expr_t>(ConstT(c), expr);
2127 }
2128
2129 template <typename T>
2131 Expr< PowerOp< Expr<T>, ConstExpr<typename Expr<T>::value_type> > >
2132 pow (const Expr<T>& expr,
2133 const typename Expr<T>::value_type& c)
2134 {
2136 typedef PowerOp< Expr<T>, ConstT > expr_t;
2137
2138 return Expr<expr_t>(expr, ConstT(c));
2139 }
2140
2141 template <typename T>
2144 pow (const typename Expr<T>::scalar_type& c,
2145 const Expr<T>& expr)
2146 {
2149
2150 return Expr<expr_t>(ConstT(c), expr);
2151 }
2152
2153 template <typename T>
2156 pow (const Expr<T>& expr,
2157 const typename Expr<T>::scalar_type& c)
2158 {
2160 typedef PowerOp< Expr<T>, ConstT > expr_t;
2161
2162 return Expr<expr_t>(expr, ConstT(c));
2163 }
2164 }
2165}
2166
2167//--------------------------if_then_else operator -----------------------
2168// Can't use the above macros because it is a ternary operator (sort of).
2169// Also, relies on C++11
2170
2171#ifdef HAVE_SACADO_CXX11
2172
2173namespace Sacado {
2174 namespace Fad {
2175
2176 template <typename CondT, typename ExprT1, typename ExprT2>
2177 class IfThenElseOp {};
2178
2179 template <typename CondT, typename ExprT1, typename ExprT2>
2180 struct ExprSpec< IfThenElseOp< CondT, ExprT1, ExprT2 > > {
2181 typedef typename ExprSpec<ExprT1>::type type;
2182 };
2183
2184 template <typename CondT, typename ExprT1, typename ExprT2>
2185 class Expr< IfThenElseOp< CondT, ExprT1, ExprT2 >,ExprSpecDefault > {
2186
2187 public:
2188
2189 typedef typename ExprT1::value_type value_type_1;
2190 typedef typename ExprT2::value_type value_type_2;
2191 typedef typename Sacado::Promote<value_type_1,
2192 value_type_2>::type value_type;
2193
2194 typedef typename ExprT1::scalar_type scalar_type_1;
2195 typedef typename ExprT2::scalar_type scalar_type_2;
2196 typedef typename Sacado::Promote<scalar_type_1,
2197 scalar_type_2>::type scalar_type;
2198
2199 typedef typename ExprT1::base_expr_type base_expr_type_1;
2200 typedef typename ExprT2::base_expr_type base_expr_type_2;
2201 typedef typename Sacado::Promote<base_expr_type_1,
2202 base_expr_type_2>::type base_expr_type;
2203
2205 Expr(const CondT& cond_, const ExprT1& expr1_, const ExprT2& expr2_) :
2206 cond(cond_), expr1(expr1_), expr2(expr2_) {}
2207
2209 int size() const {
2210 int sz1 = expr1.size(), sz2 = expr2.size();
2211 return sz1 > sz2 ? sz1 : sz2;
2212 }
2213
2215 bool hasFastAccess() const {
2216 return expr1.hasFastAccess() && expr2.hasFastAccess();
2217 }
2218
2220 bool isPassive() const {
2221 return expr1.isPassive() && expr2.isPassive();
2222 }
2223
2225 bool updateValue() const {
2226 return expr1.updateValue() && expr2.updateValue();
2227 }
2228
2230 void cache() const {}
2231
2233 const value_type val() const {
2235 return if_then_else( cond, expr1.val(), expr2.val() );
2236 }
2237
2239 const value_type dx(int i) const {
2241 return if_then_else( cond, expr1.dx(i), expr2.dx(i) );
2242 }
2243
2245 const value_type fastAccessDx(int i) const {
2247 return if_then_else( cond, expr1.fastAccessDx(i), expr2.fastAccessDx(i) );
2248 }
2249
2250 protected:
2251
2252 const CondT& cond;
2253 const ExprT1& expr1;
2254 const ExprT2& expr2;
2255
2256 };
2257
2258 template <typename CondT, typename ExprT1, typename T2>
2259 struct ExprSpec< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> > > {
2260 typedef typename ExprSpec<ExprT1>::type type;
2261 };
2262
2263 template <typename CondT, typename ExprT1, typename T2>
2264 class Expr< IfThenElseOp< CondT, ExprT1, ConstExpr<T2> >,ExprSpecDefault > {
2265
2266 public:
2267
2268 typedef ConstExpr<T2> ConstT;
2269 typedef ConstExpr<T2> ExprT2;
2270 typedef typename ExprT1::value_type value_type_1;
2271 typedef typename ExprT2::value_type value_type_2;
2272 typedef typename Sacado::Promote<value_type_1,
2273 value_type_2>::type value_type;
2274
2275 typedef typename ExprT1::scalar_type scalar_type_1;
2276 typedef typename ExprT2::scalar_type scalar_type_2;
2277 typedef typename Sacado::Promote<scalar_type_1,
2278 scalar_type_2>::type scalar_type;
2279
2280 typedef typename ExprT1::base_expr_type base_expr_type_1;
2281 typedef typename ExprT2::base_expr_type base_expr_type_2;
2282 typedef typename Sacado::Promote<base_expr_type_1,
2283 base_expr_type_2>::type base_expr_type;
2284
2286 Expr(const CondT& cond_, const ExprT1& expr1_, const ConstT& c_) :
2287 cond(cond_), expr1(expr1_), c(c_) {}
2288
2290 int size() const {
2291 return expr1.size();
2292 }
2293
2295 bool hasFastAccess() const {
2296 return expr1.hasFastAccess();
2297 }
2298
2300 bool isPassive() const {
2301 return expr1.isPassive();
2302 }
2303
2305 bool updateValue() const { return expr1.updateValue(); }
2306
2308 void cache() const {}
2309
2311 const value_type val() const {
2313 return if_then_else( cond, expr1.val(), c.val() );
2314 }
2315
2317 const value_type dx(int i) const {
2319 return if_then_else( cond, expr1.dx(i), value_type(0.0) );
2320 }
2321
2323 const value_type fastAccessDx(int i) const {
2325 return if_then_else( cond, expr1.fastAccessDx(i), value_type(0.0) );
2326 }
2327
2328 protected:
2329
2330 const CondT& cond;
2331 const ExprT1& expr1;
2332 ConstT c;
2333 };
2334
2335 template <typename CondT, typename T1, typename ExprT2>
2336 struct ExprSpec< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 > > {
2337 typedef typename ExprSpec<ExprT2>::type type;
2338 };
2339
2340 template <typename CondT, typename T1, typename ExprT2>
2341 class Expr< IfThenElseOp< CondT, ConstExpr<T1>, ExprT2 >,ExprSpecDefault > {
2342
2343 public:
2344
2345 typedef ConstExpr<T1> ConstT;
2346 typedef ConstExpr<T1> ExprT1;
2347 typedef typename ExprT1::value_type value_type_1;
2348 typedef typename ExprT2::value_type value_type_2;
2349 typedef typename Sacado::Promote<value_type_1,
2350 value_type_2>::type value_type;
2351
2352 typedef typename ExprT1::scalar_type scalar_type_1;
2353 typedef typename ExprT2::scalar_type scalar_type_2;
2354 typedef typename Sacado::Promote<scalar_type_1,
2355 scalar_type_2>::type scalar_type;
2356
2357 typedef typename ExprT1::base_expr_type base_expr_type_1;
2358 typedef typename ExprT2::base_expr_type base_expr_type_2;
2359 typedef typename Sacado::Promote<base_expr_type_1,
2360 base_expr_type_2>::type base_expr_type;
2361
2363 Expr(const CondT& cond_, const ConstT& c_, const ExprT2& expr2_) :
2364 cond(cond_), c(c_), expr2(expr2_) {}
2365
2367 int size() const {
2368 return expr2.size();
2369 }
2370
2372 bool hasFastAccess() const {
2373 return expr2.hasFastAccess();
2374 }
2375
2377 bool isPassive() const {
2378 return expr2.isPassive();
2379 }
2380
2382 bool updateValue() const { return expr2.updateValue(); }
2383
2385 void cache() const {}
2386
2388 const value_type val() const {
2390 return if_then_else( cond, c.val(), expr2.val() );
2391 }
2392
2394 const value_type dx(int i) const {
2396 return if_then_else( cond, value_type(0.0), expr2.dx(i) );
2397 }
2398
2400 const value_type fastAccessDx(int i) const {
2402 return if_then_else( cond, value_type(0.0), expr2.fastAccessDx(i) );
2403 }
2404
2405 protected:
2406
2407 const CondT& cond;
2408 ConstT c;
2409 const ExprT2& expr2;
2410 };
2411
2412 template <typename CondT, typename T1, typename T2>
2414 typename mpl::enable_if_c< IsFadExpr<T1>::value && IsFadExpr<T2>::value &&
2416 Expr< IfThenElseOp< CondT, T1, T2 > >
2417 >::type
2418 if_then_else (const CondT& cond, const T1& expr1, const T2& expr2)
2419 {
2420 typedef IfThenElseOp< CondT, T1, T2 > expr_t;
2421
2422 return Expr<expr_t>(cond, expr1, expr2);
2423 }
2424
2425 template <typename CondT, typename T>
2427 Expr< IfThenElseOp< CondT, Expr<T>, Expr<T> > >
2428 if_then_else (const CondT& cond, const Expr<T>& expr1, const Expr<T>& expr2)
2429 {
2430 typedef IfThenElseOp< CondT, Expr<T>, Expr<T> > expr_t;
2431
2432 return Expr<expr_t>(cond, expr1, expr2);
2433 }
2434
2435 template <typename CondT, typename T>
2437 Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::value_type>,
2438 Expr<T> > >
2439 if_then_else (const CondT& cond, const typename Expr<T>::value_type& c,
2440 const Expr<T>& expr)
2441 {
2442 typedef ConstExpr<typename Expr<T>::value_type> ConstT;
2443 typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2444
2445 return Expr<expr_t>(cond, ConstT(c), expr);
2446 }
2447
2448 template <typename CondT, typename T>
2450 Expr< IfThenElseOp< CondT, Expr<T>,
2451 ConstExpr<typename Expr<T>::value_type> > >
2452 if_then_else (const CondT& cond, const Expr<T>& expr,
2453 const typename Expr<T>::value_type& c)
2454 {
2455 typedef ConstExpr<typename Expr<T>::value_type> ConstT;
2456 typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2457
2458 return Expr<expr_t>(cond, expr, ConstT(c));
2459 }
2460
2461 template <typename CondT, typename T>
2463 typename mpl::disable_if<
2464 mpl::is_same< typename Expr<T>::value_type,
2465 typename Expr<T>::scalar_type>,
2466 Expr< IfThenElseOp< CondT, ConstExpr<typename Expr<T>::scalar_type>,
2467 Expr<T> > >
2468 >::type
2469 if_then_else (const CondT& cond, const typename Expr<T>::scalar_type& c,
2470 const Expr<T>& expr)
2471 {
2472 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
2473 typedef IfThenElseOp< CondT, ConstT, Expr<T> > expr_t;
2474
2475 return Expr<expr_t>(cond, ConstT(c), expr);
2476 }
2477
2478 template <typename CondT, typename T>
2480 typename mpl::disable_if<
2481 mpl::is_same< typename Expr<T>::value_type,
2482 typename Expr<T>::scalar_type>,
2483 Expr< IfThenElseOp< CondT, Expr<T>,
2484 ConstExpr<typename Expr<T>::scalar_type> > >
2485 >::type
2486 if_then_else (const CondT& cond, const Expr<T>& expr,
2487 const typename Expr<T>::scalar_type& c)
2488 {
2489 typedef ConstExpr<typename Expr<T>::scalar_type> ConstT;
2490 typedef IfThenElseOp< CondT, Expr<T>, ConstT > expr_t;
2491
2492 return Expr<expr_t>(cond, expr, ConstT(c));
2493 }
2494 }
2495}
2496
2497#endif
2498
2499//-------------------------- Relational Operators -----------------------
2500
2501#ifdef HAVE_SACADO_CXX11
2502
2503namespace Sacado {
2504 namespace Fad {
2505 template <typename T1, typename T2 = T1>
2506 struct ConditionalReturnType {
2507 typedef decltype( std::declval<T1>() == std::declval<T2>() ) type;
2508 };
2509 }
2510}
2511
2512#define FAD_RELOP_MACRO(OP) \
2513namespace Sacado { \
2514 namespace Fad { \
2515 template <typename ExprT1, typename ExprT2> \
2516 SACADO_INLINE_FUNCTION \
2517 typename ConditionalReturnType<typename Expr<ExprT1>::value_type, \
2518 typename Expr<ExprT2>::value_type>::type \
2519 operator OP (const Expr<ExprT1>& expr1, \
2520 const Expr<ExprT2>& expr2) \
2521 { \
2522 return expr1.val() OP expr2.val(); \
2523 } \
2524 \
2525 template <typename ExprT2> \
2526 SACADO_INLINE_FUNCTION \
2527 typename ConditionalReturnType<typename Expr<ExprT2>::value_type>::type \
2528 operator OP (const typename Expr<ExprT2>::value_type& a, \
2529 const Expr<ExprT2>& expr2) \
2530 { \
2531 return a OP expr2.val(); \
2532 } \
2533 \
2534 template <typename ExprT1> \
2535 SACADO_INLINE_FUNCTION \
2536 typename ConditionalReturnType<typename Expr<ExprT1>::value_type>::type \
2537 operator OP (const Expr<ExprT1>& expr1, \
2538 const typename Expr<ExprT1>::value_type& b) \
2539 { \
2540 return expr1.val() OP b; \
2541 } \
2542 } \
2543}
2544
2545#else
2546
2547#define FAD_RELOP_MACRO(OP) \
2548namespace Sacado { \
2549 namespace Fad { \
2550 template <typename ExprT1, typename ExprT2> \
2551 SACADO_INLINE_FUNCTION \
2552 bool \
2553 operator OP (const Expr<ExprT1>& expr1, \
2554 const Expr<ExprT2>& expr2) \
2555 { \
2556 return expr1.val() OP expr2.val(); \
2557 } \
2558 \
2559 template <typename ExprT2> \
2560 SACADO_INLINE_FUNCTION \
2561 bool \
2562 operator OP (const typename Expr<ExprT2>::value_type& a, \
2563 const Expr<ExprT2>& expr2) \
2564 { \
2565 return a OP expr2.val(); \
2566 } \
2567 \
2568 template <typename ExprT1> \
2569 SACADO_INLINE_FUNCTION \
2570 bool \
2571 operator OP (const Expr<ExprT1>& expr1, \
2572 const typename Expr<ExprT1>::value_type& b) \
2573 { \
2574 return expr1.val() OP b; \
2575 } \
2576 } \
2577}
2578
2579#endif
2580
2587FAD_RELOP_MACRO(<<=)
2588FAD_RELOP_MACRO(>>=)
2591
2592#undef FAD_RELOP_MACRO
2593
2594namespace Sacado {
2595
2596 namespace Fad {
2597
2598 template <typename ExprT>
2600 bool operator ! (const Expr<ExprT>& expr)
2601 {
2602 return ! expr.val();
2603 }
2604
2605 } // namespace Fad
2606
2607} // namespace Sacado
2608
2609//-------------------------- Boolean Operators -----------------------
2610namespace Sacado {
2611
2612 namespace Fad {
2613
2614 template <typename ExprT>
2616 bool toBool(const Expr<ExprT>& x) {
2617 bool is_zero = (x.val() == 0.0);
2618 for (int i=0; i<x.size(); i++)
2619 is_zero = is_zero && (x.dx(i) == 0.0);
2620 return !is_zero;
2621 }
2622
2623 } // namespace Fad
2624
2625} // namespace Sacado
2626
2627#define FAD_BOOL_MACRO(OP) \
2628namespace Sacado { \
2629 namespace Fad { \
2630 template <typename ExprT1, typename ExprT2> \
2631 SACADO_INLINE_FUNCTION \
2632 bool \
2633 operator OP (const Expr<ExprT1>& expr1, \
2634 const Expr<ExprT2>& expr2) \
2635 { \
2636 return toBool(expr1) OP toBool(expr2); \
2637 } \
2638 \
2639 template <typename ExprT2> \
2640 SACADO_INLINE_FUNCTION \
2641 bool \
2642 operator OP (const typename Expr<ExprT2>::value_type& a, \
2643 const Expr<ExprT2>& expr2) \
2644 { \
2645 return a OP toBool(expr2); \
2646 } \
2647 \
2648 template <typename ExprT1> \
2649 SACADO_INLINE_FUNCTION \
2650 bool \
2651 operator OP (const Expr<ExprT1>& expr1, \
2652 const typename Expr<ExprT1>::value_type& b) \
2653 { \
2654 return toBool(expr1) OP b; \
2655 } \
2656 } \
2657}
2658
2661
2662#undef FAD_BOOL_MACRO
2663
2664//-------------------------- I/O Operators -----------------------
2665
2666namespace Sacado {
2667
2668 namespace Fad {
2669
2670 template <typename ExprT>
2671 std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
2672 os << x.val() << " [";
2673
2674 for (int i=0; i< x.size(); i++) {
2675 os << " " << x.dx(i);
2676 }
2677
2678 os << " ]";
2679 return os;
2680 }
2681
2682 } // namespace Fad
2683
2684} // namespace Sacado
2685
2686#endif // SACADO_FAD_OPS_HPP
#define SACADO_INLINE_FUNCTION
expr expr expr bar false
expr true
expr expr dx(i)
cbrt(expr.val())
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
log(expr.val())
#define FAD_BINARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
expr expr SinOp
expr expr SinhOp
expr expr SqrtOp
expr1 expr1 expr1 c expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr1 expr2 expr1 expr2 expr1 Atan2Op
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
expr1 expr1 expr1 c expr1 expr2 expr1 expr2 expr1 DivisionOp
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
expr expr ACosOp
expr expr ATanOp
expr expr ACoshOp
expr expr ASinOp
log10(expr.val())
expr expr TanhOp
expr expr TanOp
expr expr CoshOp
#define FAD_RELOP_MACRO(OP)
expr expr ASinhOp
exp(expr.val())
expr expr Log10Op
#define FAD_BOOL_MACRO(OP)
expr expr expr ExpOp
expr2 expr2 c c c c MaxOp
expr expr dx(i)
expr expr AbsOp
expr expr ATanhOp
SubtractionOp
expr expr CosOp
UnaryMinusOp
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
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
#define SACADO_FAD_OP_ENABLE_SCALAR_EXPR(OP)
#define SACADO_FAD_OP_ENABLE_EXPR_SCALAR(OP)
#define T1(r, f)
#define T2(r, f)
Fad specializations for Teuchos::BLAS wrappers.
Constant expression template.
T2 base_expr_type
Typename of base-expressions.
ScalarType< value_type >::type scalar_type
Typename of scalar's (which may be different from ConstT)
T2 value_type
Typename of argument values.
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
SACADO_INLINE_FUNCTION Expr(const ExprT1 &expr1_, const ExprT2 &expr2_)
Sacado::Promote< base_expr_type_1, base_expr_type_2 >::type base_expr_type
Wrapper for a generic expression template.
SimpleFad< ValueT > log(const SimpleFad< ValueT > &a)
std::ostream & operator<<(std::ostream &os, const Expr< ExprT > &x)
SimpleFad< ValueT > sqrt(const SimpleFad< ValueT > &a)
SimpleFad< ValueT > operator*(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
SACADO_INLINE_FUNCTION Expr< SafeSqrtOp< Expr< T > > > safe_sqrt(const Expr< T > &)
SACADO_INLINE_FUNCTION bool toBool(const Expr< ExprT > &x)
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)
SACADO_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
SACADO_INLINE_FUNCTION T if_then_else(const Cond cond, const T &a, const T &b)
Base template specification for Promote.