Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_ELRCacheFad_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_ELRCACHEFAD_SFAD_HPP
53 #define SACADO_ELRCACHEFAD_SFAD_HPP
54 
58 #include "Sacado_mpl_range_c.hpp"
59 #include "Sacado_mpl_for_each.hpp"
60 
61 namespace Sacado {
62 
63  namespace ELRCacheFad {
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 "ELRCacheFad::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 "ELRCacheFad::SFad() Error: Supplied derivative dimension does not match compile time length.";
146  if (i >= Num)
147  throw "ELRCacheFad::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 "ELRCacheFad::SFad() Error: Attempt to assign with incompatible sizes";
169 #endif
170 
171  x.cache();
172 
173  this->val() = x.val();
174 
175  // Compute partials
176  LocalAccumOp< Expr<S> > op(x);
177 
178  // Compute each tangent direction
179  for(int i=0; i<Num; ++i) {
180  op.t = T(0.);
181  op.i = i;
182 
183  // Automatically unrolled loop that computes
184  // for (int j=0; j<N; j++)
185  // op.t += op.partials[j] * x.getTangent<j>(i);
187 
188  dx_[i] = op.t;
189 
190  }
191  }
192 
195  ~Expr() {}
196 
198 
205  void diff(const int ith, const int n) {
206 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
207  if (n != Num)
208  throw "ELRCacheFad::diff() Error: Supplied derivative dimension does not match compile time length.";
209 #endif
210 
211  ss_array<T>::zero(dx_, Num);
212  dx_[ith] = T(1.);
213  }
214 
216 
221  void resize(int sz) {
222 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
223  if (sz != Num)
224  throw "ELRCacheFad::resize() Error: Cannot resize fixed derivative array dimension";
225 #endif
226  }
227 
229 
234  void expand(int sz) { resize(sz); }
235 
238  void zero() { ss_array<T>::zero(dx_, Num); }
239 
242  void setUpdateValue(bool update_val) { }
243 
246  bool updateValue() const { return true; }
247 
250  void cache() const {}
251 
253  template <typename S>
255  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
256  typedef IsEqual<value_type> IE;
257  if (x.size() != this->size()) return false;
258  bool eq = IE::eval(x.val(), this->val());
259  for (int i=0; i<this->size(); i++)
260  eq = eq && IE::eval(x.dx(i), this->dx(i));
261  return eq;
262  }
263 
265 
270 
273  const T& val() const { return val_;}
274 
277  T& val() { return val_;}
278 
280 
285 
288  int size() const { return Num;}
289 
295  int availableSize() const { return Num; }
296 
299  bool hasFastAccess() const { return true; }
300 
303  bool isPassive() const { return false; }
304 
307  void setIsConstant(bool is_const) {}
308 
311  const T* dx() const { return &(dx_[0]);}
312 
315  const T& dx(int i) const { return dx_[i]; }
316 
319  T& fastAccessDx(int i) { return dx_[i];}
320 
323  const T& fastAccessDx(int i) const { return dx_[i];}
324 
327  void computePartials(const T& bar, value_type partials[]) const {
328  partials[0] = bar;
329  }
330 
333  void getTangents(int i, value_type dots[]) const {
334  dots[0] = this->dx_[i];
335  }
336 
338  template <int Arg>
340  bool isActive() const { return true; }
341 
343  template <int Arg>
345  T getTangent(int i) const { return this->dx_[i]; }
346 
349  const value_type* getDx(int j) const { return this->dx(); }
350 
352 
357 
359  template <typename S>
361  SACADO_ENABLE_VALUE_FUNC(Expr&) operator=(const S& v) {
362  val_ = v;
363  ss_array<T>::zero(dx_, Num);
364  return *this;
365  }
366 
369  Expr& operator=(const Expr& x) {
370  if (this != &x) {
371  // Copy value
372  val_ = x.val_;
373 
374  // Copy dx_
375  //ss_array<T>::copy(x.dx_, dx_, Num);
376  for (int i=0; i<Num; i++)
377  dx_[i] = x.dx_[i];
378  }
379  return *this;
380  }
381 
383  template <typename S>
385  SACADO_ENABLE_EXPR_FUNC(Expr&) operator=(const Expr<S>& x) {
386 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
387  if (x.size() != Num)
388  throw "ELRCacheFad::operator=() Error: Attempt to assign with incompatible sizes";
389 #endif
390 
391  x.cache();
392 
393  // Compute partials
394  LocalAccumOp< Expr<S> > op(x);
395 
396  // Compute each tangent direction
397  for(int i=0; i<Num; ++i) {
398  op.t = T(0.);
399  op.i = i;
400 
401  // Automatically unrolled loop that computes
402  // for (int j=0; j<N; j++)
403  // op.t += op.partials[j] * x.getTangent<j>(i);
405 
406  dx_[i] = op.t;
407  }
408 
409  this->val() = x.val();
410 
411  return *this;
412  }
413 
415 
420 
422  template <typename S>
424  SACADO_ENABLE_VALUE_FUNC(Expr&) operator += (const S& v) {
425  this->val() += v;
426  return *this;
427  }
428 
430  template <typename S>
432  SACADO_ENABLE_VALUE_FUNC(Expr&) operator -= (const S& v) {
433  this->val() -= v;
434  return *this;
435  }
436 
438  template <typename S>
440  SACADO_ENABLE_VALUE_FUNC(Expr&) operator *= (const S& v) {
441  this->val() *= v;
442  for (int i=0; i<Num; ++i)
443  dx_[i] *= v;
444  return *this;
445  }
446 
448  template <typename S>
450  SACADO_ENABLE_VALUE_FUNC(Expr&) operator /= (const S& v) {
451  this->val() /= v;
452  for (int i=0; i<Num; ++i)
453  dx_[i] /= v;
454  return *this;
455  }
456 
458  template <typename S>
460  SACADO_ENABLE_EXPR_FUNC(Expr&) operator += (const Expr<S>& x) {
461 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
462  if (x.size() != Num)
463  throw "ELRCacheFad::operator+=() Error: Attempt to assign with incompatible sizes";
464 #endif
465 
466  x.cache();
467 
468  // Compute partials
469  LocalAccumOp< Expr<S> > op(x);
470 
471  // Compute each tangent direction
472  for(int i=0; i<Num; ++i) {
473  op.t = T(0.);
474  op.i = i;
475 
476  // Automatically unrolled loop that computes
477  // for (int j=0; j<N; j++)
478  // op.t += op.partials[j] * x.getTangent<j>(i);
480 
481  dx_[i] += op.t;
482  }
483 
484  this->val() += x.val();
485 
486  return *this;
487  }
488 
490  template <typename S>
492  SACADO_ENABLE_EXPR_FUNC(Expr&) operator -= (const Expr<S>& x) {
493 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
494  if (x.size() != Num)
495  throw "ELRCacheFad::operator-=() Error: Attempt to assign with incompatible sizes";
496 #endif
497 
498  x.cache();
499 
500  // Compute partials
501  LocalAccumOp< Expr<S> > op(x);
502 
503  // Compute each tangent direction
504  for(int i=0; i<Num; ++i) {
505  op.t = T(0.);
506  op.i = i;
507 
508  // Automatically unrolled loop that computes
509  // for (int j=0; j<N; j++)
510  // op.t += op.partials[j] * x.getTangent<j>(i);
512 
513  dx_[i] -= op.t;
514  }
515 
516  this->val() -= x.val();
517 
518  return *this;
519  }
520 
522  template <typename S>
524  SACADO_ENABLE_EXPR_FUNC(Expr&) operator *= (const Expr<S>& x) {
525  x.cache();
526 
527  T xval = x.val();
528 
529 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
530  if (x.size() != Num)
531  throw "ELRCacheFad::operator*=() Error: Attempt to assign with incompatible sizes";
532 #endif
533 
534  // Compute partials
535  LocalAccumOp< Expr<S> > op(x);
536 
537  // Compute each tangent direction
538  for(int i=0; i<Num; ++i) {
539  op.t = T(0.);
540  op.i = i;
541 
542  // Automatically unrolled loop that computes
543  // for (int j=0; j<N; j++)
544  // op.t += op.partials[j] * x.getTangent<j>(i);
546 
547  dx_[i] = val_ * op.t + dx_[i] * xval;
548  }
549 
550  // Compute value
551  val_ *= xval;
552 
553  return *this;
554  }
555 
557  template <typename S>
559  SACADO_ENABLE_EXPR_FUNC(Expr&) operator /= (const Expr<S>& x) {
560  x.cache();
561 
562  T xval = x.val();
563 
564 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
565  if (x.size() != Num)
566  throw "ELRCacheFad::operator/=() Error: Attempt to assign with incompatible sizes";
567 #endif
568 
569  // Compute partials
570  LocalAccumOp< Expr<S> > op(x);
571 
572  T xval2 = xval*xval;
573 
574  // Compute each tangent direction
575  for(int i=0; i<Num; ++i) {
576  op.t = T(0.);
577  op.i = i;
578 
579  // Automatically unrolled loop that computes
580  // for (int j=0; j<N; j++)
581  // op.t += op.partials[j] * x.getTangent<j>(i);
583 
584  dx_[i] = (dx_[i] * xval - val_ * op.t) / xval2;
585  }
586 
587  // Compute value
588  val_ /= xval;
589 
590  return *this;
591  }
592 
594 
595  protected:
596 
598  T dx_[Num];
599 
602 
603  // Functor for mpl::for_each to compute the local accumulation
604  // of a tangent derivative
605  //
606  // We use getTangent<>() to get dx components from expression
607  // arguments instead of getting the argument directly or extracting
608  // the dx array due to striding in ViewFad (or could use striding
609  // directly here if we need to get dx array).
610  template <typename ExprT>
611  struct LocalAccumOp {
612  typedef typename ExprT::value_type value_type;
613  static const int N = ExprT::num_args;
614  const ExprT& x;
615  mutable value_type t;
616  value_type partials[N];
617  int i;
619  LocalAccumOp(const ExprT& x_) : x(x_) {
620  x.computePartials(value_type(1.), partials);
621  }
622  template <typename ArgT>
624  void operator () (ArgT arg) const {
625  const int Arg = ArgT::value;
626  if (x.template isActive<Arg>())
627  t += partials[Arg] * x.template getTangent<Arg>(i);
628  }
629  };
630 
631  }; // class Expr<SFadExprTag>
632 
633  } // namespace ELRCacheFad
634 
635 } // namespace Sacado
636 
637 #define FAD_NS ELRCacheFad
638 #include "Sacado_Fad_SFad_tmpl.hpp"
639 #undef FAD_NS
640 
643 
644 #endif // SACADO_ELRCACHEFAD_SFAD_HPP
SACADO_INLINE_FUNCTION Expr(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
void f()
SACADO_INLINE_FUNCTION Expr(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
SACADO_INLINE_FUNCTION const value_type * getDx(int j) const
Get dx array.
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
SACADO_INLINE_FUNCTION void cache() const
Cache values.
SACADO_INLINE_FUNCTION const T & fastAccessDx(int i) const
Returns derivative component i without bounds checking.
expr expr dx(i)
#define SACADO_ENABLE_VALUE_CTOR_DECL
SACADO_INLINE_FUNCTION SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr< S > &x) const
Returns whether two Fad objects have the same values.
SACADO_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
SFad< value_type, Num > base_expr_type
Typename of base-expressions.
SACADO_INLINE_FUNCTION Expr()
Default constructor.
expr bar
SACADO_INLINE_FUNCTION void computePartials(const T &bar, value_type partials[]) const
Return partials w.r.t. arguments.
#define SACADO_ENABLE_EXPR_CTOR_DECL
SACADO_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
SACADO_INLINE_FUNCTION void expand(int sz)
Expand derivative array to size sz.
SACADO_INLINE_FUNCTION T & val()
Returns value.
SACADO_INLINE_FUNCTION void resize(int sz)
Resize derivative array to length sz.
SACADO_INLINE_FUNCTION Expr(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
SACADO_INLINE_FUNCTION Expr(const Expr &x)
Copy constructor.
expr val()
#define T
Definition: Sacado_rad.hpp:573
SACADO_INLINE_FUNCTION void zero()
Zero out the derivative array.
SACADO_INLINE_FUNCTION void diff(const int ith, const int n)
Set Fad object as the ith independent variable.
#define SACADO_ENABLE_VALUE_FUNC(RETURN_TYPE)
Base template specification for testing equivalence.
SACADO_INLINE_FUNCTION const T & val() const
Returns value.
static double x_
Wrapper for a generic expression template.
SACADO_INLINE_FUNCTION int size() const
Returns number of derivative components.
RemoveConst< T >::type value_type
Typename of values.
SACADO_INLINE_FUNCTION Expr(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
const int N
SACADO_INLINE_FUNCTION bool isActive() const
Return whether argument is active.
#define SACADO_ENABLE_EXPR_FUNC(RETURN_TYPE)
SACADO_INLINE_FUNCTION const T & dx(int i) const
Returns derivative component i with bounds checking.
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors...
int value
SACADO_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
SACADO_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
SACADO_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
static SACADO_INLINE_FUNCTION void zero(T *dest, int sz)
Zero out array dest of length sz.
Initialize the derivative array.
SACADO_INLINE_FUNCTION T & fastAccessDx(int i)
Returns derivative component i without bounds checking.
SACADO_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
#define SACADO_INLINE_FUNCTION
SACADO_INLINE_FUNCTION void getTangents(int i, value_type dots[]) const
Return tangent component i of arguments.
SACADO_INLINE_FUNCTION const T * dx() const
Returns derivative array.
SACADO_INLINE_FUNCTION T getTangent(int i) const
Return tangent component i of argument Arg.
A tag for specializing Expr for SFad expressions.