Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_LFad_LogicalSparseOps.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_LFAD_LOGICALSPARSEOPS_HPP
53#define SACADO_LFAD_LOGICALSPARSEOPS_HPP
54
55#include <cmath>
56#include <algorithm> // for std::min and std::max
57#include <ostream> // for std::ostream
58
59#define FAD_UNARYOP_MACRO(OPNAME,OP,VALUE) \
60namespace Sacado { \
61 namespace LFad { \
62 \
63 template <typename ExprT> \
64 class OP {}; \
65 \
66 template <typename ExprT> \
67 class Expr< OP<ExprT> > { \
68 public: \
69 \
70 typedef typename ExprT::value_type value_type; \
71 typedef typename ExprT::scalar_type scalar_type; \
72 typedef typename ExprT::base_expr_type base_expr_type; \
73 typedef typename ExprT::logical_type logical_type; \
74 \
75 Expr(const ExprT& expr_) : expr(expr_) {} \
76 \
77 int size() const { return expr.size(); } \
78 \
79 bool hasFastAccess() const { return expr.hasFastAccess(); } \
80 \
81 bool isPassive() const { return expr.isPassive();} \
82 \
83 value_type val() const { \
84 return VALUE; \
85 } \
86 \
87 logical_type dx(int i) const { \
88 return expr.dx(i); \
89 } \
90 \
91 logical_type fastAccessDx(int i) const { \
92 return expr.fastAccessDx(i); \
93 } \
94 \
95 protected: \
96 \
97 const ExprT& expr; \
98 }; \
99 \
100 template <typename T> \
101 inline Expr< OP< Expr<T> > > \
102 OPNAME (const Expr<T>& expr) \
103 { \
104 typedef OP< Expr<T> > expr_t; \
105 \
106 return Expr<expr_t>(expr); \
107 } \
108 } \
109}
110
112 UnaryPlusOp,
113 expr.val())
114FAD_UNARYOP_MACRO(operator-,
116 -expr.val())
119 std::exp(expr.val()))
122 std::log(expr.val()))
125 std::log10(expr.val()))
128 std::sqrt(expr.val()))
131 std::cos(expr.val()))
134 std::sin(expr.val()))
137 std::tan(expr.val()))
140 std::acos(expr.val()))
143 std::asin(expr.val()))
146 std::atan(expr.val()))
149 std::cosh(expr.val()))
152 std::sinh(expr.val()))
155 std::tanh(expr.val()))
158 acosh(expr.val()))
161 asinh(expr.val()))
164 atanh(expr.val()))
167 std::abs(expr.val()))
170 std::fabs(expr.val()))
173 std::cbrt(expr.val()))
174
175#undef FAD_UNARYOP_MACRO
176
177#define FAD_BINARYOP_MACRO(OPNAME,OP,VALUE,DX,FASTACCESSDX,VAL_CONST_DX_1,VAL_CONST_DX_2,CONST_DX_1,CONST_DX_2,CONST_FASTACCESSDX_1,CONST_FASTACCESSDX_2) \
178namespace Sacado { \
179 namespace LFad { \
180 \
181 template <typename ExprT1, typename ExprT2> \
182 class OP {}; \
183 \
184 template <typename T1, typename T2> \
185 class Expr< OP< Expr<T1>, Expr<T2> > > { \
186 \
187 public: \
188 \
189 typedef Expr<T1> ExprT1; \
190 typedef Expr<T2> ExprT2; \
191 typedef typename ExprT1::value_type value_type_1; \
192 typedef typename ExprT2::value_type value_type_2; \
193 typedef typename Sacado::Promote<value_type_1, \
194 value_type_2>::type value_type; \
195 \
196 typedef typename ExprT1::scalar_type scalar_type_1; \
197 typedef typename ExprT2::scalar_type scalar_type_2; \
198 typedef typename Sacado::Promote<scalar_type_1, \
199 scalar_type_2>::type scalar_type; \
200 \
201 typedef typename ExprT1::base_expr_type base_expr_type_1; \
202 typedef typename ExprT2::base_expr_type base_expr_type_2; \
203 typedef typename Sacado::Promote<base_expr_type_1, \
204 base_expr_type_2>::type base_expr_type; \
205 \
206 typedef typename ExprT1::logical_type logical_type; \
207 \
208 Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
209 expr1(expr1_), expr2(expr2_) {} \
210 \
211 int size() const { \
212 int sz1 = expr1.size(), sz2 = expr2.size(); \
213 return sz1 > sz2 ? sz1 : sz2; \
214 } \
215 \
216 bool hasFastAccess() const { \
217 return expr1.hasFastAccess() && expr2.hasFastAccess(); \
218 } \
219 \
220 bool isPassive() const { \
221 return expr1.isPassive() && expr2.isPassive(); \
222 } \
223 \
224 const value_type val() const { \
225 return VALUE; \
226 } \
227 \
228 const logical_type dx(int i) const { \
229 return DX; \
230 } \
231 \
232 const logical_type fastAccessDx(int i) const { \
233 return FASTACCESSDX; \
234 } \
235 \
236 protected: \
237 \
238 const ExprT1& expr1; \
239 const ExprT2& expr2; \
240 \
241 }; \
242 \
243 template <typename T1> \
244 class Expr< OP< Expr<T1>, typename Expr<T1>::value_type> > { \
245 \
246 public: \
247 \
248 typedef Expr<T1> ExprT1; \
249 typedef typename ExprT1::value_type value_type; \
250 typedef typename ExprT1::scalar_type scalar_type; \
251 typedef typename ExprT1::base_expr_type base_expr_type; \
252 typedef typename ExprT1::value_type ConstT; \
253 typedef typename ExprT1::logical_type logical_type; \
254 \
255 Expr(const ExprT1& expr1_, const ConstT& c_) : \
256 expr1(expr1_), c(c_) {} \
257 \
258 int size() const { \
259 return expr1.size(); \
260 } \
261 \
262 bool hasFastAccess() const { \
263 return expr1.hasFastAccess(); \
264 } \
265 \
266 bool isPassive() const { \
267 return expr1.isPassive(); \
268 } \
269 \
270 const value_type val() const { \
271 return VAL_CONST_DX_2; \
272 } \
273 \
274 const logical_type dx(int i) const { \
275 return CONST_DX_2; \
276 } \
277 \
278 const logical_type fastAccessDx(int i) const { \
279 return CONST_FASTACCESSDX_2; \
280 } \
281 \
282 protected: \
283 \
284 const ExprT1& expr1; \
285 const ConstT& c; \
286 }; \
287 \
288 template <typename T2> \
289 class Expr< OP< typename Expr<T2>::value_type, Expr<T2> > > { \
290 \
291 public: \
292 \
293 typedef Expr<T2> ExprT2; \
294 typedef typename ExprT2::value_type value_type; \
295 typedef typename ExprT2::scalar_type scalar_type; \
296 typedef typename ExprT2::base_expr_type base_expr_type; \
297 typedef typename ExprT2::value_type ConstT; \
298 typedef typename ExprT2::logical_type logical_type; \
299 \
300 Expr(const ConstT& c_, const ExprT2& expr2_) : \
301 c(c_), expr2(expr2_) {} \
302 \
303 int size() const { \
304 return expr2.size(); \
305 } \
306 \
307 bool hasFastAccess() const { \
308 return expr2.hasFastAccess(); \
309 } \
310 \
311 bool isPassive() const { \
312 return expr2.isPassive(); \
313 } \
314 \
315 const value_type val() const { \
316 return VAL_CONST_DX_1; \
317 } \
318 \
319 const logical_type dx(int i) const { \
320 return CONST_DX_1; \
321 } \
322 \
323 const logical_type fastAccessDx(int i) const { \
324 return CONST_FASTACCESSDX_1; \
325 } \
326 \
327 protected: \
328 \
329 const ConstT& c; \
330 const ExprT2& expr2; \
331 }; \
332 \
333 template <typename T1, typename T2> \
334 inline Expr< OP< Expr<T1>, Expr<T2> > > \
335 OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2) \
336 { \
337 typedef OP< Expr<T1>, Expr<T2> > expr_t; \
338 \
339 return Expr<expr_t>(expr1, expr2); \
340 } \
341 \
342 template <typename T> \
343 inline Expr< OP< Expr<T>, Expr<T> > > \
344 OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
345 { \
346 typedef OP< Expr<T>, Expr<T> > expr_t; \
347 \
348 return Expr<expr_t>(expr1, expr2); \
349 } \
350 \
351 template <typename T> \
352 inline Expr< OP< typename Expr<T>::value_type, Expr<T> > > \
353 OPNAME (const typename Expr<T>::value_type& c, \
354 const Expr<T>& expr) \
355 { \
356 typedef typename Expr<T>::value_type ConstT; \
357 typedef OP< ConstT, Expr<T> > expr_t; \
358 \
359 return Expr<expr_t>(c, expr); \
360 } \
361 \
362 template <typename T> \
363 inline Expr< OP< Expr<T>, typename Expr<T>::value_type > > \
364 OPNAME (const Expr<T>& expr, \
365 const typename Expr<T>::value_type& c) \
366 { \
367 typedef typename Expr<T>::value_type ConstT; \
368 typedef OP< Expr<T>, ConstT > expr_t; \
369 \
370 return Expr<expr_t>(expr, c); \
371 } \
372 } \
373}
374
375
376FAD_BINARYOP_MACRO(operator+,
378 expr1.val() + expr2.val(),
379 expr1.dx(i) || expr2.dx(i),
380 expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
381 c + expr2.val(),
382 expr1.val() + c,
383 expr2.dx(i),
384 expr1.dx(i),
385 expr2.fastAccessDx(i),
386 expr1.fastAccessDx(i))
387FAD_BINARYOP_MACRO(operator-,
389 expr1.val() - expr2.val(),
390 expr1.dx(i) || expr2.dx(i),
391 expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
392 c - expr2.val(),
393 expr1.val() - c,
394 expr2.dx(i),
395 expr1.dx(i),
396 expr2.fastAccessDx(i),
397 expr1.fastAccessDx(i))
398FAD_BINARYOP_MACRO(operator*,
400 expr1.val() * expr2.val(),
401 expr1.dx(i) || expr2.dx(i),
402 expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
403 c * expr2.val(),
404 expr1.val() * c,
405 expr2.dx(i),
406 expr1.dx(i),
407 expr2.fastAccessDx(i),
408 expr1.fastAccessDx(i))
409FAD_BINARYOP_MACRO(operator/,
411 expr1.val() / expr2.val(),
412 expr1.dx(i) || expr2.dx(i),
413 expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
414 c / expr2.val(),
415 expr1.val() / c,
416 expr2.dx(i),
417 expr1.dx(i),
418 expr2.fastAccessDx(i),
419 expr1.fastAccessDx(i))
422 std::atan2(expr1.val(), expr2.val()),
423 expr1.dx(i) || expr2.dx(i),
424 expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
425 std::atan2(c, expr2.val()),
426 std::atan2(expr1.val(), c),
427 expr2.dx(i),
428 expr1.dx(i),
429 expr2.fastAccessDx(i),
430 expr1.fastAccessDx(i))
433 std::pow(expr1.val(), expr2.val()),
434 expr1.dx(i) || expr2.dx(i),
435 expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
436 std::pow(c, expr2.val()),
437 std::pow(expr1.val(), c),
438 expr2.dx(i),
439 expr1.dx(i),
440 expr2.fastAccessDx(i),
441 expr1.fastAccessDx(i))
444 std::max(expr1.val(), expr2.val()),
445 expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
446 expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) :
447 expr2.fastAccessDx(i),
448 std::max(c, expr2.val()),
449 std::max(expr1.val(), c),
450 c >= expr2.val() ? logical_type(0) : expr2.dx(i),
451 expr1.val() >= c ? expr1.dx(i) : logical_type(0),
452 c >= expr2.val() ? logical_type(0) : expr2.fastAccessDx(i),
453 expr1.val() >= c ? expr1.fastAccessDx(i) : logical_type(0))
455 MinOp,
456 std::min(expr1.val(), expr2.val()),
457 expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
458 expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) :
459 expr2.fastAccessDx(i),
460 std::min(c, expr2.val()),
461 std::min(expr1.val(), c),
462 c <= expr2.val() ? logical_type(0) : expr2.dx(i),
463 expr1.val() <= c ? expr1.dx(i) : logical_type(0),
464 c <= expr2.val() ? logical_type(0) : expr2.fastAccessDx(i),
465 expr1.val() <= c ? expr1.fastAccessDx(i) : logical_type(0))
466
467#undef FAD_BINARYOP_MACRO
468
469//-------------------------- Relational Operators -----------------------
470
471#define FAD_RELOP_MACRO(OP) \
472namespace Sacado { \
473 namespace LFad { \
474 template <typename ExprT1, typename ExprT2> \
475 inline bool \
476 operator OP (const Expr<ExprT1>& expr1, \
477 const Expr<ExprT2>& expr2) \
478 { \
479 return expr1.val() OP expr2.val(); \
480 } \
481 \
482 template <typename ExprT2> \
483 inline bool \
484 operator OP (const typename Expr<ExprT2>::value_type& a, \
485 const Expr<ExprT2>& expr2) \
486 { \
487 return a OP expr2.val(); \
488 } \
489 \
490 template <typename ExprT1> \
491 inline bool \
492 operator OP (const Expr<ExprT1>& expr1, \
493 const typename Expr<ExprT1>::value_type& b) \
494 { \
495 return expr1.val() OP b; \
496 } \
497 } \
498}
499
510
511#undef FAD_RELOP_MACRO
512
513namespace Sacado {
514
515 namespace LFad {
516
517 template <typename ExprT>
518 inline bool operator ! (const Expr<ExprT>& expr)
519 {
520 return ! expr.val();
521 }
522
523 } // namespace LFad
524
525} // namespace Sacado
526
527//-------------------------- Boolean Operators -----------------------
528namespace Sacado {
529
530 namespace LFad {
531
532 template <typename ExprT>
533 bool toBool(const Expr<ExprT>& x) {
534 bool is_zero = (x.val() == 0.0);
535 for (int i=0; i<x.size(); i++)
536 is_zero = is_zero && (x.dx(i) == 0.0);
537 return !is_zero;
538 }
539
540 } // namespace Fad
541
542} // namespace Sacado
543
544#define FAD_BOOL_MACRO(OP) \
545namespace Sacado { \
546 namespace LFad { \
547 template <typename ExprT1, typename ExprT2> \
548 inline bool \
549 operator OP (const Expr<ExprT1>& expr1, \
550 const Expr<ExprT2>& expr2) \
551 { \
552 return toBool(expr1) OP toBool(expr2); \
553 } \
554 \
555 template <typename ExprT2> \
556 inline bool \
557 operator OP (const typename Expr<ExprT2>::value_type& a, \
558 const Expr<ExprT2>& expr2) \
559 { \
560 return a OP toBool(expr2); \
561 } \
562 \
563 template <typename ExprT1> \
564 inline bool \
565 operator OP (const Expr<ExprT1>& expr1, \
566 const typename Expr<ExprT1>::value_type& b) \
567 { \
568 return toBool(expr1) OP b; \
569 } \
570 } \
571}
572
575
576#undef FAD_BOOL_MACRO
577
578//-------------------------- I/O Operators -----------------------
579
580namespace Sacado {
581
582 namespace LFad {
583
584 template <typename ExprT>
585 std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
586 os << x.val() << " [";
587
588 for (int i=0; i< x.size(); i++) {
589 os << " " << x.dx(i);
590 }
591
592 os << " ]";
593 return os;
594 }
595
596 } // namespace LFad
597
598} // namespace Sacado
599
600
601#endif // SACADO_LFAD_LOGICALSPARSEOPS_HPP
fabs(expr.val())
asinh(expr.val())
log(expr.val())
tan(expr.val())
abs(expr.val())
expr expr SinOp
expr expr SinhOp
cos(expr.val())
expr expr SqrtOp
expr expr dx(i)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
expr expr ACosOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
expr expr ATanOp
cosh(expr.val())
expr expr ACoshOp
acos(expr.val())
sin(expr.val())
sinh(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 DivisionOp
expr expr ASinOp
cbrt(expr.val())
log10(expr.val())
expr expr TanhOp
expr expr TanOp
expr expr CoshOp
#define FAD_BINARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, CDX1, CDX2, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
#define FAD_RELOP_MACRO(OP)
expr expr ASinhOp
exp(expr.val())
atan(expr.val())
expr expr Log10Op
#define FAD_BOOL_MACRO(OP)
acosh(expr.val())
expr expr expr ExpOp
expr expr AbsOp
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
atan2(expr1.val(), expr2.val())
sqrt(expr.val())
expr expr ATanhOp
tanh(expr.val())
atanh(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 Atan2Op
expr expr CosOp
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
asin(expr.val())
#define FAD_UNARYOP_MACRO(OPNAME, OP, VALUE)
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
#define FAD_RELOP_MACRO(OP)
#define FAD_BINARYOP_MACRO(OPNAME, OP, VALUE, DX, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
#define FAD_BOOL_MACRO(OP)
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 PowerOp
adouble max(const adouble &a, const adouble &b)
adouble min(const adouble &a, const adouble &b)
Wrapper for a generic expression template.
SACADO_INLINE_FUNCTION bool toBool(const Expr< T > &xx)
Namespace for logical forward-mode AD classes.
std::ostream & operator<<(std::ostream &os, const ParameterLibraryBase< FamilyType, EntryType > &pl)
Taylor< T > max(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
PowExprType< Expr< T1 >, Expr< T2 > >::expr_type pow(const Expr< T1 > &expr1, const Expr< T2 > &expr2)