Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
57namespace Sacado {
58
60 namespace Fad {
61
62#ifndef SACADO_FAD_DERIV_LOOP
63#if defined(SACADO_VIEW_CUDA_HIERARCHICAL_DFAD) && !defined(SACADO_DISABLE_CUDA_IN_KOKKOS) && defined(__CUDA_ARCH__)
64#define SACADO_FAD_DERIV_LOOP(I,SZ) for (int I=threadIdx.x; I<SZ; I+=blockDim.x)
65#else
66#define SACADO_FAD_DERIV_LOOP(I,SZ) for (int I=0; I<SZ; ++I)
67#endif
68#endif
69
70#ifndef SACADO_FAD_THREAD_SINGLE
71#if (defined(SACADO_VIEW_CUDA_HIERARCHICAL) || defined(SACADO_VIEW_CUDA_HIERARCHICAL_DFAD)) && !defined(SACADO_DISABLE_CUDA_IN_KOKKOS) && defined(__CUDA_ARCH__)
72#define SACADO_FAD_THREAD_SINGLE if (threadIdx.x == 0)
73#else
74#define SACADO_FAD_THREAD_SINGLE /* */
75#endif
76#endif
77
79
84 template <typename T, typename Storage>
85 class GeneralFad : public Storage {
86
87 public:
88
91
94
99
102 GeneralFad() : Storage(T(0.)) {}
103
105
108 template <typename S>
111 Storage(x) {}
112
114
118 GeneralFad(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) :
119 Storage(sz, x, zero_out) {}
120
122
128 GeneralFad(const int sz, const int i, const T & x) :
129 Storage(sz, x, InitDerivArray) {
130 this->fastAccessDx(i)=1.;
131 }
132
135 GeneralFad(const Storage& s) : Storage(s) {}
136
140 Storage(x) {}
141
143 template <typename S>
146 Storage(x.size(), T(0.), NoInitDerivArray)
147 {
148 const int sz = x.size();
149
150 if (sz) {
151 if (x.hasFastAccess())
153 this->fastAccessDx(i) = x.fastAccessDx(i);
154 else
156 this->fastAccessDx(i) = x.dx(i);
157 }
158
159 this->val() = x.val();
160 }
161
165
167
174 void diff(const int ith, const int n) {
175 if (this->size() != n)
176 this->resize(n);
177
178 this->zero();
179 this->fastAccessDx(ith) = T(1.);
180 }
181
184 void setUpdateValue(bool update_val) { }
185
188 bool updateValue() const { return true; }
189
192 void cache() const {}
193
195 template <typename S>
197 SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
198 typedef IsEqual<value_type> IE;
199 if (x.size() != this->size()) return false;
200 bool eq = IE::eval(x.val(), this->val());
201 const int sz = this->size();
203 eq = eq && IE::eval(x.dx(i), this->dx(i));
204 return eq;
205 }
206
208
213
219 int availableSize() const { return this->length(); }
220
223 bool hasFastAccess() const { return this->size()!=0; }
224
227 bool isPassive() const { return this->size()==0; }
228
231 void setIsConstant(bool is_const) {
232 if (is_const && this->size()!=0)
233 this->resize(0);
234 }
235
237
242
244 template <typename S>
246 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator=(const S& v) {
247 this->val() = v;
248 if (this->size()) this->resize(0);
249 return *this;
250 }
251
255 operator=(const GeneralFad& x) {
256 // Copy value & dx_
257 Storage::operator=(x);
258 return *this;
259 }
260
262 template <typename S>
264 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator=(const Expr<S>& x) {
265 const int xsz = x.size();
266
267 if (xsz != this->size())
268 this->resizeAndZero(xsz);
269
270 const int sz = this->size();
271
272 // For ViewStorage, the resize above may not in fact resize the
273 // derivative array, so it is possible that sz != xsz at this point.
274 // The only valid use case here is sz > xsz == 0, so we use sz in the
275 // assignment below
276
277 if (sz) {
278 if (x.hasFastAccess()) {
280 this->fastAccessDx(i) = x.fastAccessDx(i);
281 }
282 else
284 this->fastAccessDx(i) = x.dx(i);
285 }
286
287 this->val() = x.val();
288
289 return *this;
290 }
291
293
298
300 template <typename S>
302 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
303 this->val() += v;
304 return *this;
305 }
306
308 template <typename S>
310 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator -= (const S& v) {
311 this->val() -= v;
312 return *this;
313 }
314
316 template <typename S>
318 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator *= (const S& v) {
319 const int sz = this->size();
320 this->val() *= v;
322 this->fastAccessDx(i) *= v;
323 return *this;
324 }
325
327 template <typename S>
329 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
330 const int sz = this->size();
331 this->val() /= v;
333 this->fastAccessDx(i) /= v;
334 return *this;
335 }
336
339 GeneralFad& operator += (const GeneralFad& x) {
340 const int xsz = x.size(), sz = this->size();
341
342#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
343 if ((xsz != sz) && (xsz != 0) && (sz != 0))
344 throw "Fad Error: Attempt to assign with incompatible sizes";
345#endif
346
347 if (xsz) {
348 if (sz) {
350 this->fastAccessDx(i) += x.fastAccessDx(i);
351 }
352 else {
353 this->resizeAndZero(xsz);
355 this->fastAccessDx(i) = x.fastAccessDx(i);
356 }
357 }
358
359 this->val() += x.val();
360
361 return *this;
362 }
363
366 GeneralFad& operator -= (const GeneralFad& x) {
367 const int xsz = x.size(), sz = this->size();
368
369#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
370 if ((xsz != sz) && (xsz != 0) && (sz != 0))
371 throw "Fad Error: Attempt to assign with incompatible sizes";
372#endif
373
374 if (xsz) {
375 if (sz) {
377 this->fastAccessDx(i) -= x.fastAccessDx(i);
378 }
379 else {
380 this->resizeAndZero(xsz);
382 this->fastAccessDx(i) = -x.fastAccessDx(i);
383 }
384 }
385
386 this->val() -= x.val();
387
388
389 return *this;
390 }
391
394 GeneralFad& operator *= (const GeneralFad& x) {
395 const int xsz = x.size(), sz = this->size();
396 T xval = x.val();
397 T v = this->val();
398
399#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
400 if ((xsz != sz) && (xsz != 0) && (sz != 0))
401 throw "Fad Error: Attempt to assign with incompatible sizes";
402#endif
403
404 if (xsz) {
405 if (sz) {
407 this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
408 }
409 else {
410 this->resizeAndZero(xsz);
412 this->fastAccessDx(i) = v*x.fastAccessDx(i);
413 }
414 }
415 else {
416 if (sz) {
418 this->fastAccessDx(i) *= xval;
419 }
420 }
421
422 this->val() *= xval;
423
424 return *this;
425 }
426
429 GeneralFad& operator /= (const GeneralFad& x) {
430 const int xsz = x.size(), sz = this->size();
431 T xval = x.val();
432 T v = this->val();
433
434#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
435 if ((xsz != sz) && (xsz != 0) && (sz != 0))
436 throw "Fad Error: Attempt to assign with incompatible sizes";
437#endif
438
439 if (xsz) {
440 if (sz) {
442 this->fastAccessDx(i) =
443 ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
444 }
445 else {
446 this->resizeAndZero(xsz);
448 this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
449 }
450 }
451 else {
452 if (sz) {
454 this->fastAccessDx(i) /= xval;
455 }
456 }
457
458 this->val() /= xval;
459
460 return *this;
461 }
462
464 template <typename S>
466 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator += (const Expr<S>& x) {
467 const int xsz = x.size(), sz = this->size();
468
469#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
470 if ((xsz != sz) && (xsz != 0) && (sz != 0))
471 throw "Fad Error: Attempt to assign with incompatible sizes";
472#endif
473
474 if (xsz) {
475 if (sz) {
476 if (x.hasFastAccess())
478 this->fastAccessDx(i) += x.fastAccessDx(i);
479 else
481 this->fastAccessDx(i) += x.dx(i);
482 }
483 else {
484 this->resizeAndZero(xsz);
485 if (x.hasFastAccess())
487 this->fastAccessDx(i) = x.fastAccessDx(i);
488 else
490 this->fastAccessDx(i) = x.dx(i);
491 }
492 }
493
494 this->val() += x.val();
495
496 return *this;
497 }
498
500 template <typename S>
502 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
503 const int xsz = x.size(), sz = this->size();
504
505#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
506 if ((xsz != sz) && (xsz != 0) && (sz != 0))
507 throw "Fad Error: Attempt to assign with incompatible sizes";
508#endif
509
510 if (xsz) {
511 if (sz) {
512 if (x.hasFastAccess())
514 this->fastAccessDx(i) -= x.fastAccessDx(i);
515 else
517 this->fastAccessDx(i) -= x.dx(i);
518 }
519 else {
520 this->resizeAndZero(xsz);
521 if (x.hasFastAccess())
523 this->fastAccessDx(i) = -x.fastAccessDx(i);
524 else
526 this->fastAccessDx(i) = -x.dx(i);
527 }
528 }
529
530 this->val() -= x.val();
531
532
533 return *this;
534 }
535
537 template <typename S>
539 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
540 const int xsz = x.size(), sz = this->size();
541 T xval = x.val();
542 T v = this->val();
543
544#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
545 if ((xsz != sz) && (xsz != 0) && (sz != 0))
546 throw "Fad Error: Attempt to assign with incompatible sizes";
547#endif
548
549 if (xsz) {
550 if (sz) {
551 if (x.hasFastAccess())
553 this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
554 else
556 this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
557 }
558 else {
559 this->resizeAndZero(xsz);
560 if (x.hasFastAccess())
562 this->fastAccessDx(i) = v*x.fastAccessDx(i);
563 else
565 this->fastAccessDx(i) = v*x.dx(i);
566 }
567 }
568 else {
569 if (sz) {
571 this->fastAccessDx(i) *= xval;
572 }
573 }
574
575 this->val() *= xval;
576
577 return *this;
578 }
579
581 template <typename S>
583 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
584 const int xsz = x.size(), sz = this->size();
585 T xval = x.val();
586 T v = this->val();
587
588#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
589 if ((xsz != sz) && (xsz != 0) && (sz != 0))
590 throw "Fad Error: Attempt to assign with incompatible sizes";
591#endif
592
593 if (xsz) {
594 if (sz) {
595 if (x.hasFastAccess())
597 this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
598 else
600 this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
601 }
602 else {
603 this->resizeAndZero(xsz);
604 if (x.hasFastAccess())
606 this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
607 else
609 this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
610 }
611 }
612 else {
613 if (sz) {
615 this->fastAccessDx(i) /= xval;
616 }
617 }
618
619 this->val() /= xval;
620
621 return *this;
622 }
623
625
626 }; // class GeneralFad
627
628 } // namespace Fad
629
630} // namespace Sacado
631
632#endif // SACADO_FAD_GENERALFAD_HPP
#define SACADO_INLINE_FUNCTION
#define SACADO_FAD_DERIV_LOOP(I, SZ)
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
#define SACADO_ENABLE_VALUE_FUNC(RETURN_TYPE)
#define SACADO_ENABLE_VALUE_CTOR_DECL
#define SACADO_ENABLE_EXPR_FUNC(RETURN_TYPE)
#define SACADO_ENABLE_EXPR_CTOR_DECL
#define T
Definition: Sacado_rad.hpp:573
Wrapper for a generic expression template.
Forward-mode AD class templated on the storage for the derivative array.
RemoveConst< T >::type value_type
Typename of values.
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.
SACADO_INLINE_FUNCTION ~GeneralFad()
Destructor.
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 int availableSize() const
Returns number of derivative components that can be stored without reallocation.
SACADO_INLINE_FUNCTION GeneralFad()
Default constructor.
SACADO_INLINE_FUNCTION GeneralFad(const Storage &s)
Constructor with supplied storage s.
SACADO_INLINE_FUNCTION void cache() const
Cache values.
SACADO_INLINE_FUNCTION GeneralFad(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
SACADO_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
SACADO_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
SACADO_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
SACADO_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
SACADO_INLINE_FUNCTION GeneralFad(const GeneralFad &x)
Copy constructor.
SACADO_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
SACADO_INLINE_FUNCTION GeneralFad(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
ScalarType< value_type >::type scalar_type
Typename of scalar's (which may be different from T)
SACADO_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.
@ InitDerivArray
Initialize the derivative array.
@ NoInitDerivArray
Do not initialize the derivative array.
Base template specification for testing equivalence.