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()))
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) \
178 namespace 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 
376 FAD_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))
387 FAD_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))
398 FAD_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))
409 FAD_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) \
472 namespace 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 
500 FAD_RELOP_MACRO(==)
501 FAD_RELOP_MACRO(!=)
504 FAD_RELOP_MACRO(<=)
505 FAD_RELOP_MACRO(>=)
506 FAD_RELOP_MACRO(<<=)
507 FAD_RELOP_MACRO(>>=)
510 
511 #undef FAD_RELOP_MACRO
512 
513 namespace 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 -----------------------
528 namespace 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) \
545 namespace 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 
573 FAD_BOOL_MACRO(&&)
574 FAD_BOOL_MACRO(||)
575 
576 #undef FAD_BOOL_MACRO
577 
578 //-------------------------- I/O Operators -----------------------
579 
580 namespace 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
cbrt(expr.val())
expr expr SinOp
expr2 expr1 expr2 expr2 c *expr2 c *expr1 c *expr2 c *expr1 MaxOp
asinh(expr.val())
SACADO_INLINE_FUNCTION bool operator!(const Expr< ExprT > &expr)
#define FAD_UNARYOP_MACRO(OPNAME, OP, USING, VALUE, DX, FASTACCESSDX)
asin(expr.val())
cosh(expr.val())
expr expr dx(i)
abs(expr.val())
atanh(expr.val())
expr expr CoshOp
expr expr ATanhOp
expr expr TanhOp
expr expr SqrtOp
expr expr ASinhOp
atan(expr.val())
expr1 expr1 expr2 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 MultiplicationOp
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
SACADO_INLINE_FUNCTION bool toBool(const Expr< T > &xx)
acosh(expr.val())
acos(expr.val())
SimpleFad< ValueT > max(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
expr expr ASinOp
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)
#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