Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_CacheFad_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_CACHEFAD_GENERALFAD_HPP
53 #define SACADO_CACHEFAD_GENERALFAD_HPP
54 
56 
57 namespace Sacado {
58 
60  namespace CacheFad {
61 
63 
74  template <typename T, typename Storage>
75  class GeneralFad : public Storage {
76 
77  public:
78 
80  typedef typename RemoveConst<T>::type value_type;
81 
84 
89 
92  GeneralFad() : Storage(T(0.)) {}
93 
95 
98  template <typename S>
101  Storage(x) {}
102 
104 
108  GeneralFad(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) :
109  Storage(sz, x, zero_out) {}
110 
112 
118  GeneralFad(const int sz, const int i, const T & x) :
119  Storage(sz, x, InitDerivArray) {
120  this->fastAccessDx(i)=1.;
121  }
122 
125  GeneralFad(const Storage& s) : Storage(s) {}
126 
129  GeneralFad(const GeneralFad& x) :
130  Storage(x) {}
131 
133  template <typename S>
136  Storage(x.size(), T(0.), NoInitDerivArray) {
137  x.cache();
138 
139  const int sz = x.size();
140 
141  this->val() = x.val();
142 
143  if (sz) {
144  if (x.hasFastAccess())
145  for(int i=0; i<sz; ++i)
146  this->fastAccessDx(i) = x.fastAccessDx(i);
147  else
148  for(int i=0; i<sz; ++i)
149  this->fastAccessDx(i) = x.dx(i);
150  }
151  }
152 
156 
158 
165  void diff(const int ith, const int n) {
166  if (this->size() != n)
167  this->resize(n);
168 
169  this->zero();
170  this->fastAccessDx(ith) = T(1.);
171  }
172 
175  void setUpdateValue(bool update_val) { }
176 
179  bool updateValue() const { return true; }
180 
183  void cache() const {}
184 
186  template <typename S>
188  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
189  typedef IsEqual<value_type> IE;
190  if (x.size() != this->size()) return false;
191  bool eq = IE::eval(x.val(), this->val());
192  for (int i=0; i<this->size(); i++)
193  eq = eq && IE::eval(x.dx(i), this->dx(i));
194  return eq;
195  }
196 
198 
203 
209  int availableSize() const { return this->length(); }
210 
213  bool hasFastAccess() const { return this->size()!=0;}
214 
217  bool isPassive() const { return this->size()==0;}
218 
221  void setIsConstant(bool is_const) {
222  if (is_const && this->size()!=0)
223  this->resize(0);
224  }
225 
227 
232 
234  template <typename S>
236  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator=(const S& v) {
237  this->val() = v;
238  if (this->size()) this->resize(0);
239  return *this;
240  }
241 
244  GeneralFad&
245  operator=(const GeneralFad& x) {
246  // Copy val_ and dx_
247  Storage::operator=(x);
248  return *this;
249  }
250 
252  template <typename S>
254  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator=(const Expr<S>& x) {
255  x.cache();
256 
257  const int xsz = x.size();
258 
259  if (xsz != this->size())
260  this->resizeAndZero(xsz);
261 
262  const int sz = this->size();
263 
264  // For ViewStorage, the resize above may not in fact resize the
265  // derivative array, so it is possible that sz != xsz at this point.
266  // The only valid use case here is sz > xsz == 0, so we use sz in the
267  // assignment below
268 
269  if (sz) {
270  if (x.hasFastAccess())
271  for(int i=0; i<sz; ++i)
272  this->fastAccessDx(i) = x.fastAccessDx(i);
273  else
274  for(int i=0; i<sz; ++i)
275  this->fastAccessDx(i) = x.dx(i);
276  }
277 
278  this->val() = x.val();
279 
280  return *this;
281  }
282 
284 
289 
291  template <typename S>
293  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
294  this->val() += v;
295  return *this;
296  }
297 
299  template <typename S>
301  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator -= (const S& v) {
302  this->val() -= v;
303  return *this;
304  }
305 
307  template <typename S>
309  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator *= (const S& v) {
310  const int sz = this->size();
311  this->val() *= v;
312  for (int i=0; i<sz; ++i)
313  this->fastAccessDx(i) *= v;
314  return *this;
315  }
316 
318  template <typename S>
320  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
321  const int sz = this->size();
322  this->val() /= v;
323  for (int i=0; i<sz; ++i)
324  this->fastAccessDx(i) /= v;
325  return *this;
326  }
327 
330  GeneralFad& operator += (const GeneralFad& x) {
331  const int xsz = x.size(), sz = this->size();
332 
333 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
334  if ((xsz != sz) && (xsz != 0) && (sz != 0))
335  throw "Fad Error: Attempt to assign with incompatible sizes";
336 #endif
337 
338  if (xsz) {
339  if (sz) {
340  for (int i=0; i<sz; ++i)
341  this->fastAccessDx(i) += x.fastAccessDx(i);
342  }
343  else {
344  this->resizeAndZero(xsz);
345  for (int i=0; i<xsz; ++i)
346  this->fastAccessDx(i) = x.fastAccessDx(i);
347  }
348  }
349 
350  this->val() += x.val();
351 
352  return *this;
353  }
354 
357  GeneralFad& operator -= (const GeneralFad& x) {
358  const int xsz = x.size(), sz = this->size();
359 
360 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
361  if ((xsz != sz) && (xsz != 0) && (sz != 0))
362  throw "Fad Error: Attempt to assign with incompatible sizes";
363 #endif
364 
365  if (xsz) {
366  if (sz) {
367  for(int i=0; i<sz; ++i)
368  this->fastAccessDx(i) -= x.fastAccessDx(i);
369  }
370  else {
371  this->resizeAndZero(xsz);
372  for(int i=0; i<xsz; ++i)
373  this->fastAccessDx(i) = -x.fastAccessDx(i);
374  }
375  }
376 
377  this->val() -= x.val();
378 
379 
380  return *this;
381  }
382 
385  GeneralFad& operator *= (const GeneralFad& x) {
386  const int xsz = x.size(), sz = this->size();
387  T xval = x.val();
388  T v = this->val();
389 
390 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
391  if ((xsz != sz) && (xsz != 0) && (sz != 0))
392  throw "Fad Error: Attempt to assign with incompatible sizes";
393 #endif
394 
395  if (xsz) {
396  if (sz) {
397  for(int i=0; i<sz; ++i)
398  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
399  }
400  else {
401  this->resizeAndZero(xsz);
402  for(int i=0; i<xsz; ++i)
403  this->fastAccessDx(i) = v*x.fastAccessDx(i);
404  }
405  }
406  else {
407  if (sz) {
408  for (int i=0; i<sz; ++i)
409  this->fastAccessDx(i) *= xval;
410  }
411  }
412 
413  this->val() *= xval;
414 
415  return *this;
416  }
417 
420  GeneralFad& operator /= (const GeneralFad& x) {
421  const int xsz = x.size(), sz = this->size();
422  T xval = x.val();
423  T v = this->val();
424 
425 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
426  if ((xsz != sz) && (xsz != 0) && (sz != 0))
427  throw "Fad Error: Attempt to assign with incompatible sizes";
428 #endif
429 
430  if (xsz) {
431  if (sz) {
432  for(int i=0; i<sz; ++i)
433  this->fastAccessDx(i) =
434  ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
435  }
436  else {
437  this->resizeAndZero(xsz);
438  for(int i=0; i<xsz; ++i)
439  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
440  }
441  }
442  else {
443  if (sz) {
444  for (int i=0; i<sz; ++i)
445  this->fastAccessDx(i) /= xval;
446  }
447  }
448 
449  this->val() /= xval;
450 
451  return *this;
452  }
453 
455  template <typename S>
457  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator += (const Expr<S>& x) {
458  x.cache();
459 
460  const int xsz = x.size(), sz = this->size();
461 
462 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
463  if ((xsz != sz) && (xsz != 0) && (sz != 0))
464  throw "Fad Error: Attempt to assign with incompatible sizes";
465 #endif
466 
467  if (xsz) {
468  if (sz) {
469  if (x.hasFastAccess())
470  for (int i=0; i<sz; ++i)
471  this->fastAccessDx(i) += x.fastAccessDx(i);
472  else
473  for (int i=0; i<sz; ++i)
474  this->fastAccessDx(i) += x.dx(i);
475  }
476  else {
477  this->resizeAndZero(xsz);
478  if (x.hasFastAccess())
479  for (int i=0; i<xsz; ++i)
480  this->fastAccessDx(i) = x.fastAccessDx(i);
481  else
482  for (int i=0; i<xsz; ++i)
483  this->fastAccessDx(i) = x.dx(i);
484  }
485  }
486 
487  this->val() += x.val();
488 
489  return *this;
490  }
491 
493  template <typename S>
495  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
496  x.cache();
497 
498  const int xsz = x.size(), sz = this->size();
499 
500 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
501  if ((xsz != sz) && (xsz != 0) && (sz != 0))
502  throw "Fad Error: Attempt to assign with incompatible sizes";
503 #endif
504 
505  if (xsz) {
506  if (sz) {
507  if (x.hasFastAccess())
508  for(int i=0; i<sz; ++i)
509  this->fastAccessDx(i) -= x.fastAccessDx(i);
510  else
511  for (int i=0; i<sz; ++i)
512  this->fastAccessDx(i) -= x.dx(i);
513  }
514  else {
515  this->resizeAndZero(xsz);
516  if (x.hasFastAccess())
517  for(int i=0; i<xsz; ++i)
518  this->fastAccessDx(i) = -x.fastAccessDx(i);
519  else
520  for (int i=0; i<xsz; ++i)
521  this->fastAccessDx(i) = -x.dx(i);
522  }
523  }
524 
525  this->val() -= x.val();
526 
527 
528  return *this;
529  }
530 
532  template <typename S>
534  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
535  x.cache();
536 
537  const int xsz = x.size(), sz = this->size();
538  T xval = x.val();
539 
540 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
541  if ((xsz != sz) && (xsz != 0) && (sz != 0))
542  throw "Fad Error: Attempt to assign with incompatible sizes";
543 #endif
544 
545  if (xsz) {
546  if (sz) {
547  if (x.hasFastAccess())
548  for(int i=0; i<sz; ++i)
549  this->fastAccessDx(i) = this->val() * x.fastAccessDx(i) + this->fastAccessDx(i) * xval;
550  else
551  for (int i=0; i<sz; ++i)
552  this->fastAccessDx(i) = this->val() * x.dx(i) + this->fastAccessDx(i) * xval;
553  }
554  else {
555  this->resizeAndZero(xsz);
556  if (x.hasFastAccess())
557  for(int i=0; i<xsz; ++i)
558  this->fastAccessDx(i) = this->val() * x.fastAccessDx(i);
559  else
560  for (int i=0; i<xsz; ++i)
561  this->fastAccessDx(i) = this->val() * x.dx(i);
562  }
563  }
564  else {
565  if (sz) {
566  for (int i=0; i<sz; ++i)
567  this->fastAccessDx(i) *= xval;
568  }
569  }
570 
571  this->val() *= xval;
572 
573  return *this;
574  }
575 
577  template <typename S>
579  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
580  x.cache();
581 
582  const int xsz = x.size(), sz = this->size();
583  T xval = x.val();
584 
585 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
586  if ((xsz != sz) && (xsz != 0) && (sz != 0))
587  throw "Fad Error: Attempt to assign with incompatible sizes";
588 #endif
589 
590  if (xsz) {
591  if (sz) {
592  if (x.hasFastAccess())
593  for(int i=0; i<sz; ++i)
594  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - this->val()*x.fastAccessDx(i) )/ (xval*xval);
595  else
596  for (int i=0; i<sz; ++i)
597  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - this->val()*x.dx(i) )/ (xval*xval);
598  }
599  else {
600  this->resizeAndZero(xsz);
601  if (x.hasFastAccess())
602  for(int i=0; i<xsz; ++i)
603  this->fastAccessDx(i) = - this->val()*x.fastAccessDx(i) / (xval*xval);
604  else
605  for (int i=0; i<xsz; ++i)
606  this->fastAccessDx(i) = -this->val() * x.dx(i) / (xval*xval);
607  }
608  }
609  else {
610  if (sz) {
611  for (int i=0; i<sz; ++i)
612  this->fastAccessDx(i) /= xval;
613  }
614  }
615 
616  this->val() /= xval;
617 
618  return *this;
619  }
620 
622 
623  }; // class GeneralFad
624 
625  } // namespace CacheFad
626 
627 } // namespace Sacado
628 
629 #endif // SACADO_CACHEFAD_GENERALFAD_HPP
KOKKOS_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
KOKKOS_INLINE_FUNCTION GeneralFad(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
expr expr dx(i)
Forward-mode AD class templated on the storage for the derivative array.
#define SACADO_ENABLE_VALUE_CTOR_DECL
RemoveConst< T >::type value_type
Typename of values.
KOKKOS_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
KOKKOS_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
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(const Storage &s)
Constructor with supplied storage s.
#define SACADO_ENABLE_EXPR_CTOR_DECL
KOKKOS_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
KOKKOS_INLINE_FUNCTION GeneralFad(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
KOKKOS_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
KOKKOS_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
expr val()
#define KOKKOS_INLINE_FUNCTION
#define T
Definition: Sacado_rad.hpp:573
KOKKOS_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
#define SACADO_ENABLE_VALUE_FUNC(RETURN_TYPE)
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
Base template specification for testing equivalence.
KOKKOS_INLINE_FUNCTION ~GeneralFad()
Destructor.
Do not initialize the derivative array.
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors...
KOKKOS_INLINE_FUNCTION GeneralFad(const GeneralFad &x)
Copy constructor.
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
Initialize the derivative array.
KOKKOS_INLINE_FUNCTION void cache() const
Cache values.
KOKKOS_INLINE_FUNCTION GeneralFad(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
KOKKOS_INLINE_FUNCTION GeneralFad()
Default constructor.
KOKKOS_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
Wrapper for a generic expression template.