Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
248  void cache() const {}
249 
251  template <typename S>
253  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
254  typedef IsEqual<value_type> IE;
255  if (x.size() != this->size()) return false;
256  bool eq = IE::eval(x.val(), this->val());
257  for (int i=0; i<this->size(); i++)
258  eq = eq && IE::eval(x.dx(i), this->dx(i));
259  return eq;
260  }
261 
263 
268 
271  const T& val() const { return val_;}
272 
275  T& val() { return val_;}
276 
278 
283 
286  int size() const { return Num;}
287 
293  int availableSize() const { return Num; }
294 
297  bool hasFastAccess() const { return true; }
298 
301  bool isPassive() const { return false; }
302 
305  void setIsConstant(bool is_const) {}
306 
309  const T* dx() const { return &(dx_[0]);}
310 
313  const T& dx(int i) const { return dx_[i]; }
314 
317  T& fastAccessDx(int i) { return dx_[i];}
318 
321  const T& fastAccessDx(int i) const { return dx_[i];}
322 
325  void computePartials(const T& bar, value_type partials[]) const {
326  partials[0] = bar;
327  }
328 
331  void getTangents(int i, value_type dots[]) const {
332  dots[0] = this->dx_[i];
333  }
334 
336  template <int Arg>
338  bool isActive() const { return true; }
339 
341  template <int Arg>
343  const T& getTangent(int i) const { return this->dx_[i]; }
344 
347  const value_type* getDx(int j) const { return this->dx(); }
348 
350 
355 
357  template <typename S>
359  SACADO_ENABLE_VALUE_FUNC(Expr&) operator=(const S& v) {
360  val_ = v;
361  ss_array<T>::zero(dx_, Num);
362  return *this;
363  }
364 
367  Expr& operator=(const Expr& x) {
368  if (this != &x) {
369  // Copy value
370  val_ = x.val_;
371 
372  // Copy dx_
373  //ss_array<T>::copy(x.dx_, dx_, Num);
374  for (int i=0; i<Num; i++)
375  dx_[i] = x.dx_[i];
376  }
377  return *this;
378  }
379 
381  template <typename S>
383  SACADO_ENABLE_EXPR_FUNC(Expr&) operator=(const Expr<S>& x) {
384 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
385  if (x.size() != Num)
386  throw "SELRFad::operator=() Error: Attempt to assign with incompatible sizes";
387 #endif
388 
389  // Compute partials
390  LocalAccumOp< Expr<S> > op(x);
391 
392  // Compute each tangent direction
393  for(int i=0; i<Num; ++i) {
394  op.t = T(0.);
395  op.i = i;
396 
397  // Automatically unrolled loop that computes
398  // for (int j=0; j<N; j++)
399  // op.t += op.partials[j] * x.getTangent<j>(i);
401 
402  dx_[i] = op.t;
403  }
404 
405  // Compute value
406  val_ = x.val();
407 
408  return *this;
409  }
410 
412 
417 
419  template <typename S>
421  SACADO_ENABLE_VALUE_FUNC(Expr&) operator += (const S& v) {
422  this->val() += v;
423  return *this;
424  }
425 
427  template <typename S>
429  SACADO_ENABLE_VALUE_FUNC(Expr&) operator -= (const S& v) {
430  this->val() -= v;
431  return *this;
432  }
433 
435  template <typename S>
437  SACADO_ENABLE_VALUE_FUNC(Expr&) operator *= (const S& v) {
438  this->val() *= v;
439  for (int i=0; i<Num; ++i)
440  dx_[i] *= v;
441  return *this;
442  }
443 
445  template <typename S>
447  SACADO_ENABLE_VALUE_FUNC(Expr&) operator /= (const S& v) {
448  this->val() /= v;
449  for (int i=0; i<Num; ++i)
450  dx_[i] /= v;
451  return *this;
452  }
453 
455  template <typename S>
457  SACADO_ENABLE_EXPR_FUNC(Expr&) operator += (const Expr<S>& x) {
458 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
459  if (x.size() != Num)
460  throw "SELRFad::operator+=() Error: Attempt to assign with incompatible sizes";
461 #endif
462 
463  // Compute partials
464  LocalAccumOp< Expr<S> > op(x);
465 
466  // Compute each tangent direction
467  for(int i=0; i<Num; ++i) {
468  op.t = T(0.);
469  op.i = i;
470 
471  // Automatically unrolled loop that computes
472  // for (int j=0; j<N; j++)
473  // op.t += op.partials[j] * x.getTangent<j>(i);
475 
476  dx_[i] += op.t;
477  }
478 
479  // Compute value
480  val_ += x.val();
481 
482  return *this;
483  }
484 
486  template <typename S>
488  SACADO_ENABLE_EXPR_FUNC(Expr&) operator -= (const Expr<S>& x) {
489 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
490  if (x.size() != Num)
491  throw "SELRFad::operator-=() Error: Attempt to assign with incompatible sizes";
492 #endif
493 
494  // Compute partials
495  LocalAccumOp< Expr<S> > op(x);
496 
497  // Compute each tangent direction
498  for(int i=0; i<Num; ++i) {
499  op.t = T(0.);
500  op.i = i;
501 
502  // Automatically unrolled loop that computes
503  // for (int j=0; j<N; j++)
504  // op.t += op.partials[j] * x.getTangent<j>(i);
506 
507  dx_[i] -= op.t;
508  }
509 
510  // Compute value
511  val_ -= x.val();
512 
513  return *this;
514  }
515 
517  template <typename S>
519  SACADO_ENABLE_EXPR_FUNC(Expr&) operator *= (const Expr<S>& x) {
520  T xval = x.val();
521 
522 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
523  if (x.size() != Num)
524  throw "SELRFad::operator*=() Error: Attempt to assign with incompatible sizes";
525 #endif
526 
527  // Compute partials
528  LocalAccumOp< Expr<S> > op(x);
529 
530  // Compute each tangent direction
531  for(int i=0; i<Num; ++i) {
532  op.t = T(0.);
533  op.i = i;
534 
535  // Automatically unrolled loop that computes
536  // for (int j=0; j<N; j++)
537  // op.t += op.partials[j] * x.getTangent<j>(i);
539 
540  dx_[i] = val_ * op.t + dx_[i] * xval;
541  }
542 
543  // Compute value
544  val_ *= xval;
545 
546  return *this;
547  }
548 
550  template <typename S>
552  SACADO_ENABLE_EXPR_FUNC(Expr&) operator /= (const Expr<S>& x) {
553  T xval = x.val();
554 
555 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
556  if (x.size() != Num)
557  throw "SELRFad::operator/=() Error: Attempt to assign with incompatible sizes";
558 #endif
559 
560  // Compute partials
561  LocalAccumOp< Expr<S> > op(x);
562 
563  T xval2 = xval*xval;
564 
565  // Compute each tangent direction
566  for(int i=0; i<Num; ++i) {
567  op.t = T(0.);
568  op.i = i;
569 
570  // Automatically unrolled loop that computes
571  // for (int j=0; j<N; j++)
572  // op.t += op.partials[j] * x.getTangent<j>(i);
574 
575  dx_[i] = (dx_[i] * xval - val_ * op.t) / xval2;
576  }
577 
578  // Compute value
579  val_ /= xval;
580 
581  return *this;
582  }
583 
585 
586  protected:
587 
589  T dx_[Num];
590 
593 
594  // Functor for mpl::for_each to compute the local accumulation
595  // of a tangent derivative
596  //
597  // We use getTangent<>() to get dx components from expression
598  // arguments instead of getting the argument directly or extracting
599  // the dx array due to striding in ViewFad (or could use striding
600  // directly here if we need to get dx array).
601  template <typename ExprT>
602  struct LocalAccumOp {
603  typedef typename ExprT::value_type value_type;
604  static const int N = ExprT::num_args;
605  const ExprT& x;
606  mutable value_type t;
607  value_type partials[N];
608  int i;
610  LocalAccumOp(const ExprT& x_) : x(x_) {
611  x.computePartials(value_type(1.), partials);
612  }
614  void getTangents(int i_) { i = i_; }
615  template <typename ArgT>
617  void operator () (ArgT arg) const {
618  const int Arg = ArgT::value;
619  if (x.template isActive<Arg>())
620  t += partials[Arg] * x.template getTangent<Arg>(i);
621  }
622  };
623 
624  }; // class Expr<SFadExprTag>
625 
626  } // namespace ELRFad
627 
628 } // namespace Sacado
629 
630 #define FAD_NS ELRFad
631 #include "Sacado_Fad_SFad_tmpl.hpp"
632 #undef FAD_NS
633 
634 #include "Sacado_ELRFad_ViewFad.hpp"
635 #include "Sacado_ELRFad_Ops.hpp"
636 
637 #endif // SACADO_ELRFAD_SFAD_HPP
void f()
KOKKOS_INLINE_FUNCTION Expr(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
expr expr dx(i)
KOKKOS_INLINE_FUNCTION int size() const
Returns number of derivative components.
SFad< value_type, Num > base_expr_type
Typename of base-expressions.
KOKKOS_INLINE_FUNCTION T & val()
Returns value.
#define SACADO_ENABLE_VALUE_CTOR_DECL
RemoveConst< T >::type value_type
Typename of values.
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 Expr(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
KOKKOS_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
KOKKOS_INLINE_FUNCTION const T * dx() const
Returns derivative array.
expr val()
#define KOKKOS_INLINE_FUNCTION
#define T
Definition: Sacado_rad.hpp:573
KOKKOS_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
KOKKOS_INLINE_FUNCTION T & fastAccessDx(int i)
Returns derivative component i without bounds checking.
#define SACADO_ENABLE_VALUE_FUNC(RETURN_TYPE)
KOKKOS_INLINE_FUNCTION bool isActive() const
Return whether argument is active.
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 const T & fastAccessDx(int i) const
Returns derivative component i without bounds checking.
KOKKOS_INLINE_FUNCTION ~Expr()
Destructor.
const int N
KOKKOS_INLINE_FUNCTION void diff(const int ith, const int n)
Set Fad object as the ith independent variable.
#define SACADO_ENABLE_EXPR_FUNC(RETURN_TYPE)
KOKKOS_INLINE_FUNCTION const T & val() const
Returns value.
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 const T & dx(int i) const
Returns derivative component i with bounds checking.
KOKKOS_INLINE_FUNCTION void cache() const
Cache values.
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.
KOKKOS_INLINE_FUNCTION const value_type * getDx(int j) const
Get dx array.
KOKKOS_INLINE_FUNCTION const T & getTangent(int i) const
Return tangent component i of argument Arg.
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 void computePartials(const T &bar, value_type partials[]) const
Return partials w.r.t. arguments.
KOKKOS_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
KOKKOS_INLINE_FUNCTION Expr(const Expr &x)
Copy constructor.
A tag for specializing Expr for SFad expressions.
KOKKOS_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
KOKKOS_INLINE_FUNCTION void getTangents(int i, value_type dots[]) const
Return tangent component i of arguments.