Sacado Package Browser (Single Doxygen Collection)  Version of the Day
Sacado_Fad_GeneralFad.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_FAD_GENERALFAD_HPP
53 #define SACADO_FAD_GENERALFAD_HPP
54 
56 
57 namespace Sacado {
58 
60  namespace Fad {
61 
63 
68  template <typename T, typename Storage>
69  class GeneralFad : public Storage {
70 
71  public:
72 
74  typedef typename RemoveConst<T>::type value_type;
75 
78 
83 
86  GeneralFad() : Storage(T(0.)) {}
87 
89 
92  template <typename S>
95  Storage(x) {}
96 
98 
102  GeneralFad(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) :
103  Storage(sz, x, zero_out) {}
104 
106 
112  GeneralFad(const int sz, const int i, const T & x) :
113  Storage(sz, x) {
114  this->fastAccessDx(i)=1.;
115  }
116 
119  GeneralFad(const Storage& s) : Storage(s) {}
120 
123  GeneralFad(const GeneralFad& x) :
124  Storage(x) {}
125 
127  template <typename S>
130  Storage(x.size(), T(0.))
131  {
132  const int sz = x.size();
133 
134  if (sz) {
135  if (x.hasFastAccess())
136  for(int i=0; i<sz; ++i)
137  this->fastAccessDx(i) = x.fastAccessDx(i);
138  else
139  for(int i=0; i<sz; ++i)
140  this->fastAccessDx(i) = x.dx(i);
141  }
142 
143  this->val() = x.val();
144  }
145 
149 
151 
158  void diff(const int ith, const int n) {
159  if (this->size() != n)
160  this->resize(n);
161 
162  this->zero();
163  this->fastAccessDx(ith) = T(1.);
164  }
165 
168  void setUpdateValue(bool update_val) { }
169 
172  bool updateValue() const { return true; }
173 
175  template <typename S>
177  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
178  typedef IsEqual<value_type> IE;
179  if (x.size() != this->size()) return false;
180  bool eq = IE::eval(x.val(), this->val());
181  for (int i=0; i<this->size(); i++)
182  eq = eq && IE::eval(x.dx(i), this->dx(i));
183  return eq;
184  }
185 
187 
192 
198  int availableSize() const { return this->length(); }
199 
202  bool hasFastAccess() const { return this->size()!=0; }
203 
206  bool isPassive() const { return this->size()==0; }
207 
210  void setIsConstant(bool is_const) {
211  if (is_const && this->size()!=0)
212  this->resize(0);
213  }
214 
216 
221 
223  template <typename S>
225  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator=(const S& v) {
226  this->val() = v;
227  if (this->size()) this->resize(0);
228  return *this;
229  }
230 
233  GeneralFad&
234  operator=(const GeneralFad& x) {
235  // Copy value & dx_
236  Storage::operator=(x);
237  return *this;
238  }
239 
241  template <typename S>
243  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator=(const Expr<S>& x) {
244  const int xsz = x.size();
245 
246  if (xsz != this->size())
247  this->resizeAndZero(xsz);
248 
249  const int sz = this->size();
250 
251  // For ViewStorage, the resize above may not in fact resize the
252  // derivative array, so it is possible that sz != xsz at this point.
253  // The only valid use case here is sz > xsz == 0, so we use sz in the
254  // assignment below
255 
256  if (sz) {
257  if (x.hasFastAccess())
258  for(int i=0; i<sz; ++i)
259  this->fastAccessDx(i) = x.fastAccessDx(i);
260  else
261  for(int i=0; i<sz; ++i)
262  this->fastAccessDx(i) = x.dx(i);
263  }
264 
265  this->val() = x.val();
266 
267  return *this;
268  }
269 
271 
276 
278  template <typename S>
280  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
281  this->val() += v;
282  return *this;
283  }
284 
286  template <typename S>
288  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator -= (const S& v) {
289  this->val() -= v;
290  return *this;
291  }
292 
294  template <typename S>
296  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator *= (const S& v) {
297  const int sz = this->size();
298  this->val() *= v;
299  for (int i=0; i<sz; ++i)
300  this->fastAccessDx(i) *= v;
301  return *this;
302  }
303 
305  template <typename S>
307  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
308  const int sz = this->size();
309  this->val() /= v;
310  for (int i=0; i<sz; ++i)
311  this->fastAccessDx(i) /= v;
312  return *this;
313  }
314 
317  GeneralFad& operator += (const GeneralFad& x) {
318  const int xsz = x.size(), sz = this->size();
319 
320 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
321  if ((xsz != sz) && (xsz != 0) && (sz != 0))
322  throw "Fad Error: Attempt to assign with incompatible sizes";
323 #endif
324 
325  if (xsz) {
326  if (sz) {
327  for (int i=0; i<sz; ++i)
328  this->fastAccessDx(i) += x.fastAccessDx(i);
329  }
330  else {
331  this->resizeAndZero(xsz);
332  for (int i=0; i<xsz; ++i)
333  this->fastAccessDx(i) = x.fastAccessDx(i);
334  }
335  }
336 
337  this->val() += x.val();
338 
339  return *this;
340  }
341 
344  GeneralFad& operator -= (const GeneralFad& x) {
345  const int xsz = x.size(), sz = this->size();
346 
347 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
348  if ((xsz != sz) && (xsz != 0) && (sz != 0))
349  throw "Fad Error: Attempt to assign with incompatible sizes";
350 #endif
351 
352  if (xsz) {
353  if (sz) {
354  for(int i=0; i<sz; ++i)
355  this->fastAccessDx(i) -= x.fastAccessDx(i);
356  }
357  else {
358  this->resizeAndZero(xsz);
359  for(int i=0; i<xsz; ++i)
360  this->fastAccessDx(i) = -x.fastAccessDx(i);
361  }
362  }
363 
364  this->val() -= x.val();
365 
366 
367  return *this;
368  }
369 
372  GeneralFad& operator *= (const GeneralFad& x) {
373  const int xsz = x.size(), sz = this->size();
374  T xval = x.val();
375  T v = this->val();
376 
377 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
378  if ((xsz != sz) && (xsz != 0) && (sz != 0))
379  throw "Fad Error: Attempt to assign with incompatible sizes";
380 #endif
381 
382  if (xsz) {
383  if (sz) {
384  for(int i=0; i<sz; ++i)
385  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
386  }
387  else {
388  this->resizeAndZero(xsz);
389  for(int i=0; i<xsz; ++i)
390  this->fastAccessDx(i) = v*x.fastAccessDx(i);
391  }
392  }
393  else {
394  if (sz) {
395  for (int i=0; i<sz; ++i)
396  this->fastAccessDx(i) *= xval;
397  }
398  }
399 
400  this->val() *= xval;
401 
402  return *this;
403  }
404 
407  GeneralFad& operator /= (const GeneralFad& x) {
408  const int xsz = x.size(), sz = this->size();
409  T xval = x.val();
410  T v = this->val();
411 
412 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
413  if ((xsz != sz) && (xsz != 0) && (sz != 0))
414  throw "Fad Error: Attempt to assign with incompatible sizes";
415 #endif
416 
417  if (xsz) {
418  if (sz) {
419  for(int i=0; i<sz; ++i)
420  this->fastAccessDx(i) =
421  ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
422  }
423  else {
424  this->resizeAndZero(xsz);
425  for(int i=0; i<xsz; ++i)
426  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
427  }
428  }
429  else {
430  if (sz) {
431  for (int i=0; i<sz; ++i)
432  this->fastAccessDx(i) /= xval;
433  }
434  }
435 
436  this->val() /= xval;
437 
438  return *this;
439  }
440 
442  template <typename S>
444  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator += (const Expr<S>& x) {
445  const int xsz = x.size(), sz = this->size();
446 
447 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
448  if ((xsz != sz) && (xsz != 0) && (sz != 0))
449  throw "Fad Error: Attempt to assign with incompatible sizes";
450 #endif
451 
452  if (xsz) {
453  if (sz) {
454  if (x.hasFastAccess())
455  for (int i=0; i<sz; ++i)
456  this->fastAccessDx(i) += x.fastAccessDx(i);
457  else
458  for (int i=0; i<sz; ++i)
459  this->fastAccessDx(i) += x.dx(i);
460  }
461  else {
462  this->resizeAndZero(xsz);
463  if (x.hasFastAccess())
464  for (int i=0; i<xsz; ++i)
465  this->fastAccessDx(i) = x.fastAccessDx(i);
466  else
467  for (int i=0; i<xsz; ++i)
468  this->fastAccessDx(i) = x.dx(i);
469  }
470  }
471 
472  this->val() += x.val();
473 
474  return *this;
475  }
476 
478  template <typename S>
480  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
481  const int xsz = x.size(), sz = this->size();
482 
483 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
484  if ((xsz != sz) && (xsz != 0) && (sz != 0))
485  throw "Fad Error: Attempt to assign with incompatible sizes";
486 #endif
487 
488  if (xsz) {
489  if (sz) {
490  if (x.hasFastAccess())
491  for(int i=0; i<sz; ++i)
492  this->fastAccessDx(i) -= x.fastAccessDx(i);
493  else
494  for (int i=0; i<sz; ++i)
495  this->fastAccessDx(i) -= x.dx(i);
496  }
497  else {
498  this->resizeAndZero(xsz);
499  if (x.hasFastAccess())
500  for(int i=0; i<xsz; ++i)
501  this->fastAccessDx(i) = -x.fastAccessDx(i);
502  else
503  for (int i=0; i<xsz; ++i)
504  this->fastAccessDx(i) = -x.dx(i);
505  }
506  }
507 
508  this->val() -= x.val();
509 
510 
511  return *this;
512  }
513 
515  template <typename S>
517  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
518  const int xsz = x.size(), sz = this->size();
519  T xval = x.val();
520  T v = this->val();
521 
522 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
523  if ((xsz != sz) && (xsz != 0) && (sz != 0))
524  throw "Fad Error: Attempt to assign with incompatible sizes";
525 #endif
526 
527  if (xsz) {
528  if (sz) {
529  if (x.hasFastAccess())
530  for(int i=0; i<sz; ++i)
531  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
532  else
533  for (int i=0; i<sz; ++i)
534  this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
535  }
536  else {
537  this->resizeAndZero(xsz);
538  if (x.hasFastAccess())
539  for(int i=0; i<xsz; ++i)
540  this->fastAccessDx(i) = v*x.fastAccessDx(i);
541  else
542  for (int i=0; i<xsz; ++i)
543  this->fastAccessDx(i) = v*x.dx(i);
544  }
545  }
546  else {
547  if (sz) {
548  for (int i=0; i<sz; ++i)
549  this->fastAccessDx(i) *= xval;
550  }
551  }
552 
553  this->val() *= xval;
554 
555  return *this;
556  }
557 
559  template <typename S>
561  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
562  const int xsz = x.size(), sz = this->size();
563  T xval = x.val();
564  T v = this->val();
565 
566 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
567  if ((xsz != sz) && (xsz != 0) && (sz != 0))
568  throw "Fad Error: Attempt to assign with incompatible sizes";
569 #endif
570 
571  if (xsz) {
572  if (sz) {
573  if (x.hasFastAccess())
574  for(int i=0; i<sz; ++i)
575  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
576  else
577  for (int i=0; i<sz; ++i)
578  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
579  }
580  else {
581  this->resizeAndZero(xsz);
582  if (x.hasFastAccess())
583  for(int i=0; i<xsz; ++i)
584  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
585  else
586  for (int i=0; i<xsz; ++i)
587  this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
588  }
589  }
590  else {
591  if (sz) {
592  for (int i=0; i<sz; ++i)
593  this->fastAccessDx(i) /= xval;
594  }
595  }
596 
597  this->val() /= xval;
598 
599  return *this;
600  }
601 
603 
604  }; // class GeneralFad
605 
606  } // namespace Fad
607 
608 } // namespace Sacado
609 
610 #endif // SACADO_FAD_GENERALFAD_HPP
Wrapper for a generic expression template.
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 ~GeneralFad()
Destructor.
expr val()
KOKKOS_INLINE_FUNCTION GeneralFad(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
#define SACADO_ENABLE_VALUE_CTOR_DECL
KOKKOS_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
RemoveConst< T >::type value_type
Typename of values.
#define SACADO_ENABLE_EXPR_CTOR_DECL
KOKKOS_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
KOKKOS_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
KOKKOS_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
#define KOKKOS_INLINE_FUNCTION
#define T
Definition: Sacado_rad.hpp:573
KOKKOS_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
KOKKOS_INLINE_FUNCTION GeneralFad()
Default constructor.
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
#define SACADO_ENABLE_VALUE_FUNC(RETURN_TYPE)
KOKKOS_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.
KOKKOS_INLINE_FUNCTION GeneralFad(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors...
KOKKOS_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
KOKKOS_INLINE_FUNCTION GeneralFad(const Storage &s)
Constructor with supplied storage s.
Initialize the derivative array.
expr expr dx(i)
KOKKOS_INLINE_FUNCTION GeneralFad(const GeneralFad &x)
Copy constructor.
KOKKOS_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
Forward-mode AD class templated on the storage for the derivative array.
KOKKOS_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.