Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_ELRFad_SFad.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_ELRFAD_SFAD_HPP
53#define SACADO_ELRFAD_SFAD_HPP
54
60
61namespace Sacado {
62
63 namespace ELRFad {
64
66 template <typename T, int Num>
67 struct SFadExprTag {};
68
69 // Forward declaration
70 template <typename T, int Num> class SFad;
71
79 template <typename T, int Num>
80 class Expr< SFadExprTag<T,Num> > {
81
82 public:
83
86
89
92
94 static const int num_args = 1;
95
97 static const bool is_linear = true;
98
103
106 Expr() : val_( T(0.)) { ss_array<T>::zero(dx_, Num); }
107
109
112 template <typename S>
115 val_(x) {
116 ss_array<T>::zero(dx_, Num);
117 }
118
120
124 Expr(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) : val_(x) {
125#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
126 if (sz != Num)
127 throw "SELRFad::SFad() Error: Supplied derivative dimension does not match compile time length.";
128#endif
129
130 if (zero_out == InitDerivArray)
131 ss_array<T>::zero(dx_, Num);
132 }
133
135
141 Expr(const int sz, const int i, const T & x) :
142 val_(x) {
143#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
144 if (sz != Num)
145 throw "SELRFad::SFad() Error: Supplied derivative dimension does not match compile time length.";
146 if (i >= Num)
147 throw "SELRFad::SFad() Error: Invalid derivative index.";
148#endif
149
150 ss_array<T>::zero(dx_, Num);
151 dx_[i]=1.;
152 }
153
156 Expr(const Expr& x) : val_(x.val_) {
157 //ss_array<T>::copy(x.dx_, dx_, Num);
158 for (int i=0; i<Num; i++)
159 dx_[i] = x.dx_[i];
160 }
161
163 template <typename S>
166#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
167 if (x.size() != Num)
168 throw "SELRFad::SFad() Error: Attempt to assign with incompatible sizes";
169#endif
170
171 // Compute partials
172 LocalAccumOp< Expr<S> > op(x);
173
174 // Compute each tangent direction
175 for(int i=0; i<Num; ++i) {
176 op.t = T(0.);
177 op.i = i;
178
179 // Automatically unrolled loop that computes
180 // for (int j=0; j<N; j++)
181 // op.t += op.partials[j] * x.getTangent<j>(i);
183
184 dx_[i] = op.t;
185
186 }
187
188 this->val() = x.val();
189 }
190
193 ~Expr() {}
194
196
203 void diff(const int ith, const int n) {
204#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
205 if (n != Num)
206 throw "SELRFad::diff() Error: Supplied derivative dimension does not match compile time length.";
207#endif
208
209 ss_array<T>::zero(dx_, Num);
210 dx_[ith] = T(1.);
211 }
212
214
219 void resize(int sz) {
220#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
221 if (sz != Num)
222 throw "SELRFad::resize() Error: Cannot resize fixed derivative array dimension";
223#endif
224 }
225
227
232 void expand(int sz) { resize(sz); }
233
236 void zero() { ss_array<T>::zero(dx_, Num); }
237
240 void setUpdateValue(bool update_val) { }
241
244 bool updateValue() const { return true; }
245
248 void cache() const {}
249
251 template <typename S>
253 SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
254 typedef IsEqual<value_type> IE;
255 if (x.size() != this->size()) return false;
256 bool eq = IE::eval(x.val(), this->val());
257 for (int i=0; i<this->size(); i++)
258 eq = eq && IE::eval(x.dx(i), this->dx(i));
259 return eq;
260 }
261
263
268
271 const T& val() const { return val_;}
272
275 T& val() { return val_;}
276
278
283
286 int size() const { return Num;}
287
293 int availableSize() const { return Num; }
294
297 bool hasFastAccess() const { return true; }
298
301 bool isPassive() const { return false; }
302
305 void setIsConstant(bool is_const) {}
306
309 const T* dx() const { return &(dx_[0]);}
310
313 const T& dx(int i) const { return dx_[i]; }
314
317 T& fastAccessDx(int i) { return dx_[i];}
318
321 const T& fastAccessDx(int i) const { return dx_[i];}
322
325 void computePartials(const T& bar, value_type partials[]) const {
326 partials[0] = bar;
327 }
328
331 void getTangents(int i, value_type dots[]) const {
332 dots[0] = this->dx_[i];
333 }
334
336 template <int Arg>
338 bool isActive() const { return true; }
339
341 template <int Arg>
343 const T& getTangent(int i) const { return this->dx_[i]; }
344
347 const value_type* getDx(int j) const { return this->dx(); }
348
350
355
357 template <typename S>
359 SACADO_ENABLE_VALUE_FUNC(Expr&) operator=(const S& v) {
360 val_ = v;
361 ss_array<T>::zero(dx_, Num);
362 return *this;
363 }
364
367 Expr& operator=(const Expr& x) {
368 if (this != &x) {
369 // Copy value
370 val_ = x.val_;
371
372 // Copy dx_
373 //ss_array<T>::copy(x.dx_, dx_, Num);
374 for (int i=0; i<Num; i++)
375 dx_[i] = x.dx_[i];
376 }
377 return *this;
378 }
379
381 template <typename S>
383 SACADO_ENABLE_EXPR_FUNC(Expr&) operator=(const Expr<S>& x) {
384#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
385 if (x.size() != Num)
386 throw "SELRFad::operator=() Error: Attempt to assign with incompatible sizes";
387#endif
388
389 // Compute partials
390 LocalAccumOp< Expr<S> > op(x);
391
392 // Compute each tangent direction
393 for(int i=0; i<Num; ++i) {
394 op.t = T(0.);
395 op.i = i;
396
397 // Automatically unrolled loop that computes
398 // for (int j=0; j<N; j++)
399 // op.t += op.partials[j] * x.getTangent<j>(i);
401
402 dx_[i] = op.t;
403 }
404
405 // Compute value
406 val_ = x.val();
407
408 return *this;
409 }
410
412
417
419 template <typename S>
421 SACADO_ENABLE_VALUE_FUNC(Expr&) operator += (const S& v) {
422 this->val() += v;
423 return *this;
424 }
425
427 template <typename S>
429 SACADO_ENABLE_VALUE_FUNC(Expr&) operator -= (const S& v) {
430 this->val() -= v;
431 return *this;
432 }
433
435 template <typename S>
437 SACADO_ENABLE_VALUE_FUNC(Expr&) operator *= (const S& v) {
438 this->val() *= v;
439 for (int i=0; i<Num; ++i)
440 dx_[i] *= v;
441 return *this;
442 }
443
445 template <typename S>
447 SACADO_ENABLE_VALUE_FUNC(Expr&) operator /= (const S& v) {
448 this->val() /= v;
449 for (int i=0; i<Num; ++i)
450 dx_[i] /= v;
451 return *this;
452 }
453
455 template <typename S>
457 SACADO_ENABLE_EXPR_FUNC(Expr&) operator += (const Expr<S>& x) {
458#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
459 if (x.size() != Num)
460 throw "SELRFad::operator+=() Error: Attempt to assign with incompatible sizes";
461#endif
462
463 // Compute partials
464 LocalAccumOp< Expr<S> > op(x);
465
466 // Compute each tangent direction
467 for(int i=0; i<Num; ++i) {
468 op.t = T(0.);
469 op.i = i;
470
471 // Automatically unrolled loop that computes
472 // for (int j=0; j<N; j++)
473 // op.t += op.partials[j] * x.getTangent<j>(i);
475
476 dx_[i] += op.t;
477 }
478
479 // Compute value
480 val_ += x.val();
481
482 return *this;
483 }
484
486 template <typename S>
488 SACADO_ENABLE_EXPR_FUNC(Expr&) operator -= (const Expr<S>& x) {
489#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
490 if (x.size() != Num)
491 throw "SELRFad::operator-=() Error: Attempt to assign with incompatible sizes";
492#endif
493
494 // Compute partials
495 LocalAccumOp< Expr<S> > op(x);
496
497 // Compute each tangent direction
498 for(int i=0; i<Num; ++i) {
499 op.t = T(0.);
500 op.i = i;
501
502 // Automatically unrolled loop that computes
503 // for (int j=0; j<N; j++)
504 // op.t += op.partials[j] * x.getTangent<j>(i);
506
507 dx_[i] -= op.t;
508 }
509
510 // Compute value
511 val_ -= x.val();
512
513 return *this;
514 }
515
517 template <typename S>
519 SACADO_ENABLE_EXPR_FUNC(Expr&) operator *= (const Expr<S>& x) {
520 T xval = x.val();
521
522#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
523 if (x.size() != Num)
524 throw "SELRFad::operator*=() Error: Attempt to assign with incompatible sizes";
525#endif
526
527 // Compute partials
528 LocalAccumOp< Expr<S> > op(x);
529
530 // Compute each tangent direction
531 for(int i=0; i<Num; ++i) {
532 op.t = T(0.);
533 op.i = i;
534
535 // Automatically unrolled loop that computes
536 // for (int j=0; j<N; j++)
537 // op.t += op.partials[j] * x.getTangent<j>(i);
539
540 dx_[i] = val_ * op.t + dx_[i] * xval;
541 }
542
543 // Compute value
544 val_ *= xval;
545
546 return *this;
547 }
548
550 template <typename S>
552 SACADO_ENABLE_EXPR_FUNC(Expr&) operator /= (const Expr<S>& x) {
553 T xval = x.val();
554
555#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
556 if (x.size() != Num)
557 throw "SELRFad::operator/=() Error: Attempt to assign with incompatible sizes";
558#endif
559
560 // Compute partials
561 LocalAccumOp< Expr<S> > op(x);
562
563 T xval2 = xval*xval;
564
565 // Compute each tangent direction
566 for(int i=0; i<Num; ++i) {
567 op.t = T(0.);
568 op.i = i;
569
570 // Automatically unrolled loop that computes
571 // for (int j=0; j<N; j++)
572 // op.t += op.partials[j] * x.getTangent<j>(i);
574
575 dx_[i] = (dx_[i] * xval - val_ * op.t) / xval2;
576 }
577
578 // Compute value
579 val_ /= xval;
580
581 return *this;
582 }
583
585
586 protected:
587
589 T dx_[Num];
590
593
594 // Functor for mpl::for_each to compute the local accumulation
595 // of a tangent derivative
596 //
597 // We use getTangent<>() to get dx components from expression
598 // arguments instead of getting the argument directly or extracting
599 // the dx array due to striding in ViewFad (or could use striding
600 // directly here if we need to get dx array).
601 template <typename ExprT>
602 struct LocalAccumOp {
603 typedef typename ExprT::value_type value_type;
604 static const int N = ExprT::num_args;
605 const ExprT& x;
606 mutable value_type t;
607 value_type partials[N];
608 int i;
610 LocalAccumOp(const ExprT& x_) : x(x_) {
611 x.computePartials(value_type(1.), partials);
612 }
614 void getTangents(int i_) { i = i_; }
615 template <typename ArgT>
617 void operator () (ArgT arg) const {
618 const int Arg = ArgT::value;
619 if (x.template isActive<Arg>())
620 t += partials[Arg] * x.template getTangent<Arg>(i);
621 }
622 };
623
624 }; // class Expr<SFadExprTag>
625
626 } // namespace ELRFad
627
628} // namespace Sacado
629
630#define FAD_NS ELRFad
632#undef FAD_NS
633
635#include "Sacado_ELRFad_Ops.hpp"
636
637#endif // SACADO_ELRFAD_SFAD_HPP
#define SACADO_INLINE_FUNCTION
expr expr dx(i)
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
const int N
SACADO_INLINE_FUNCTION Expr(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
SACADO_INLINE_FUNCTION const T & fastAccessDx(int i) const
Returns derivative component i without bounds checking.
SFad< value_type, Num > base_expr_type
Typename of base-expressions.
SACADO_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
SACADO_INLINE_FUNCTION int size() const
Returns number of derivative components.
SACADO_INLINE_FUNCTION const T * dx() const
Returns derivative array.
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 Expr()
Default constructor.
SACADO_INLINE_FUNCTION const T & dx(int i) const
Returns derivative component i with bounds checking.
SACADO_INLINE_FUNCTION void cache() const
Cache values.
SACADO_INLINE_FUNCTION Expr(const Expr &x)
Copy constructor.
SACADO_INLINE_FUNCTION void zero()
Zero out the derivative array.
SACADO_INLINE_FUNCTION T & fastAccessDx(int i)
Returns derivative component i without bounds checking.
SACADO_INLINE_FUNCTION ~Expr()
Destructor.
SACADO_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
SACADO_INLINE_FUNCTION Expr(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
SACADO_INLINE_FUNCTION void getTangents(int i, value_type dots[]) const
Return tangent component i of arguments.
SACADO_INLINE_FUNCTION void expand(int sz)
Expand derivative array to size sz.
ScalarType< value_type >::type scalar_type
Typename of scalar's (which may be different from T)
SACADO_INLINE_FUNCTION const T & getTangent(int i) const
Return tangent component i of argument Arg.
SACADO_INLINE_FUNCTION SACADO_ENABLE_VALUE_FUNC(Expr &) operator
Assignment operator with constant right-hand-side.
SACADO_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
SACADO_INLINE_FUNCTION void resize(int sz)
Resize derivative array to length sz.
SACADO_INLINE_FUNCTION void diff(const int ith, const int n)
Set Fad object as the ith independent variable.
SACADO_INLINE_FUNCTION const T & val() const
Returns value.
SACADO_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
SACADO_INLINE_FUNCTION const value_type * getDx(int j) const
Get dx array.
SACADO_INLINE_FUNCTION T & val()
Returns value.
SACADO_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
SACADO_INLINE_FUNCTION void computePartials(const T &bar, value_type partials[]) const
Return partials w.r.t. arguments.
SACADO_INLINE_FUNCTION Expr(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
SACADO_INLINE_FUNCTION Expr(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
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 bool isActive() const
Return whether argument is active.
Wrapper for a generic expression template.
static double x_
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors.
@ InitDerivArray
Initialize the derivative array.
A tag for specializing Expr for SFad expressions.
Base template specification for testing equivalence.
static SACADO_INLINE_FUNCTION void zero(T *dest, int sz)
Zero out array dest of length sz.