Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_Fad_GeneralFad.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Sacado Package
4 //
5 // Copyright 2006 NTESS and the Sacado contributors.
6 // SPDX-License-Identifier: LGPL-2.1-or-later
7 //
8 // ***********************************************************************
9 //
10 // The forward-mode AD classes in Sacado are a derivative work of the
11 // expression template classes in the Fad package by Nicolas Di Cesare.
12 // The following banner is included in the original Fad source code:
13 //
14 // ************ DO NOT REMOVE THIS BANNER ****************
15 //
16 // Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
17 // http://www.ann.jussieu.fr/~dicesare
18 //
19 // CEMRACS 98 : C++ courses,
20 // templates : new C++ techniques
21 // for scientific computing
22 //
23 //********************************************************
24 //
25 // A short implementation ( not all operators and
26 // functions are overloaded ) of 1st order Automatic
27 // Differentiation in forward mode (FAD) using
28 // EXPRESSION TEMPLATES.
29 //
30 //********************************************************
31 // @HEADER
32 
33 #ifndef SACADO_FAD_GENERALFAD_HPP
34 #define SACADO_FAD_GENERALFAD_HPP
35 
37 
38 namespace Sacado {
39 
41  namespace Fad {
42 
43 #ifndef SACADO_FAD_DERIV_LOOP
44 #if defined(SACADO_VIEW_CUDA_HIERARCHICAL_DFAD) && !defined(SACADO_DISABLE_CUDA_IN_KOKKOS) && defined(__CUDA_ARCH__)
45 #define SACADO_FAD_DERIV_LOOP(I,SZ) for (int I=threadIdx.x; I<SZ; I+=blockDim.x)
46 #else
47 #define SACADO_FAD_DERIV_LOOP(I,SZ) for (int I=0; I<SZ; ++I)
48 #endif
49 #endif
50 
51 #ifndef SACADO_FAD_THREAD_SINGLE
52 #if (defined(SACADO_VIEW_CUDA_HIERARCHICAL) || defined(SACADO_VIEW_CUDA_HIERARCHICAL_DFAD)) && !defined(SACADO_DISABLE_CUDA_IN_KOKKOS) && defined(__CUDA_ARCH__)
53 #define SACADO_FAD_THREAD_SINGLE if (threadIdx.x == 0)
54 #else
55 #define SACADO_FAD_THREAD_SINGLE /* */
56 #endif
57 #endif
58 
60 
65  template <typename T, typename Storage>
66  class GeneralFad : public Storage {
67 
68  public:
69 
71  typedef typename RemoveConst<T>::type value_type;
72 
75 
80 
83  GeneralFad() : Storage(T(0.)) {}
84 
86 
89  template <typename S>
92  Storage(x) {}
93 
95 
99  GeneralFad(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) :
100  Storage(sz, x, zero_out) {}
101 
103 
109  GeneralFad(const int sz, const int i, const T & x) :
110  Storage(sz, x, InitDerivArray) {
111  this->fastAccessDx(i)=1.;
112  }
113 
116  GeneralFad(const Storage& s) : Storage(s) {}
117 
121  Storage(x) {}
122 
124  template <typename S>
127  Storage(x.size(), T(0.), NoInitDerivArray)
128  {
129  const int sz = x.size();
130 
131  if (sz) {
132  if (x.hasFastAccess())
134  this->fastAccessDx(i) = x.fastAccessDx(i);
135  else
137  this->fastAccessDx(i) = x.dx(i);
138  }
139 
140  this->val() = x.val();
141  }
142 
146 
148 
155  void diff(const int ith, const int n) {
156  if (this->size() != n)
157  this->resize(n);
158 
159  this->zero();
160  this->fastAccessDx(ith) = T(1.);
161  }
162 
165  void setUpdateValue(bool update_val) { }
166 
169  bool updateValue() const { return true; }
170 
173  void cache() const {}
174 
176  template <typename S>
178  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
179  typedef IsEqual<value_type> IE;
180  if (x.size() != this->size()) return false;
181  bool eq = IE::eval(x.val(), this->val());
182  const int sz = this->size();
184  eq = eq && IE::eval(x.dx(i), this->dx(i));
185  return eq;
186  }
187 
189 
194 
200  int availableSize() const { return this->length(); }
201 
204  bool hasFastAccess() const { return this->size()!=0; }
205 
208  bool isPassive() const { return this->size()==0; }
209 
212  void setIsConstant(bool is_const) {
213  if (is_const && this->size()!=0)
214  this->resize(0);
215  }
216 
218 
223 
225  template <typename S>
227  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator=(const S& v) {
228  this->val() = v;
229  if (this->size()) this->resize(0);
230  return *this;
231  }
232 
235  GeneralFad&
236  operator=(const GeneralFad& x) {
237  // Copy value & dx_
238  Storage::operator=(x);
239  return *this;
240  }
241 
243  template <typename S>
245  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator=(const Expr<S>& x) {
246  const int xsz = x.size();
247 
248  if (xsz != this->size())
249  this->resizeAndZero(xsz);
250 
251  const int sz = this->size();
252 
253  // For ViewStorage, the resize above may not in fact resize the
254  // derivative array, so it is possible that sz != xsz at this point.
255  // The only valid use case here is sz > xsz == 0, so we use sz in the
256  // assignment below
257 
258  if (sz) {
259  if (x.hasFastAccess()) {
261  this->fastAccessDx(i) = x.fastAccessDx(i);
262  }
263  else
265  this->fastAccessDx(i) = x.dx(i);
266  }
267 
268  this->val() = x.val();
269 
270  return *this;
271  }
272 
274 
279 
281  template <typename S>
283  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
284  this->val() += v;
285  return *this;
286  }
287 
289  template <typename S>
291  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator -= (const S& v) {
292  this->val() -= v;
293  return *this;
294  }
295 
297  template <typename S>
299  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator *= (const S& v) {
300  const int sz = this->size();
301  this->val() *= v;
303  this->fastAccessDx(i) *= v;
304  return *this;
305  }
306 
308  template <typename S>
310  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
311  const int sz = this->size();
312  this->val() /= v;
314  this->fastAccessDx(i) /= v;
315  return *this;
316  }
317 
320  GeneralFad& operator += (const GeneralFad& x) {
321  const int xsz = x.size(), sz = this->size();
322 
323 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
324  if ((xsz != sz) && (xsz != 0) && (sz != 0))
325  throw "Fad Error: Attempt to assign with incompatible sizes";
326 #endif
327 
328  if (xsz) {
329  if (sz) {
331  this->fastAccessDx(i) += x.fastAccessDx(i);
332  }
333  else {
334  this->resizeAndZero(xsz);
335  SACADO_FAD_DERIV_LOOP(i,xsz)
336  this->fastAccessDx(i) = x.fastAccessDx(i);
337  }
338  }
339 
340  this->val() += x.val();
341 
342  return *this;
343  }
344 
347  GeneralFad& operator -= (const GeneralFad& x) {
348  const int xsz = x.size(), sz = this->size();
349 
350 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
351  if ((xsz != sz) && (xsz != 0) && (sz != 0))
352  throw "Fad Error: Attempt to assign with incompatible sizes";
353 #endif
354 
355  if (xsz) {
356  if (sz) {
358  this->fastAccessDx(i) -= x.fastAccessDx(i);
359  }
360  else {
361  this->resizeAndZero(xsz);
362  SACADO_FAD_DERIV_LOOP(i,xsz)
363  this->fastAccessDx(i) = -x.fastAccessDx(i);
364  }
365  }
366 
367  this->val() -= x.val();
368 
369 
370  return *this;
371  }
372 
375  GeneralFad& operator *= (const GeneralFad& x) {
376  const int xsz = x.size(), sz = this->size();
377  T xval = x.val();
378  T v = this->val();
379 
380 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
381  if ((xsz != sz) && (xsz != 0) && (sz != 0))
382  throw "Fad Error: Attempt to assign with incompatible sizes";
383 #endif
384 
385  if (xsz) {
386  if (sz) {
388  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
389  }
390  else {
391  this->resizeAndZero(xsz);
392  SACADO_FAD_DERIV_LOOP(i,xsz)
393  this->fastAccessDx(i) = v*x.fastAccessDx(i);
394  }
395  }
396  else {
397  if (sz) {
399  this->fastAccessDx(i) *= xval;
400  }
401  }
402 
403  this->val() *= xval;
404 
405  return *this;
406  }
407 
410  GeneralFad& operator /= (const GeneralFad& x) {
411  const int xsz = x.size(), sz = this->size();
412  T xval = x.val();
413  T v = this->val();
414 
415 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
416  if ((xsz != sz) && (xsz != 0) && (sz != 0))
417  throw "Fad Error: Attempt to assign with incompatible sizes";
418 #endif
419 
420  if (xsz) {
421  if (sz) {
423  this->fastAccessDx(i) =
424  ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
425  }
426  else {
427  this->resizeAndZero(xsz);
428  SACADO_FAD_DERIV_LOOP(i,xsz)
429  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
430  }
431  }
432  else {
433  if (sz) {
435  this->fastAccessDx(i) /= xval;
436  }
437  }
438 
439  this->val() /= xval;
440 
441  return *this;
442  }
443 
445  template <typename S>
447  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator += (const Expr<S>& x) {
448  const int xsz = x.size(), sz = this->size();
449 
450 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
451  if ((xsz != sz) && (xsz != 0) && (sz != 0))
452  throw "Fad Error: Attempt to assign with incompatible sizes";
453 #endif
454 
455  if (xsz) {
456  if (sz) {
457  if (x.hasFastAccess())
459  this->fastAccessDx(i) += x.fastAccessDx(i);
460  else
462  this->fastAccessDx(i) += x.dx(i);
463  }
464  else {
465  this->resizeAndZero(xsz);
466  if (x.hasFastAccess())
467  SACADO_FAD_DERIV_LOOP(i,xsz)
468  this->fastAccessDx(i) = x.fastAccessDx(i);
469  else
470  SACADO_FAD_DERIV_LOOP(i,xsz)
471  this->fastAccessDx(i) = x.dx(i);
472  }
473  }
474 
475  this->val() += x.val();
476 
477  return *this;
478  }
479 
481  template <typename S>
483  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
484  const int xsz = x.size(), sz = this->size();
485 
486 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
487  if ((xsz != sz) && (xsz != 0) && (sz != 0))
488  throw "Fad Error: Attempt to assign with incompatible sizes";
489 #endif
490 
491  if (xsz) {
492  if (sz) {
493  if (x.hasFastAccess())
495  this->fastAccessDx(i) -= x.fastAccessDx(i);
496  else
498  this->fastAccessDx(i) -= x.dx(i);
499  }
500  else {
501  this->resizeAndZero(xsz);
502  if (x.hasFastAccess())
503  SACADO_FAD_DERIV_LOOP(i,xsz)
504  this->fastAccessDx(i) = -x.fastAccessDx(i);
505  else
506  SACADO_FAD_DERIV_LOOP(i,xsz)
507  this->fastAccessDx(i) = -x.dx(i);
508  }
509  }
510 
511  this->val() -= x.val();
512 
513 
514  return *this;
515  }
516 
518  template <typename S>
520  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
521  const int xsz = x.size(), sz = this->size();
522  T xval = x.val();
523  T v = this->val();
524 
525 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
526  if ((xsz != sz) && (xsz != 0) && (sz != 0))
527  throw "Fad Error: Attempt to assign with incompatible sizes";
528 #endif
529 
530  if (xsz) {
531  if (sz) {
532  if (x.hasFastAccess())
534  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
535  else
537  this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
538  }
539  else {
540  this->resizeAndZero(xsz);
541  if (x.hasFastAccess())
542  SACADO_FAD_DERIV_LOOP(i,xsz)
543  this->fastAccessDx(i) = v*x.fastAccessDx(i);
544  else
545  SACADO_FAD_DERIV_LOOP(i,xsz)
546  this->fastAccessDx(i) = v*x.dx(i);
547  }
548  }
549  else {
550  if (sz) {
552  this->fastAccessDx(i) *= xval;
553  }
554  }
555 
556  this->val() *= xval;
557 
558  return *this;
559  }
560 
562  template <typename S>
564  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
565  const int xsz = x.size(), sz = this->size();
566  T xval = x.val();
567  T v = this->val();
568 
569 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
570  if ((xsz != sz) && (xsz != 0) && (sz != 0))
571  throw "Fad Error: Attempt to assign with incompatible sizes";
572 #endif
573 
574  if (xsz) {
575  if (sz) {
576  if (x.hasFastAccess())
578  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
579  else
581  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
582  }
583  else {
584  this->resizeAndZero(xsz);
585  if (x.hasFastAccess())
586  SACADO_FAD_DERIV_LOOP(i,xsz)
587  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
588  else
589  SACADO_FAD_DERIV_LOOP(i,xsz)
590  this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
591  }
592  }
593  else {
594  if (sz) {
596  this->fastAccessDx(i) /= xval;
597  }
598  }
599 
600  this->val() /= xval;
601 
602  return *this;
603  }
604 
606 
607  }; // class GeneralFad
608 
609  } // namespace Fad
610 
611 } // namespace Sacado
612 
613 #endif // SACADO_FAD_GENERALFAD_HPP
Wrapper for a generic expression template.
SACADO_INLINE_FUNCTION GeneralFad(const GeneralFad &x)
Copy constructor.
expr expr dx(i)
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
#define SACADO_ENABLE_VALUE_CTOR_DECL
SACADO_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
#define SACADO_FAD_DERIV_LOOP(I, SZ)
SACADO_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
RemoveConst< T >::type value_type
Typename of values.
#define SACADO_ENABLE_EXPR_CTOR_DECL
SACADO_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
SACADO_INLINE_FUNCTION GeneralFad(const Storage &s)
Constructor with supplied storage s.
expr val()
#define T
Definition: Sacado_rad.hpp:553
#define SACADO_ENABLE_VALUE_FUNC(RETURN_TYPE)
SACADO_INLINE_FUNCTION GeneralFad(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
Base template specification for testing equivalence.
SACADO_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
Do not initialize the derivative array.
SACADO_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
SACADO_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
SACADO_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors...
SACADO_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
SACADO_INLINE_FUNCTION GeneralFad(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
Initialize the derivative array.
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 void cache() const
Cache values.
SACADO_INLINE_FUNCTION ~GeneralFad()
Destructor.
#define SACADO_INLINE_FUNCTION
SACADO_INLINE_FUNCTION GeneralFad(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
SACADO_INLINE_FUNCTION GeneralFad()
Default constructor.
Forward-mode AD class templated on the storage for the derivative array.