Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
57namespace Sacado {
58
60 namespace CacheFad {
61
63
74 template <typename T, typename Storage>
75 class GeneralFad : public Storage {
76
77 public:
78
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
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
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
#define SACADO_INLINE_FUNCTION
expr expr dx(i)
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.
SACADO_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
SACADO_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
SACADO_INLINE_FUNCTION GeneralFad(const Storage &s)
Constructor with supplied storage s.
SACADO_INLINE_FUNCTION GeneralFad(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
SACADO_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
SACADO_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
SACADO_INLINE_FUNCTION GeneralFad(const GeneralFad &x)
Copy constructor.
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 GeneralFad()
Default constructor.
SACADO_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
SACADO_INLINE_FUNCTION GeneralFad(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
RemoveConst< T >::type value_type
Typename of values.
SACADO_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
SACADO_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
SACADO_INLINE_FUNCTION void cache() const
Cache values.
SACADO_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
SACADO_INLINE_FUNCTION ~GeneralFad()
Destructor.
ScalarType< value_type >::type scalar_type
Typename of scalar's (which may be different from T)
SACADO_INLINE_FUNCTION SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr< S > &x) const
Returns whether two Fad objects have the same values.
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.