Stokhos 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_MP_Vector.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 // The forward-mode AD classes in Sacado are a derivative work of the
10 // expression template classes in the Fad package by Nicolas Di Cesare.
11 // The following banner is included in the original Fad source code:
12 //
13 // ************ DO NOT REMOVE THIS BANNER ****************
14 //
15 // Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
16 // http://www.ann.jussieu.fr/~dicesare
17 //
18 // CEMRACS 98 : C++ courses,
19 // templates : new C++ techniques
20 // for scientific computing
21 //
22 //********************************************************
23 //
24 // A short implementation ( not all operators and
25 // functions are overloaded ) of 1st order Automatic
26 // Differentiation in forward mode (FAD) using
27 // EXPRESSION TEMPLATES.
28 //
29 //********************************************************
30 // @HEADER
31 
32 #ifndef SACADO_FAD_GENERALFAD_MP_VECTOR_HPP
33 #define SACADO_FAD_GENERALFAD_MP_VECTOR_HPP
34 
35 #include "Sacado_Fad_GeneralFad.hpp"
37 
38 namespace Stokhos {
39  template <typename Ord, typename Val, int Num, typename Dev>
40  class StaticFixedStorage;
41 }
42 
43 namespace Sacado {
44 
45  namespace MP {
46  template <typename S> class Vector;
47  }
48 
49  namespace Fad {
50 
51  template <typename Ord, typename Val, int VecNum, typename Dev,
52  typename Storage>
53  struct ExprSpec< GeneralFad< Sacado::MP::Vector< Stokhos::StaticFixedStorage<Ord,Val,VecNum,Dev> >, Storage > > {
55  };
56 
58 
63  template <typename Ord, typename Val, int VecNum, typename Dev,
64  typename Storage>
65  class GeneralFad< Sacado::MP::Vector< Stokhos::StaticFixedStorage<Ord,Val,VecNum,Dev> >, Storage> : public Storage {
66 
67  public:
68 
70 
72  typedef typename RemoveConst<T>::type value_type;
73 
75  typedef typename ScalarType<value_type>::type scalar_type;
76 
77  typedef typename value_type::value_type val_type;
78 
83 
85  KOKKOS_INLINE_FUNCTION
86  GeneralFad() : Storage(T(0.)) {}
87 
89 
92  template <typename S>
93  KOKKOS_INLINE_FUNCTION
94  GeneralFad(const S& x, SACADO_ENABLE_VALUE_CTOR_DECL) :
95  Storage(x) {}
96 
98 
101  KOKKOS_INLINE_FUNCTION
102  GeneralFad(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) :
103  Storage(sz, x, zero_out) {}
104 
106 
111  KOKKOS_INLINE_FUNCTION
112  GeneralFad(const int sz, const int i, const T & x) :
113  Storage(sz, x) {
114  this->fastAccessDx(i)=1.;
115  }
116 
118  KOKKOS_INLINE_FUNCTION
119  GeneralFad(const Storage& s) : Storage(s) {}
120 
122  KOKKOS_INLINE_FUNCTION
123  GeneralFad(const GeneralFad& x) :
124  Storage(x) {}
125 
127  template <typename S>
128  KOKKOS_INLINE_FUNCTION
129  GeneralFad(const Expr<S>& x, SACADO_ENABLE_EXPR_CTOR_DECL) :
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  for (int j=0; j<VecNum; ++j)
138  fastAccessDx(i,j) = x.fastAccessDx(i,j);
139  else
140  for(int i=0; i<sz; ++i)
141  for (int j=0; j<VecNum; ++j)
142  fastAccessDx(i,j) = x.dx(i,j);
143  }
144 
145  for (int j=0; j<VecNum; ++j)
146  val(j) = x.val(j);
147  }
148 
150  KOKKOS_INLINE_FUNCTION
152 
154 
160  KOKKOS_INLINE_FUNCTION
161  void diff(const int ith, const int n) {
162  if (this->size() != n)
163  this->resize(n);
164 
165  this->zero();
166  this->fastAccessDx(ith) = T(1.);
167  }
168 
170  KOKKOS_INLINE_FUNCTION
171  void setUpdateValue(bool update_val) { }
172 
174  KOKKOS_INLINE_FUNCTION
175  bool updateValue() const { return true; }
176 
178  template <typename S>
179  KOKKOS_INLINE_FUNCTION
180  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
181  typedef IsEqual<value_type> IE;
182  if (x.size() != this->size()) return false;
183  bool eq = IE::eval(x.val(), this->val());
184  for (int i=0; i<this->size(); i++)
185  eq = eq && IE::eval(x.dx(i), this->dx(i));
186  return eq;
187  }
188 
190 
195 
197  KOKKOS_INLINE_FUNCTION
198  const T& val() const { return Storage::val();}
199 
201  KOKKOS_INLINE_FUNCTION
202  T& val() { return Storage::val();}
203 
205  KOKKOS_INLINE_FUNCTION
206  const val_type& val(int j) const { return Storage::val().fastAccessCoeff(j);}
207 
209  KOKKOS_INLINE_FUNCTION
210  val_type& val(int j) { return Storage::val().fastAccessCoeff(j);}
211 
213 
218 
223  KOKKOS_INLINE_FUNCTION
224  int availableSize() const { return this->length(); }
225 
227  KOKKOS_INLINE_FUNCTION
228  bool hasFastAccess() const { return this->size()!=0; }
229 
231  KOKKOS_INLINE_FUNCTION
232  bool isPassive() const { return this->size()==0; }
233 
235  KOKKOS_INLINE_FUNCTION
236  void setIsConstant(bool is_const) {
237  if (is_const && this->size()!=0)
238  this->resize(0);
239  }
240 
242  KOKKOS_INLINE_FUNCTION
243  const T* dx() const { return this->dx_; }
244 
246  KOKKOS_INLINE_FUNCTION
247  T dx(int i) const { return this->size() ? this->dx_[i] : T(0.); }
248 
250  KOKKOS_INLINE_FUNCTION
251  T& fastAccessDx(int i) { return this->dx_[i];}
252 
254  KOKKOS_INLINE_FUNCTION
255  const T& fastAccessDx(int i) const { return this->dx_[i];}
256 
258  KOKKOS_INLINE_FUNCTION
259  val_type dx(int i, int j) const { return this->size() ? this->dx_[i].fastAccessCoeff(j) : val_type(0.0); }
260 
262  KOKKOS_INLINE_FUNCTION
263  val_type& fastAccessDx(int i, int j) {
264  return this->dx_[i].fastAccessCoeff(j);
265  }
266 
268  KOKKOS_INLINE_FUNCTION
269  const val_type& fastAccessDx(int i, int j) const {
270  return this->dx_[i].fastAccessCoeff(j);
271  }
272 
274 
279 
281  template <typename S>
282  KOKKOS_INLINE_FUNCTION
283  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator=(const S& v) {
284  this->val() = v;
285  if (this->size()) this->resize(0);
286  return *this;
287  }
288 
290  KOKKOS_INLINE_FUNCTION
291  GeneralFad&
292  operator=(const GeneralFad& x) {
293  // Copy value & dx_
294  Storage::operator=(x);
295  return *this;
296  }
297 
299  template <typename S>
300  KOKKOS_INLINE_FUNCTION
301  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator=(const Expr<S>& x) {
302  const int xsz = x.size();
303 
304  if (xsz != this->size())
305  this->resizeAndZero(xsz);
306 
307  const int sz = this->size();
308 
309  // For ViewStorage, the resize above may not in fact resize the
310  // derivative array, so it is possible that sz != xsz at this point.
311  // The only valid use case here is sz > xsz == 0, so we use sz in the
312  // assignment below
313 
314  if (sz) {
315  if (x.hasFastAccess())
316  for(int i=0; i<sz; ++i)
317  for (int j=0; j<VecNum; ++j)
318  fastAccessDx(i,j) = x.fastAccessDx(i,j);
319  else
320  for(int i=0; i<sz; ++i)
321  for (int j=0; j<VecNum; ++j)
322  fastAccessDx(i,j) = x.dx(i,j);
323  }
324 
325  for (int j=0; j<VecNum; ++j)
326  val(j) = x.val(j);
327 
328  return *this;
329  }
330 
332 
337 
339  template <typename S>
340  KOKKOS_INLINE_FUNCTION
341  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
342  this->val() += v;
343  return *this;
344  }
345 
347  template <typename S>
348  KOKKOS_INLINE_FUNCTION
349  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator -= (const S& v) {
350  this->val() -= v;
351  return *this;
352  }
353 
355  template <typename S>
356  KOKKOS_INLINE_FUNCTION
357  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator *= (const S& v) {
358  const int sz = this->size();
359  this->val() *= v;
360  for (int i=0; i<sz; ++i)
361  this->fastAccessDx(i) *= v;
362  return *this;
363  }
364 
366  template <typename S>
367  KOKKOS_INLINE_FUNCTION
368  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
369  const int sz = this->size();
370  this->val() /= v;
371  for (int i=0; i<sz; ++i)
372  this->fastAccessDx(i) /= v;
373  return *this;
374  }
375 
377  KOKKOS_INLINE_FUNCTION
378  GeneralFad& operator += (const GeneralFad& x) {
379  const int xsz = x.size(), sz = this->size();
380 
381 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
382  if ((xsz != sz) && (xsz != 0) && (sz != 0))
383  throw "Fad Error: Attempt to assign with incompatible sizes";
384 #endif
385 
386  if (xsz) {
387  if (sz) {
388  for (int i=0; i<sz; ++i)
389  this->fastAccessDx(i) += x.fastAccessDx(i);
390  }
391  else {
392  this->resizeAndZero(xsz);
393  for (int i=0; i<xsz; ++i)
394  this->fastAccessDx(i) = x.fastAccessDx(i);
395  }
396  }
397 
398  this->val() += x.val();
399  return *this;
400  }
401 
403  KOKKOS_INLINE_FUNCTION
404  GeneralFad& operator -= (const GeneralFad& x) {
405  const int xsz = x.size(), sz = this->size();
406 
407 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
408  if ((xsz != sz) && (xsz != 0) && (sz != 0))
409  throw "Fad Error: Attempt to assign with incompatible sizes";
410 #endif
411 
412  if (xsz) {
413  if (sz) {
414  for(int i=0; i<sz; ++i)
415  this->fastAccessDx(i) -= x.fastAccessDx(i);
416  }
417  else {
418  this->resizeAndZero(xsz);
419  for(int i=0; i<xsz; ++i)
420  this->fastAccessDx(i) = -x.fastAccessDx(i);
421  }
422  }
423 
424  this->val() -= x.val();
425 
426  return *this;
427  }
428 
430  KOKKOS_INLINE_FUNCTION
431  GeneralFad& operator *= (const GeneralFad& x) {
432  const int xsz = x.size(), sz = this->size();
433  T xval = x.val();
434  T v = this->val();
435 
436 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
437  if ((xsz != sz) && (xsz != 0) && (sz != 0))
438  throw "Fad Error: Attempt to assign with incompatible sizes";
439 #endif
440 
441  if (xsz) {
442  if (sz) {
443  for(int i=0; i<sz; ++i)
444  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
445  }
446  else {
447  this->resizeAndZero(xsz);
448  for(int i=0; i<xsz; ++i)
449  this->fastAccessDx(i) = v*x.fastAccessDx(i);
450  }
451  }
452  else {
453  if (sz) {
454  for (int i=0; i<sz; ++i)
455  this->fastAccessDx(i) *= xval;
456  }
457  }
458 
459  this->val() *= xval;
460 
461  return *this;
462  }
463 
465  KOKKOS_INLINE_FUNCTION
466  GeneralFad& operator /= (const GeneralFad& x) {
467  const int xsz = x.size(), sz = this->size();
468  T xval = x.val();
469  T v = this->val();
470 
471 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
472  if ((xsz != sz) && (xsz != 0) && (sz != 0))
473  throw "Fad Error: Attempt to assign with incompatible sizes";
474 #endif
475 
476  if (xsz) {
477  if (sz) {
478  for(int i=0; i<sz; ++i)
479  this->fastAccessDx(i) =
480  ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
481  }
482  else {
483  this->resizeAndZero(xsz);
484  for(int i=0; i<xsz; ++i)
485  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
486  }
487  }
488  else {
489  if (sz) {
490  for (int i=0; i<sz; ++i)
491  this->fastAccessDx(i) /= xval;
492  }
493  }
494 
495  this->val() /= xval;
496 
497  return *this;
498  }
499 
501  template <typename S>
502  KOKKOS_INLINE_FUNCTION
503  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator += (const Expr<S>& x) {
504  const int xsz = x.size();
505  int sz = this->size();
506 
507 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
508  if ((xsz != sz) && (xsz != 0) && (sz != 0))
509  throw "Fad Error: Attempt to assign with incompatible sizes";
510 #endif
511 
512  if (xsz > sz) {
513  this->resizeAndZero(xsz);
514  sz = this->size();
515  }
516 
517  if (sz) {
518  if (x.hasFastAccess())
519  for(int i=0; i<sz; ++i)
520  for (int j=0; j<VecNum; ++j)
521  fastAccessDx(i,j) += x.fastAccessDx(i,j);
522  else
523  for(int i=0; i<sz; ++i)
524  for (int j=0; j<VecNum; ++j)
525  fastAccessDx(i,j) += x.dx(i,j);
526  }
527 
528  for (int j=0; j<VecNum; ++j)
529  val(j) += x.val(j);
530 
531  return *this;
532  }
533 
535  template <typename S>
536  KOKKOS_INLINE_FUNCTION
537  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
538  const int xsz = x.size(), sz = this->size();
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  for (int j=0; j<VecNum; ++j)
550  fastAccessDx(i,j) -= x.fastAccessDx(i,j);
551  else
552  for (int i=0; i<sz; ++i)
553  for (int j=0; j<VecNum; ++j)
554  fastAccessDx(i,j) -= x.dx(i,j);
555  }
556  else {
557  this->resizeAndZero(xsz);
558  if (x.hasFastAccess())
559  for(int i=0; i<xsz; ++i)
560  for (int j=0; j<VecNum; ++j)
561  fastAccessDx(i,j) = -x.fastAccessDx(i,j);
562  else
563  for (int i=0; i<xsz; ++i)
564  for (int j=0; j<VecNum; ++j)
565  fastAccessDx(i,j) = -x.dx(i,j);
566  }
567  }
568 
569  for (int j=0; j<VecNum; ++j)
570  val(j) -= x.val(j);
571 
572  return *this;
573  }
574 
576  template <typename S>
577  KOKKOS_INLINE_FUNCTION
578  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
579  const int xsz = x.size(), sz = this->size();
580  T xval = x.val();
581  T v = this->val();
582 
583 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
584  if ((xsz != sz) && (xsz != 0) && (sz != 0))
585  throw "Fad Error: Attempt to assign with incompatible sizes";
586 #endif
587 
588  if (xsz) {
589  if (sz) {
590  if (x.hasFastAccess())
591  for(int i=0; i<sz; ++i)
592  for (int j=0; j<VecNum; ++j)
593  fastAccessDx(i,j) = v.fastAccessCoeff(j)*x.fastAccessDx(i,j) + fastAccessDx(i,j)*xval.fastAccessCoeff(j);
594  else
595  for (int i=0; i<sz; ++i)
596  for (int j=0; j<VecNum; ++j)
597  fastAccessDx(i,j) = v.fastAccessCoeff(j)*x.dx(i,j) + fastAccessDx(i,j)*xval.fastAccessCoeff(j);
598  }
599  else {
600  this->resizeAndZero(xsz);
601  if (x.hasFastAccess())
602  for(int i=0; i<xsz; ++i)
603  for (int j=0; j<VecNum; ++j)
604  fastAccessDx(i,j) = v.fastAccessCoeff(j)*x.fastAccessDx(i,j);
605  else
606  for (int i=0; i<xsz; ++i)
607  for (int j=0; j<VecNum; ++j)
608  fastAccessDx(i,j) = v.fastAccessCoeff(j)*x.dx(i,j);
609  }
610  }
611  else {
612  if (sz) {
613  for (int i=0; i<sz; ++i)
614  for (int j=0; j<VecNum; ++j)
615  fastAccessDx(i,j) *= xval.fastAccessCoeff(j);
616  }
617  }
618 
619  for (int j=0; j<VecNum; ++j)
620  val(j) *= xval.fastAccessCoeff(j);
621 
622  return *this;
623  }
624 
626  template <typename S>
627  KOKKOS_INLINE_FUNCTION
628  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
629  const int xsz = x.size(), sz = this->size();
630  T xval = x.val();
631  T v = this->val();
632 
633 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
634  if ((xsz != sz) && (xsz != 0) && (sz != 0))
635  throw "Fad Error: Attempt to assign with incompatible sizes";
636 #endif
637 
638  if (xsz) {
639  T xval2 = xval*xval;
640  if (sz) {
641  if (x.hasFastAccess())
642  for(int i=0; i<sz; ++i)
643  for (int j=0; j<VecNum; ++j)
644  fastAccessDx(i,j) = ( fastAccessDx(i,j)*xval.fastAccessCoeff(j) - v.fastAccessCoeff(j)*x.fastAccessDx(i,j) )/ (xval2.fastAccessCoeff(j));
645  else
646  for (int i=0; i<sz; ++i)
647  for (int j=0; j<VecNum; ++j)
648  fastAccessDx(i,j) = ( fastAccessDx(i,j)*xval.fastAccessCoeff(j) - v.fastAccessCoeff(j)*x.dx(i,j) )/ (xval2.fastAccessCoeff(j));
649  }
650  else {
651  this->resizeAndZero(xsz);
652  if (x.hasFastAccess())
653  for(int i=0; i<xsz; ++i)
654  for (int j=0; j<VecNum; ++j)
655  fastAccessDx(i,j) = - v.fastAccessCoeff(j)*x.fastAccessDx(i,j) / (xval2.fastAccessCoeff(j));
656  else
657  for (int i=0; i<xsz; ++i)
658  for (int j=0; j<VecNum; ++j)
659  fastAccessDx(i,j) = -v.fastAccessCoeff(j)*x.dx(i,j) / (xval2.fastAccessCoeff(j));
660  }
661  }
662  else {
663  if (sz) {
664  for (int i=0; i<sz; ++i)
665  for (int j=0; j<VecNum; ++j)
666  this->fastAccessDx(i,j) /= xval.fastAccessCoeff(j);
667  }
668  }
669 
670  for (int j=0; j<VecNum; ++j)
671  val(j) /= xval.fastAccessCoeff(j);
672 
673  return *this;
674  }
675 
677 
678  private:
679 
680  template <typename S>
681  KOKKOS_INLINE_FUNCTION
682  void
683  fastAssign(const Expr<S>& x) {
684  const int sz = this->size();
685  for(int i=0; i<sz; ++i)
686  for (int j=0; j<VecNum; ++j)
687  fastAccessDx(i,j) = x.fastAccessDx(i,j);
688 
689  for (int j=0; j<VecNum; ++j)
690  val(j) = x.val(j);
691  }
692 
693  template <typename S>
694  KOKKOS_INLINE_FUNCTION
695  void
696  slowAssign(const Expr<S>& x) {
697  const int sz = this->size();
698  for(int i=0; i<sz; ++i)
699  for (int j=0; j<VecNum; ++j)
700  fastAccessDx(i,j) = x.dx(i,j);
701 
702  for (int j=0; j<VecNum; ++j)
703  val(j) = x.val(j);
704  }
705 
706  }; // class GeneralFad
707 
708  } // namespace Fad
709 
710 } // namespace Sacado
711 
712 #endif // SACADO_FAD_GENERALFAD_HPP
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 Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
KOKKOS_INLINE_FUNCTION T dx(int i) const
Returns derivative component i with bounds checking.
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
KOKKOS_INLINE_FUNCTION val_type & fastAccessDx(int i, int j)
Returns derivative component i without bounds checking.
KOKKOS_INLINE_FUNCTION const val_type & fastAccessDx(int i, int j) const
Returns derivative component i without bounds checking.
KOKKOS_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
KOKKOS_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
expr expr expr expr fastAccessDx(i, j)) FAD_UNARYOP_MACRO(exp
KOKKOS_INLINE_FUNCTION GeneralFad(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
KOKKOS_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
expr expr expr dx(i, j)
expr val()
KOKKOS_INLINE_FUNCTION T & fastAccessDx(int i)
Returns derivative component i without bounds checking.
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 val_type dx(int i, int j) const
Returns derivative component i with bounds checking.
KOKKOS_INLINE_FUNCTION GeneralFad(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
KOKKOS_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
KOKKOS_INLINE_FUNCTION const T & fastAccessDx(int i) const
Returns derivative component i without bounds checking.
ordinal_type size() const
Return size.