Sacado Package Browser (Single Doxygen Collection)  Version of the Day
Sacado_ELRFad_SFad.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Sacado Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25 // (etphipp@sandia.gov).
26 //
27 // ***********************************************************************
28 //
29 // The forward-mode AD classes in Sacado are a derivative work of the
30 // expression template classes in the Fad package by Nicolas Di Cesare.
31 // The following banner is included in the original Fad source code:
32 //
33 // ************ DO NOT REMOVE THIS BANNER ****************
34 //
35 // Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
36 // http://www.ann.jussieu.fr/~dicesare
37 //
38 // CEMRACS 98 : C++ courses,
39 // templates : new C++ techniques
40 // for scientific computing
41 //
42 //********************************************************
43 //
44 // A short implementation ( not all operators and
45 // functions are overloaded ) of 1st order Automatic
46 // Differentiation in forward mode (FAD) using
47 // EXPRESSION TEMPLATES.
48 //
49 //********************************************************
50 // @HEADER
51 
52 #ifndef SACADO_ELRFAD_SFAD_HPP
53 #define SACADO_ELRFAD_SFAD_HPP
54 
58 #include "Sacado_mpl_range_c.hpp"
59 #include "Sacado_mpl_for_each.hpp"
60 
61 namespace Sacado {
62 
63  namespace ELRFad {
64 
66  template <typename T, int Num>
67  struct SFadExprTag {};
68 
69  // Forward declaration
70  template <typename T, int Num> class SFad;
71 
79  template <typename T, int Num>
80  class Expr< SFadExprTag<T,Num> > {
81 
82  public:
83 
85  typedef typename RemoveConst<T>::type value_type;
86 
89 
92 
94  static const int num_args = 1;
95 
97  static const bool is_linear = true;
98 
103 
106  Expr() : val_( T(0.)) { ss_array<T>::zero(dx_, Num); }
107 
109 
112  template <typename S>
115  val_(x) {
116  ss_array<T>::zero(dx_, Num);
117  }
118 
120 
124  Expr(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) : val_(x) {
125 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
126  if (sz != Num)
127  throw "SELRFad::SFad() Error: Supplied derivative dimension does not match compile time length.";
128 #endif
129 
130  if (zero_out == InitDerivArray)
131  ss_array<T>::zero(dx_, Num);
132  }
133 
135 
141  Expr(const int sz, const int i, const T & x) :
142  val_(x) {
143 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
144  if (sz != Num)
145  throw "SELRFad::SFad() Error: Supplied derivative dimension does not match compile time length.";
146  if (i >= Num)
147  throw "SELRFad::SFad() Error: Invalid derivative index.";
148 #endif
149 
150  ss_array<T>::zero(dx_, Num);
151  dx_[i]=1.;
152  }
153 
156  Expr(const Expr& x) : val_(x.val_) {
157  //ss_array<T>::copy(x.dx_, dx_, Num);
158  for (int i=0; i<Num; i++)
159  dx_[i] = x.dx_[i];
160  }
161 
163  template <typename S>
166 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
167  if (x.size() != Num)
168  throw "SELRFad::SFad() Error: Attempt to assign with incompatible sizes";
169 #endif
170 
171  // Compute partials
172  LocalAccumOp< Expr<S> > op(x);
173 
174  // Compute each tangent direction
175  for(int i=0; i<Num; ++i) {
176  op.t = T(0.);
177  op.i = i;
178 
179  // Automatically unrolled loop that computes
180  // for (int j=0; j<N; j++)
181  // op.t += op.partials[j] * x.getTangent<j>(i);
183 
184  dx_[i] = op.t;
185 
186  }
187 
188  this->val() = x.val();
189  }
190 
193  ~Expr() {}
194 
196 
203  void diff(const int ith, const int n) {
204 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
205  if (n != Num)
206  throw "SELRFad::diff() Error: Supplied derivative dimension does not match compile time length.";
207 #endif
208 
209  ss_array<T>::zero(dx_, Num);
210  dx_[ith] = T(1.);
211  }
212 
214 
219  void resize(int sz) {
220 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
221  if (sz != Num)
222  throw "SELRFad::resize() Error: Cannot resize fixed derivative array dimension";
223 #endif
224  }
225 
227 
232  void expand(int sz) { resize(sz); }
233 
236  void zero() { ss_array<T>::zero(dx_, Num); }
237 
240  void setUpdateValue(bool update_val) { }
241 
244  bool updateValue() const { return true; }
245 
247  template <typename S>
249  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
250  typedef IsEqual<value_type> IE;
251  if (x.size() != this->size()) return false;
252  bool eq = IE::eval(x.val(), this->val());
253  for (int i=0; i<this->size(); i++)
254  eq = eq && IE::eval(x.dx(i), this->dx(i));
255  return eq;
256  }
257 
259 
264 
267  const T& val() const { return val_;}
268 
271  T& val() { return val_;}
272 
274 
279 
282  int size() const { return Num;}
283 
289  int availableSize() const { return Num; }
290 
293  bool hasFastAccess() const { return true; }
294 
297  bool isPassive() const { return false; }
298 
301  void setIsConstant(bool is_const) {}
302 
305  const T* dx() const { return &(dx_[0]);}
306 
309  const T& dx(int i) const { return dx_[i]; }
310 
313  T& fastAccessDx(int i) { return dx_[i];}
314 
317  const T& fastAccessDx(int i) const { return dx_[i];}
318 
321  void computePartials(const T& bar, value_type partials[]) const {
322  partials[0] = bar;
323  }
324 
327  void getTangents(int i, value_type dots[]) const {
328  dots[0] = this->dx_[i];
329  }
330 
332  template <int Arg>
334  bool isActive() const { return true; }
335 
337  template <int Arg>
339  const T& getTangent(int i) const { return this->dx_[i]; }
340 
343  const value_type* getDx(int j) const { return this->dx(); }
344 
346 
351 
353  template <typename S>
355  SACADO_ENABLE_VALUE_FUNC(Expr&) operator=(const S& v) {
356  val_ = v;
357  ss_array<T>::zero(dx_, Num);
358  return *this;
359  }
360 
363  Expr& operator=(const Expr& x) {
364  if (this != &x) {
365  // Copy value
366  val_ = x.val_;
367 
368  // Copy dx_
369  //ss_array<T>::copy(x.dx_, dx_, Num);
370  for (int i=0; i<Num; i++)
371  dx_[i] = x.dx_[i];
372  }
373  return *this;
374  }
375 
377  template <typename S>
379  SACADO_ENABLE_EXPR_FUNC(Expr&) operator=(const Expr<S>& x) {
380 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
381  if (x.size() != Num)
382  throw "SELRFad::operator=() Error: Attempt to assign with incompatible sizes";
383 #endif
384 
385  // Compute partials
386  LocalAccumOp< Expr<S> > op(x);
387 
388  // Compute each tangent direction
389  for(int i=0; i<Num; ++i) {
390  op.t = T(0.);
391  op.i = i;
392 
393  // Automatically unrolled loop that computes
394  // for (int j=0; j<N; j++)
395  // op.t += op.partials[j] * x.getTangent<j>(i);
397 
398  dx_[i] = op.t;
399  }
400 
401  // Compute value
402  val_ = x.val();
403 
404  return *this;
405  }
406 
408 
413 
415  template <typename S>
417  SACADO_ENABLE_VALUE_FUNC(Expr&) operator += (const S& v) {
418  this->val() += v;
419  return *this;
420  }
421 
423  template <typename S>
425  SACADO_ENABLE_VALUE_FUNC(Expr&) operator -= (const S& v) {
426  this->val() -= v;
427  return *this;
428  }
429 
431  template <typename S>
433  SACADO_ENABLE_VALUE_FUNC(Expr&) operator *= (const S& v) {
434  this->val() *= v;
435  for (int i=0; i<Num; ++i)
436  dx_[i] *= v;
437  return *this;
438  }
439 
441  template <typename S>
443  SACADO_ENABLE_VALUE_FUNC(Expr&) operator /= (const S& v) {
444  this->val() /= v;
445  for (int i=0; i<Num; ++i)
446  dx_[i] /= v;
447  return *this;
448  }
449 
451  template <typename S>
453  SACADO_ENABLE_EXPR_FUNC(Expr&) operator += (const Expr<S>& x) {
454 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
455  if (x.size() != Num)
456  throw "SELRFad::operator+=() Error: Attempt to assign with incompatible sizes";
457 #endif
458 
459  // Compute partials
460  LocalAccumOp< Expr<S> > op(x);
461 
462  // Compute each tangent direction
463  for(int i=0; i<Num; ++i) {
464  op.t = T(0.);
465  op.i = i;
466 
467  // Automatically unrolled loop that computes
468  // for (int j=0; j<N; j++)
469  // op.t += op.partials[j] * x.getTangent<j>(i);
471 
472  dx_[i] += op.t;
473  }
474 
475  // Compute value
476  val_ += x.val();
477 
478  return *this;
479  }
480 
482  template <typename S>
484  SACADO_ENABLE_EXPR_FUNC(Expr&) operator -= (const Expr<S>& x) {
485 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
486  if (x.size() != Num)
487  throw "SELRFad::operator-=() Error: Attempt to assign with incompatible sizes";
488 #endif
489 
490  // Compute partials
491  LocalAccumOp< Expr<S> > op(x);
492 
493  // Compute each tangent direction
494  for(int i=0; i<Num; ++i) {
495  op.t = T(0.);
496  op.i = i;
497 
498  // Automatically unrolled loop that computes
499  // for (int j=0; j<N; j++)
500  // op.t += op.partials[j] * x.getTangent<j>(i);
502 
503  dx_[i] -= op.t;
504  }
505 
506  // Compute value
507  val_ -= x.val();
508 
509  return *this;
510  }
511 
513  template <typename S>
515  SACADO_ENABLE_EXPR_FUNC(Expr&) operator *= (const Expr<S>& x) {
516  T xval = x.val();
517 
518 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
519  if (x.size() != Num)
520  throw "SELRFad::operator*=() Error: Attempt to assign with incompatible sizes";
521 #endif
522 
523  // Compute partials
524  LocalAccumOp< Expr<S> > op(x);
525 
526  // Compute each tangent direction
527  for(int i=0; i<Num; ++i) {
528  op.t = T(0.);
529  op.i = i;
530 
531  // Automatically unrolled loop that computes
532  // for (int j=0; j<N; j++)
533  // op.t += op.partials[j] * x.getTangent<j>(i);
535 
536  dx_[i] = val_ * op.t + dx_[i] * xval;
537  }
538 
539  // Compute value
540  val_ *= xval;
541 
542  return *this;
543  }
544 
546  template <typename S>
548  SACADO_ENABLE_EXPR_FUNC(Expr&) operator /= (const Expr<S>& x) {
549  T xval = x.val();
550 
551 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
552  if (x.size() != Num)
553  throw "SELRFad::operator/=() Error: Attempt to assign with incompatible sizes";
554 #endif
555 
556  // Compute partials
557  LocalAccumOp< Expr<S> > op(x);
558 
559  T xval2 = xval*xval;
560 
561  // Compute each tangent direction
562  for(int i=0; i<Num; ++i) {
563  op.t = T(0.);
564  op.i = i;
565 
566  // Automatically unrolled loop that computes
567  // for (int j=0; j<N; j++)
568  // op.t += op.partials[j] * x.getTangent<j>(i);
570 
571  dx_[i] = (dx_[i] * xval - val_ * op.t) / xval2;
572  }
573 
574  // Compute value
575  val_ /= xval;
576 
577  return *this;
578  }
579 
581 
582  protected:
583 
585  T dx_[Num];
586 
589 
590  // Functor for mpl::for_each to compute the local accumulation
591  // of a tangent derivative
592  //
593  // We use getTangent<>() to get dx components from expression
594  // arguments instead of getting the argument directly or extracting
595  // the dx array due to striding in ViewFad (or could use striding
596  // directly here if we need to get dx array).
597  template <typename ExprT>
598  struct LocalAccumOp {
599  typedef typename ExprT::value_type value_type;
600  static const int N = ExprT::num_args;
601  const ExprT& x;
602  mutable value_type t;
603  value_type partials[N];
604  int i;
606  LocalAccumOp(const ExprT& x_) : x(x_) {
607  x.computePartials(value_type(1.), partials);
608  }
610  void getTangents(int i_) { i = i_; }
611  template <typename ArgT>
613  void operator () (ArgT arg) const {
614  const int Arg = ArgT::value;
615  if (x.template isActive<Arg>())
616  t += partials[Arg] * x.template getTangent<Arg>(i);
617  }
618  };
619 
620  }; // class Expr<SFadExprTag>
621 
622  } // namespace ELRFad
623 
624 } // namespace Sacado
625 
626 #define FAD_NS ELRFad
627 #include "Sacado_Fad_SFad_tmpl.hpp"
628 #undef FAD_NS
629 
630 #include "Sacado_ELRFad_ViewFad.hpp"
631 #include "Sacado_ELRFad_Ops.hpp"
632 
633 #endif // SACADO_ELRFAD_SFAD_HPP
KOKKOS_INLINE_FUNCTION Expr(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
expr val()
KOKKOS_INLINE_FUNCTION const T & fastAccessDx(int i) const
Returns derivative component i without bounds checking.
SFad< value_type, Num > base_expr_type
Typename of base-expressions.
KOKKOS_INLINE_FUNCTION const value_type * getDx(int j) const
Get dx array.
KOKKOS_INLINE_FUNCTION T & val()
Returns value.
#define SACADO_ENABLE_VALUE_CTOR_DECL
RemoveConst< T >::type value_type
Typename of values.
KOKKOS_INLINE_FUNCTION const T & val() const
Returns value.
KOKKOS_INLINE_FUNCTION void expand(int sz)
Expand derivative array to size sz.
expr bar
#define SACADO_ENABLE_EXPR_CTOR_DECL
static KOKKOS_INLINE_FUNCTION void zero(T *dest, int sz)
Zero out array dest of length sz.
KOKKOS_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
KOKKOS_INLINE_FUNCTION int size() const
Returns number of derivative components.
KOKKOS_INLINE_FUNCTION SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr< S > &x) const
Returns whether two Fad objects have the same values.
KOKKOS_INLINE_FUNCTION Expr(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
#define KOKKOS_INLINE_FUNCTION
#define T
Definition: Sacado_rad.hpp:573
KOKKOS_INLINE_FUNCTION T & fastAccessDx(int i)
Returns derivative component i without bounds checking.
#define SACADO_ENABLE_VALUE_FUNC(RETURN_TYPE)
KOKKOS_INLINE_FUNCTION Expr(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
KOKKOS_INLINE_FUNCTION void resize(int sz)
Resize derivative array to length sz.
Base template specification for testing equivalence.
KOKKOS_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
KOKKOS_INLINE_FUNCTION void computePartials(const T &bar, value_type partials[]) const
Return partials w.r.t. arguments.
KOKKOS_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
KOKKOS_INLINE_FUNCTION const T & getTangent(int i) const
Return tangent component i of argument Arg.
KOKKOS_INLINE_FUNCTION const T & dx(int i) const
Returns derivative component i with bounds checking.
KOKKOS_INLINE_FUNCTION ~Expr()
Destructor.
KOKKOS_INLINE_FUNCTION void diff(const int ith, const int n)
Set Fad object as the ith independent variable.
KOKKOS_INLINE_FUNCTION void getTangents(int i, value_type dots[]) const
Return tangent component i of arguments.
#define SACADO_ENABLE_EXPR_FUNC(RETURN_TYPE)
KOKKOS_INLINE_FUNCTION const T * dx() const
Returns derivative array.
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors...
KOKKOS_INLINE_FUNCTION Expr(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
Wrapper for a generic expression template.
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
KOKKOS_INLINE_FUNCTION Expr()
Default constructor.
KOKKOS_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
KOKKOS_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
Initialize the derivative array.
KOKKOS_INLINE_FUNCTION void zero()
Zero out the derivative array.
expr expr dx(i)
KOKKOS_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
KOKKOS_INLINE_FUNCTION bool isActive() const
Return whether argument is active.
KOKKOS_INLINE_FUNCTION Expr(const Expr &x)
Copy constructor.
A tag for specializing Expr for SFad expressions.