Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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) \
60 namespace 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 
111 FAD_UNARYOP_MACRO(operator+,
112  UnaryPlusOp,
113  expr.val())
114 FAD_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()))
171 #ifdef HAVE_SACADO_CXX11
173  CbrtOp,
174  std::cbrt(expr.val()))
175 #endif
176 
177 #undef FAD_UNARYOP_MACRO
178 
179 #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) \
180 namespace Sacado { \
181  namespace LFad { \
182  \
183  template <typename ExprT1, typename ExprT2> \
184  class OP {}; \
185  \
186  template <typename T1, typename T2> \
187  class Expr< OP< Expr<T1>, Expr<T2> > > { \
188  \
189  public: \
190  \
191  typedef Expr<T1> ExprT1; \
192  typedef Expr<T2> ExprT2; \
193  typedef typename ExprT1::value_type value_type_1; \
194  typedef typename ExprT2::value_type value_type_2; \
195  typedef typename Sacado::Promote<value_type_1, \
196  value_type_2>::type value_type; \
197  \
198  typedef typename ExprT1::scalar_type scalar_type_1; \
199  typedef typename ExprT2::scalar_type scalar_type_2; \
200  typedef typename Sacado::Promote<scalar_type_1, \
201  scalar_type_2>::type scalar_type; \
202  \
203  typedef typename ExprT1::base_expr_type base_expr_type_1; \
204  typedef typename ExprT2::base_expr_type base_expr_type_2; \
205  typedef typename Sacado::Promote<base_expr_type_1, \
206  base_expr_type_2>::type base_expr_type; \
207  \
208  typedef typename ExprT1::logical_type logical_type; \
209  \
210  Expr(const ExprT1& expr1_, const ExprT2& expr2_) : \
211  expr1(expr1_), expr2(expr2_) {} \
212  \
213  int size() const { \
214  int sz1 = expr1.size(), sz2 = expr2.size(); \
215  return sz1 > sz2 ? sz1 : sz2; \
216  } \
217  \
218  bool hasFastAccess() const { \
219  return expr1.hasFastAccess() && expr2.hasFastAccess(); \
220  } \
221  \
222  bool isPassive() const { \
223  return expr1.isPassive() && expr2.isPassive(); \
224  } \
225  \
226  const value_type val() const { \
227  return VALUE; \
228  } \
229  \
230  const logical_type dx(int i) const { \
231  return DX; \
232  } \
233  \
234  const logical_type fastAccessDx(int i) const { \
235  return FASTACCESSDX; \
236  } \
237  \
238  protected: \
239  \
240  const ExprT1& expr1; \
241  const ExprT2& expr2; \
242  \
243  }; \
244  \
245  template <typename T1> \
246  class Expr< OP< Expr<T1>, typename Expr<T1>::value_type> > { \
247  \
248  public: \
249  \
250  typedef Expr<T1> ExprT1; \
251  typedef typename ExprT1::value_type value_type; \
252  typedef typename ExprT1::scalar_type scalar_type; \
253  typedef typename ExprT1::base_expr_type base_expr_type; \
254  typedef typename ExprT1::value_type ConstT; \
255  typedef typename ExprT1::logical_type logical_type; \
256  \
257  Expr(const ExprT1& expr1_, const ConstT& c_) : \
258  expr1(expr1_), c(c_) {} \
259  \
260  int size() const { \
261  return expr1.size(); \
262  } \
263  \
264  bool hasFastAccess() const { \
265  return expr1.hasFastAccess(); \
266  } \
267  \
268  bool isPassive() const { \
269  return expr1.isPassive(); \
270  } \
271  \
272  const value_type val() const { \
273  return VAL_CONST_DX_2; \
274  } \
275  \
276  const logical_type dx(int i) const { \
277  return CONST_DX_2; \
278  } \
279  \
280  const logical_type fastAccessDx(int i) const { \
281  return CONST_FASTACCESSDX_2; \
282  } \
283  \
284  protected: \
285  \
286  const ExprT1& expr1; \
287  const ConstT& c; \
288  }; \
289  \
290  template <typename T2> \
291  class Expr< OP< typename Expr<T2>::value_type, Expr<T2> > > { \
292  \
293  public: \
294  \
295  typedef Expr<T2> ExprT2; \
296  typedef typename ExprT2::value_type value_type; \
297  typedef typename ExprT2::scalar_type scalar_type; \
298  typedef typename ExprT2::base_expr_type base_expr_type; \
299  typedef typename ExprT2::value_type ConstT; \
300  typedef typename ExprT2::logical_type logical_type; \
301  \
302  Expr(const ConstT& c_, const ExprT2& expr2_) : \
303  c(c_), expr2(expr2_) {} \
304  \
305  int size() const { \
306  return expr2.size(); \
307  } \
308  \
309  bool hasFastAccess() const { \
310  return expr2.hasFastAccess(); \
311  } \
312  \
313  bool isPassive() const { \
314  return expr2.isPassive(); \
315  } \
316  \
317  const value_type val() const { \
318  return VAL_CONST_DX_1; \
319  } \
320  \
321  const logical_type dx(int i) const { \
322  return CONST_DX_1; \
323  } \
324  \
325  const logical_type fastAccessDx(int i) const { \
326  return CONST_FASTACCESSDX_1; \
327  } \
328  \
329  protected: \
330  \
331  const ConstT& c; \
332  const ExprT2& expr2; \
333  }; \
334  \
335  template <typename T1, typename T2> \
336  inline Expr< OP< Expr<T1>, Expr<T2> > > \
337  OPNAME (const Expr<T1>& expr1, const Expr<T2>& expr2) \
338  { \
339  typedef OP< Expr<T1>, Expr<T2> > expr_t; \
340  \
341  return Expr<expr_t>(expr1, expr2); \
342  } \
343  \
344  template <typename T> \
345  inline Expr< OP< Expr<T>, Expr<T> > > \
346  OPNAME (const Expr<T>& expr1, const Expr<T>& expr2) \
347  { \
348  typedef OP< Expr<T>, Expr<T> > expr_t; \
349  \
350  return Expr<expr_t>(expr1, expr2); \
351  } \
352  \
353  template <typename T> \
354  inline Expr< OP< typename Expr<T>::value_type, Expr<T> > > \
355  OPNAME (const typename Expr<T>::value_type& c, \
356  const Expr<T>& expr) \
357  { \
358  typedef typename Expr<T>::value_type ConstT; \
359  typedef OP< ConstT, Expr<T> > expr_t; \
360  \
361  return Expr<expr_t>(c, expr); \
362  } \
363  \
364  template <typename T> \
365  inline Expr< OP< Expr<T>, typename Expr<T>::value_type > > \
366  OPNAME (const Expr<T>& expr, \
367  const typename Expr<T>::value_type& c) \
368  { \
369  typedef typename Expr<T>::value_type ConstT; \
370  typedef OP< Expr<T>, ConstT > expr_t; \
371  \
372  return Expr<expr_t>(expr, c); \
373  } \
374  } \
375 }
376 
377 
378 FAD_BINARYOP_MACRO(operator+,
380  expr1.val() + expr2.val(),
381  expr1.dx(i) || expr2.dx(i),
382  expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
383  c + expr2.val(),
384  expr1.val() + c,
385  expr2.dx(i),
386  expr1.dx(i),
387  expr2.fastAccessDx(i),
388  expr1.fastAccessDx(i))
389 FAD_BINARYOP_MACRO(operator-,
391  expr1.val() - expr2.val(),
392  expr1.dx(i) || expr2.dx(i),
393  expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
394  c - expr2.val(),
395  expr1.val() - c,
396  expr2.dx(i),
397  expr1.dx(i),
398  expr2.fastAccessDx(i),
399  expr1.fastAccessDx(i))
400 FAD_BINARYOP_MACRO(operator*,
402  expr1.val() * expr2.val(),
403  expr1.dx(i) || expr2.dx(i),
404  expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
405  c * expr2.val(),
406  expr1.val() * c,
407  expr2.dx(i),
408  expr1.dx(i),
409  expr2.fastAccessDx(i),
410  expr1.fastAccessDx(i))
411 FAD_BINARYOP_MACRO(operator/,
413  expr1.val() / expr2.val(),
414  expr1.dx(i) || expr2.dx(i),
415  expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
416  c / expr2.val(),
417  expr1.val() / c,
418  expr2.dx(i),
419  expr1.dx(i),
420  expr2.fastAccessDx(i),
421  expr1.fastAccessDx(i))
424  std::atan2(expr1.val(), expr2.val()),
425  expr1.dx(i) || expr2.dx(i),
426  expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
427  std::atan2(c, expr2.val()),
428  std::atan2(expr1.val(), c),
429  expr2.dx(i),
430  expr1.dx(i),
431  expr2.fastAccessDx(i),
432  expr1.fastAccessDx(i))
435  std::pow(expr1.val(), expr2.val()),
436  expr1.dx(i) || expr2.dx(i),
437  expr1.fastAccessDx(i) || expr2.fastAccessDx(i),
438  std::pow(c, expr2.val()),
439  std::pow(expr1.val(), c),
440  expr2.dx(i),
441  expr1.dx(i),
442  expr2.fastAccessDx(i),
443  expr1.fastAccessDx(i))
446  std::max(expr1.val(), expr2.val()),
447  expr1.val() >= expr2.val() ? expr1.dx(i) : expr2.dx(i),
448  expr1.val() >= expr2.val() ? expr1.fastAccessDx(i) :
449  expr2.fastAccessDx(i),
450  std::max(c, expr2.val()),
451  std::max(expr1.val(), c),
452  c >= expr2.val() ? logical_type(0) : expr2.dx(i),
453  expr1.val() >= c ? expr1.dx(i) : logical_type(0),
454  c >= expr2.val() ? logical_type(0) : expr2.fastAccessDx(i),
455  expr1.val() >= c ? expr1.fastAccessDx(i) : logical_type(0))
457  MinOp,
458  std::min(expr1.val(), expr2.val()),
459  expr1.val() <= expr2.val() ? expr1.dx(i) : expr2.dx(i),
460  expr1.val() <= expr2.val() ? expr1.fastAccessDx(i) :
461  expr2.fastAccessDx(i),
462  std::min(c, expr2.val()),
463  std::min(expr1.val(), c),
464  c <= expr2.val() ? logical_type(0) : expr2.dx(i),
465  expr1.val() <= c ? expr1.dx(i) : logical_type(0),
466  c <= expr2.val() ? logical_type(0) : expr2.fastAccessDx(i),
467  expr1.val() <= c ? expr1.fastAccessDx(i) : logical_type(0))
468 
469 #undef FAD_BINARYOP_MACRO
470 
471 //-------------------------- Relational Operators -----------------------
472 
473 #define FAD_RELOP_MACRO(OP) \
474 namespace Sacado { \
475  namespace LFad { \
476  template <typename ExprT1, typename ExprT2> \
477  inline bool \
478  operator OP (const Expr<ExprT1>& expr1, \
479  const Expr<ExprT2>& expr2) \
480  { \
481  return expr1.val() OP expr2.val(); \
482  } \
483  \
484  template <typename ExprT2> \
485  inline bool \
486  operator OP (const typename Expr<ExprT2>::value_type& a, \
487  const Expr<ExprT2>& expr2) \
488  { \
489  return a OP expr2.val(); \
490  } \
491  \
492  template <typename ExprT1> \
493  inline bool \
494  operator OP (const Expr<ExprT1>& expr1, \
495  const typename Expr<ExprT1>::value_type& b) \
496  { \
497  return expr1.val() OP b; \
498  } \
499  } \
500 }
501 
502 FAD_RELOP_MACRO(==)
503 FAD_RELOP_MACRO(!=)
506 FAD_RELOP_MACRO(<=)
507 FAD_RELOP_MACRO(>=)
508 FAD_RELOP_MACRO(<<=)
509 FAD_RELOP_MACRO(>>=)
512 
513 #undef FAD_RELOP_MACRO
514 
515 namespace Sacado {
516 
517  namespace LFad {
518 
519  template <typename ExprT>
520  inline bool operator ! (const Expr<ExprT>& expr)
521  {
522  return ! expr.val();
523  }
524 
525  } // namespace LFad
526 
527 } // namespace Sacado
528 
529 //-------------------------- Boolean Operators -----------------------
530 namespace Sacado {
531 
532  namespace LFad {
533 
534  template <typename ExprT>
535  bool toBool(const Expr<ExprT>& x) {
536  bool is_zero = (x.val() == 0.0);
537  for (int i=0; i<x.size(); i++)
538  is_zero = is_zero && (x.dx(i) == 0.0);
539  return !is_zero;
540  }
541 
542  } // namespace Fad
543 
544 } // namespace Sacado
545 
546 #define FAD_BOOL_MACRO(OP) \
547 namespace Sacado { \
548  namespace LFad { \
549  template <typename ExprT1, typename ExprT2> \
550  inline bool \
551  operator OP (const Expr<ExprT1>& expr1, \
552  const Expr<ExprT2>& expr2) \
553  { \
554  return toBool(expr1) OP toBool(expr2); \
555  } \
556  \
557  template <typename ExprT2> \
558  inline bool \
559  operator OP (const typename Expr<ExprT2>::value_type& a, \
560  const Expr<ExprT2>& expr2) \
561  { \
562  return a OP toBool(expr2); \
563  } \
564  \
565  template <typename ExprT1> \
566  inline bool \
567  operator OP (const Expr<ExprT1>& expr1, \
568  const typename Expr<ExprT1>::value_type& b) \
569  { \
570  return toBool(expr1) OP b; \
571  } \
572  } \
573 }
574 
575 FAD_BOOL_MACRO(&&)
576 FAD_BOOL_MACRO(||)
577 
578 #undef FAD_BOOL_MACRO
579 
580 //-------------------------- I/O Operators -----------------------
581 
582 namespace Sacado {
583 
584  namespace LFad {
585 
586  template <typename ExprT>
587  std::ostream& operator << (std::ostream& os, const Expr<ExprT>& x) {
588  os << x.val() << " [";
589 
590  for (int i=0; i< x.size(); i++) {
591  os << " " << x.dx(i);
592  }
593 
594  os << " ]";
595  return os;
596  }
597 
598  } // namespace LFad
599 
600 } // namespace Sacado
601 
602 
603 #endif // SACADO_LFAD_LOGICALSPARSEOPS_HPP
cbrt(expr.val())
expr expr SinOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
asinh(expr.val())
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
asin(expr.val())
cosh(expr.val())
expr expr dx(i)
abs(expr.val())
KOKKOS_INLINE_FUNCTION bool toBool(const Expr< T > &xx)
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MinOp
atanh(expr.val())
expr expr CoshOp
expr expr ATanhOp
expr expr TanhOp
expr expr SqrtOp
expr expr ASinhOp
atan(expr.val())
KOKKOS_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
KOKKOS_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)
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 Atan2Op
expr val()
tanh(expr.val())
expr expr CosOp
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_BINARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, CDX1, CDX2, FASTACCESSDX, VAL_CONST_DX_1, VAL_CONST_DX_2, CONST_DX_1, CONST_DX_2, CONST_FASTACCESSDX_1, CONST_FASTACCESSDX_2)
expr expr ATanOp
SimpleFad< ValueT > min(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ACosOp
sqrt(expr.val())
sinh(expr.val())
tan(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 expr1 expr1 c *expr2 expr1 c *expr2 expr1 c *expr2 expr1 DivisionOp
#define FAD_RELOP_MACRO(OP)
atan2(expr1.val(), expr2.val())
sin(expr.val())
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
log(expr.val())
expr expr ACoshOp
expr expr Log10Op
expr expr SinhOp
acosh(expr.val())
acos(expr.val())
SimpleFad< ValueT > max(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ASinOp
#define FAD_BOOL_MACRO(OP)
exp(expr.val())
expr expr expr ExpOp
fabs(expr.val())
expr expr AbsOp
expr expr TanOp
log10(expr.val())
cos(expr.val())
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 PowerOp