Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_ELRFad_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_ELRFAD_GENERALFAD_HPP
53#define SACADO_ELRFAD_GENERALFAD_HPP
54
56#include "Sacado_dummy_arg.hpp"
59#include<ostream>
60
61namespace Sacado {
62
64 namespace ELRFad {
65
67
72 template <typename T, typename Storage>
73 class GeneralFad : public Storage {
74
75 public:
76
79
82
87
90 GeneralFad() : Storage(T(0.)) {}
91
93
96 template <typename S>
99 Storage(x) {}
100
102
106 GeneralFad(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) :
107 Storage(sz, x, zero_out) {}
108
110
116 GeneralFad(const int sz, const int i, const T & x) :
117 Storage(sz, x, InitDerivArray) {
118 this->fastAccessDx(i)=1.;
119 }
120
123 GeneralFad(const Storage& s) : Storage(s) {}
124
128 Storage(x) {}
129
131 template <typename S>
134 Storage(x.size(), T(0.), NoInitDerivArray) {
135 const int sz = x.size();
136 if (sz) {
137
138 if (Expr<S>::is_linear) {
139 if (x.hasFastAccess())
140 for(int i=0; i<sz; ++i)
141 this->fastAccessDx(i) = x.fastAccessDx(i);
142 else
143 for(int i=0; i<sz; ++i)
144 this->fastAccessDx(i) = x.dx(i);
145 }
146 else {
147
148 if (x.hasFastAccess()) {
149 // Compute partials
150 FastLocalAccumOp< Expr<S> > op(x);
151
152 // Compute each tangent direction
153 for(op.i=0; op.i<sz; ++op.i) {
154 op.t = T(0.);
155
156 // Automatically unrolled loop that computes
157 // for (int j=0; j<N; j++)
158 // op.t += op.partials[j] * x.getTangent<j>(i);
160
161 this->fastAccessDx(op.i) = op.t;
162 }
163 }
164 else {
165 // Compute partials
167
168 // Compute each tangent direction
169 for(op.i=0; op.i<sz; ++op.i) {
170 op.t = T(0.);
171
172 // Automatically unrolled loop that computes
173 // for (int j=0; j<N; j++)
174 // op.t += op.partials[j] * x.getTangent<j>(i);
176
177 this->fastAccessDx(op.i) = op.t;
178 }
179 }
180
181 }
182
183 }
184
185 // Compute value
186 this->val() = x.val();
187 }
188
192
194
201 void diff(const int ith, const int n) {
202 if (this->size() != n)
203 this->resize(n);
204
205 this->zero();
206 this->fastAccessDx(ith) = T(1.);
207
208 }
209
212 void setUpdateValue(bool update_val) { }
213
216 bool updateValue() const { return true; }
217
220 void cache() const {}
221
223 template <typename S>
225 SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
226 typedef IsEqual<value_type> IE;
227 if (x.size() != this->size()) return false;
228 bool eq = IE::eval(x.val(), this->val());
229 for (int i=0; i<this->size(); i++)
230 eq = eq && IE::eval(x.dx(i), this->dx(i));
231 return eq;
232 }
233
235
240
246 int availableSize() const { return this->length(); }
247
250 bool hasFastAccess() const { return this->size()!=0;}
251
254 bool isPassive() const { return this->size()==0;}
255
258 void setIsConstant(bool is_const) {
259 if (is_const && this->size()!=0)
260 this->resize(0);
261 }
262
264
269
271 template <typename S>
273 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator=(const S& v) {
274 this->val() = v;
275 if (this->size()) this->resize(0);
276 return *this;
277 }
278
282 operator=(const GeneralFad& x) {
283 // Copy value & dx_
284 Storage::operator=(x);
285 return *this;
286 }
287
289 template <typename S>
291 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator=(const Expr<S>& x) {
292 const int xsz = x.size();
293 if (xsz != this->size())
294 this->resizeAndZero(xsz);
295
296 const int sz = this->size();
297
298 // For ViewStorage, the resize above may not in fact resize the
299 // derivative array, so it is possible that sz != xsz at this point.
300 // The only valid use case here is sz > xsz == 0, so we use sz in the
301 // assignment below
302
303 if (sz) {
304
305 if (Expr<S>::is_linear) {
306 if (x.hasFastAccess())
307 for(int i=0; i<sz; ++i)
308 this->fastAccessDx(i) = x.fastAccessDx(i);
309 else
310 for(int i=0; i<sz; ++i)
311 this->fastAccessDx(i) = x.dx(i);
312 }
313 else {
314
315 if (x.hasFastAccess()) {
316 // Compute partials
317 FastLocalAccumOp< Expr<S> > op(x);
318
319 // Compute each tangent direction
320 for(op.i=0; op.i<sz; ++op.i) {
321 op.t = T(0.);
322
323 // Automatically unrolled loop that computes
324 // for (int j=0; j<N; j++)
325 // op.t += op.partials[j] * x.getTangent<j>(i);
327
328 this->fastAccessDx(op.i) = op.t;
329 }
330 }
331 else {
332 // Compute partials
333 SlowLocalAccumOp< Expr<S> > op(x);
334
335 // Compute each tangent direction
336 for(op.i=0; op.i<sz; ++op.i) {
337 op.t = T(0.);
338
339 // Automatically unrolled loop that computes
340 // for (int j=0; j<N; j++)
341 // op.t += op.partials[j] * x.getTangent<j>(i);
343
344 this->fastAccessDx(op.i) = op.t;
345 }
346 }
347 }
348 }
349
350 // Compute value
351 this->val() = x.val();
352
353 return *this;
354 }
355
357
362
364 template <typename S>
366 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
367 this->val() += v;
368 return *this;
369 }
370
372 template <typename S>
374 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator -= (const S& v) {
375 this->val() -= v;
376 return *this;
377 }
378
380 template <typename S>
382 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator *= (const S& v) {
383 const int sz = this->size();
384 this->val() *= v;
385 for (int i=0; i<sz; ++i)
386 this->fastAccessDx(i) *= v;
387 return *this;
388 }
389
391 template <typename S>
393 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
394 const int sz = this->size();
395 this->val() /= v;
396 for (int i=0; i<sz; ++i)
397 this->fastAccessDx(i) /= v;
398 return *this;
399 }
400
403 GeneralFad& operator += (const GeneralFad& x) {
404 const int xsz = x.size(), sz = this->size();
405
406#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
407 if ((xsz != sz) && (xsz != 0) && (sz != 0))
408 throw "Fad Error: Attempt to assign with incompatible sizes";
409#endif
410
411 if (xsz) {
412 if (sz) {
413 for (int i=0; i<sz; ++i)
414 this->fastAccessDx(i) += x.fastAccessDx(i);
415 }
416 else {
417 this->resizeAndZero(xsz);
418 for (int i=0; i<xsz; ++i)
419 this->fastAccessDx(i) = x.fastAccessDx(i);
420 }
421 }
422
423 this->val() += x.val();
424
425 return *this;
426 }
427
430 GeneralFad& operator -= (const GeneralFad& x) {
431 const int xsz = x.size(), sz = this->size();
432
433#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
434 if ((xsz != sz) && (xsz != 0) && (sz != 0))
435 throw "Fad Error: Attempt to assign with incompatible sizes";
436#endif
437
438 if (xsz) {
439 if (sz) {
440 for(int i=0; i<sz; ++i)
441 this->fastAccessDx(i) -= x.fastAccessDx(i);
442 }
443 else {
444 this->resizeAndZero(xsz);
445 for(int i=0; i<xsz; ++i)
446 this->fastAccessDx(i) = -x.fastAccessDx(i);
447 }
448 }
449
450 this->val() -= x.val();
451
452
453 return *this;
454 }
455
458 GeneralFad& operator *= (const GeneralFad& x) {
459 const int xsz = x.size(), sz = this->size();
460 T xval = x.val();
461 T v = this->val();
462
463#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
464 if ((xsz != sz) && (xsz != 0) && (sz != 0))
465 throw "Fad Error: Attempt to assign with incompatible sizes";
466#endif
467
468 if (xsz) {
469 if (sz) {
470 for(int i=0; i<sz; ++i)
471 this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
472 }
473 else {
474 this->resizeAndZero(xsz);
475 for(int i=0; i<xsz; ++i)
476 this->fastAccessDx(i) = v*x.fastAccessDx(i);
477 }
478 }
479 else {
480 if (sz) {
481 for (int i=0; i<sz; ++i)
482 this->fastAccessDx(i) *= xval;
483 }
484 }
485
486 this->val() *= xval;
487
488 return *this;
489 }
490
493 GeneralFad& operator /= (const GeneralFad& x) {
494 const int xsz = x.size(), sz = this->size();
495 T xval = x.val();
496 T v = this->val();
497
498#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
499 if ((xsz != sz) && (xsz != 0) && (sz != 0))
500 throw "Fad Error: Attempt to assign with incompatible sizes";
501#endif
502
503 if (xsz) {
504 if (sz) {
505 for(int i=0; i<sz; ++i)
506 this->fastAccessDx(i) =
507 ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
508 }
509 else {
510 this->resizeAndZero(xsz);
511 for(int i=0; i<xsz; ++i)
512 this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
513 }
514 }
515 else {
516 if (sz) {
517 for (int i=0; i<sz; ++i)
518 this->fastAccessDx(i) /= xval;
519 }
520 }
521
522 this->val() /= xval;
523
524 return *this;
525 }
526
528 template <typename S>
530 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator += (const Expr<S>& x) {
531 const int xsz = x.size(), sz = this->size();
532
533#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
534 if ((xsz != sz) && (xsz != 0) && (sz != 0))
535 throw "Fad Error: Attempt to assign with incompatible sizes";
536#endif
537
538 if (Expr<S>::is_linear) {
539 if (xsz) {
540 if (sz) {
541 if (x.hasFastAccess())
542 for (int i=0; i<sz; ++i)
543 this->fastAccessDx(i) += x.fastAccessDx(i);
544 else
545 for (int i=0; i<sz; ++i)
546 this->fastAccessDx(i) += x.dx(i);
547 }
548 else {
549 this->resizeAndZero(xsz);
550 if (x.hasFastAccess())
551 for (int i=0; i<xsz; ++i)
552 this->fastAccessDx(i) = x.fastAccessDx(i);
553 else
554 for (int i=0; i<xsz; ++i)
555 this->fastAccessDx(i) = x.dx(i);
556 }
557 }
558 }
559 else {
560
561 if (xsz) {
562
563 if (sz != xsz)
564 this->resizeAndZero(xsz);
565
566 if (x.hasFastAccess()) {
567 // Compute partials
568 FastLocalAccumOp< Expr<S> > op(x);
569
570 // Compute each tangent direction
571 for(op.i=0; op.i<xsz; ++op.i) {
572 op.t = T(0.);
573
574 // Automatically unrolled loop that computes
575 // for (int j=0; j<N; j++)
576 // op.t += op.partials[j] * x.getTangent<j>(i);
578
579 this->fastAccessDx(op.i) += op.t;
580 }
581 }
582 else {
583 // Compute partials
584 SlowLocalAccumOp< Expr<S> > op(x);
585
586 // Compute each tangent direction
587 for(op.i=0; op.i<xsz; ++op.i) {
588 op.t = T(0.);
589
590 // Automatically unrolled loop that computes
591 // for (int j=0; j<N; j++)
592 // op.t += op.partials[j] * x.getTangent<j>(i);
594
595 this->fastAccessDx(op.i) += op.t;
596 }
597 }
598
599 }
600
601 }
602
603 // Compute value
604 this->val() += x.val();
605
606 return *this;
607 }
608
610 template <typename S>
612 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
613 const int xsz = x.size(), sz = this->size();
614
615#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
616 if ((xsz != sz) && (xsz != 0) && (sz != 0))
617 throw "Fad Error: Attempt to assign with incompatible sizes";
618#endif
619
620 if (Expr<S>::is_linear) {
621 if (xsz) {
622 if (sz) {
623 if (x.hasFastAccess())
624 for(int i=0; i<sz; ++i)
625 this->fastAccessDx(i) -= x.fastAccessDx(i);
626 else
627 for (int i=0; i<sz; ++i)
628 this->fastAccessDx(i) -= x.dx(i);
629 }
630 else {
631 this->resizeAndZero(xsz);
632 if (x.hasFastAccess())
633 for(int i=0; i<xsz; ++i)
634 this->fastAccessDx(i) = -x.fastAccessDx(i);
635 else
636 for (int i=0; i<xsz; ++i)
637 this->fastAccessDx(i) = -x.dx(i);
638 }
639 }
640 }
641 else {
642
643 if (xsz) {
644
645 if (sz != xsz)
646 this->resizeAndZero(xsz);
647
648 if (x.hasFastAccess()) {
649 // Compute partials
650 FastLocalAccumOp< Expr<S> > op(x);
651
652 // Compute each tangent direction
653 for(op.i=0; op.i<xsz; ++op.i) {
654 op.t = T(0.);
655
656 // Automatically unrolled loop that computes
657 // for (int j=0; j<N; j++)
658 // op.t += op.partials[j] * x.getTangent<j>(i);
660
661 this->fastAccessDx(op.i) -= op.t;
662 }
663 }
664 else {
665 // Compute partials
666 SlowLocalAccumOp< Expr<S> > op(x);
667
668 // Compute each tangent direction
669 for(op.i=0; op.i<xsz; ++op.i) {
670 op.t = T(0.);
671
672 // Automatically unrolled loop that computes
673 // for (int j=0; j<N; j++)
674 // op.t += op.partials[j] * x.getTangent<j>(i);
676
677 this->fastAccessDx(op.i) -= op.t;
678 }
679 }
680 }
681
682 }
683
684 this->val() -= x.val();
685
686 return *this;
687 }
688
690 template <typename S>
692 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
693 const int xsz = x.size(), sz = this->size();
694 T xval = x.val();
695 T v = this->val();
696
697#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
698 if ((xsz != sz) && (xsz != 0) && (sz != 0))
699 throw "Fad Error: Attempt to assign with incompatible sizes";
700#endif
701
702 if (Expr<S>::is_linear) {
703 if (xsz) {
704 if (sz) {
705 if (x.hasFastAccess())
706 for(int i=0; i<sz; ++i)
707 this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
708 else
709 for (int i=0; i<sz; ++i)
710 this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
711 }
712 else {
713 this->resizeAndZero(xsz);
714 if (x.hasFastAccess())
715 for(int i=0; i<xsz; ++i)
716 this->fastAccessDx(i) = v*x.fastAccessDx(i);
717 else
718 for (int i=0; i<xsz; ++i)
719 this->fastAccessDx(i) = v*x.dx(i);
720 }
721 }
722 else {
723 if (sz) {
724 for (int i=0; i<sz; ++i)
725 this->fastAccessDx(i) *= xval;
726 }
727 }
728 }
729 else {
730
731 if (xsz) {
732
733 if (sz) {
734
735 if (x.hasFastAccess()) {
736 // Compute partials
737 FastLocalAccumOp< Expr<S> > op(x);
738
739 // Compute each tangent direction
740 for(op.i=0; op.i<xsz; ++op.i) {
741 op.t = T(0.);
742
743 // Automatically unrolled loop that computes
744 // for (int j=0; j<N; j++)
745 // op.t += op.partials[j] * x.getTangent<j>(i);
747
748 this->fastAccessDx(op.i) =
749 v * op.t + this->fastAccessDx(op.i) * xval;
750 }
751 }
752 else {
753 // Compute partials
754 SlowLocalAccumOp< Expr<S> > op(x);
755
756 // Compute each tangent direction
757 for(op.i=0; op.i<xsz; ++op.i) {
758 op.t = T(0.);
759
760 // Automatically unrolled loop that computes
761 // for (int j=0; j<N; j++)
762 // op.t += op.partials[j] * x.getTangent<j>(i);
764
765 this->fastAccessDx(op.i) =
766 v * op.t + this->fastAccessDx(op.i) * xval;
767 }
768 }
769
770 }
771
772 else {
773
774 this->resizeAndZero(xsz);
775
776 if (x.hasFastAccess()) {
777 // Compute partials
778 FastLocalAccumOp< Expr<S> > op(x);
779
780 // Compute each tangent direction
781 for(op.i=0; op.i<xsz; ++op.i) {
782 op.t = T(0.);
783
784 // Automatically unrolled loop that computes
785 // for (int j=0; j<N; j++)
786 // op.t += op.partials[j] * x.getTangent<j>(i);
788
789 this->fastAccessDx(op.i) = v * op.t;
790 }
791 }
792 else {
793 // Compute partials
794 SlowLocalAccumOp< Expr<S> > op(x);
795
796 // Compute each tangent direction
797 for(op.i=0; op.i<xsz; ++op.i) {
798 op.t = T(0.);
799
800 // Automatically unrolled loop that computes
801 // for (int j=0; j<N; j++)
802 // op.t += op.partials[j] * x.getTangent<j>(i);
804
805 this->fastAccessDx(op.i) = v * op.t;
806 }
807 }
808
809 }
810
811 }
812
813 else {
814
815 if (sz) {
816 for (int i=0; i<sz; ++i)
817 this->fastAccessDx(i) *= xval;
818 }
819
820 }
821
822 }
823
824 this->val() *= xval;
825
826 return *this;
827 }
828
830 template <typename S>
832 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
833 const int xsz = x.size(), sz = this->size();
834 T xval = x.val();
835 T v = this->val();
836
837#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
838 if ((xsz != sz) && (xsz != 0) && (sz != 0))
839 throw "Fad Error: Attempt to assign with incompatible sizes";
840#endif
841
842 if (Expr<S>::is_linear) {
843 if (xsz) {
844 if (sz) {
845 if (x.hasFastAccess())
846 for(int i=0; i<sz; ++i)
847 this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
848 else
849 for (int i=0; i<sz; ++i)
850 this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
851 }
852 else {
853 this->resizeAndZero(xsz);
854 if (x.hasFastAccess())
855 for(int i=0; i<xsz; ++i)
856 this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
857 else
858 for (int i=0; i<xsz; ++i)
859 this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
860 }
861 }
862 else {
863 if (sz) {
864 for (int i=0; i<sz; ++i)
865 this->fastAccessDx(i) /= xval;
866 }
867 }
868 }
869 else {
870
871 if (xsz) {
872
873 T xval2 = xval*xval;
874
875 if (sz) {
876
877 if (x.hasFastAccess()) {
878 // Compute partials
879 FastLocalAccumOp< Expr<S> > op(x);
880
881 // Compute each tangent direction
882 for(op.i=0; op.i<xsz; ++op.i) {
883 op.t = T(0.);
884
885 // Automatically unrolled loop that computes
886 // for (int j=0; j<N; j++)
887 // op.t += op.partials[j] * x.getTangent<j>(i);
889
890 this->fastAccessDx(op.i) =
891 (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
892 }
893 }
894 else {
895 // Compute partials
896 SlowLocalAccumOp< Expr<S> > op(x);
897
898 // Compute each tangent direction
899 for(op.i=0; op.i<xsz; ++op.i) {
900 op.t = T(0.);
901
902 // Automatically unrolled loop that computes
903 // for (int j=0; j<N; j++)
904 // op.t += op.partials[j] * x.getTangent<j>(i);
906
907 this->fastAccessDx(op.i) =
908 (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
909 }
910 }
911
912 }
913
914 else {
915
916 this->resizeAndZero(xsz);
917
918 if (x.hasFastAccess()) {
919 // Compute partials
920 FastLocalAccumOp< Expr<S> > op(x);
921
922 // Compute each tangent direction
923 for(op.i=0; op.i<xsz; ++op.i) {
924 op.t = T(0.);
925
926 // Automatically unrolled loop that computes
927 // for (int j=0; j<N; j++)
928 // op.t += op.partials[j] * x.getTangent<j>(i);
930
931 this->fastAccessDx(op.i) = (-v * op.t) / xval2;
932 }
933 }
934 else {
935 // Compute partials
936 SlowLocalAccumOp< Expr<S> > op(x);
937
938 // Compute each tangent direction
939 for(op.i=0; op.i<xsz; ++op.i) {
940 op.t = T(0.);
941
942 // Automatically unrolled loop that computes
943 // for (int j=0; j<N; j++)
944 // op.t += op.partials[j] * x.getTangent<j>(i);
946
947 this->fastAccessDx(op.i) = (-v * op.t) / xval2;
948 }
949 }
950
951 }
952
953 }
954
955 else {
956
957 if (sz) {
958 for (int i=0; i<sz; ++i)
959 this->fastAccessDx(i) /= xval;
960 }
961
962 }
963
964 }
965
966 this->val() /= xval;
967
968 return *this;
969 }
970
972
973 protected:
974
975 // Functor for mpl::for_each to compute the local accumulation
976 // of a tangent derivative
977 //
978 // We use getTangent<>() to get dx components from expression
979 // arguments instead of getting the argument directly or extracting
980 // the dx array due to striding in ViewFad (or could use striding
981 // directly here if we need to get dx array).
982 template <typename ExprT>
983 struct FastLocalAccumOp {
984 typedef typename ExprT::value_type value_type;
985 static const int N = ExprT::num_args;
986 const ExprT& x;
987 mutable value_type t;
988 value_type partials[N];
989 int i;
991 FastLocalAccumOp(const ExprT& x_) : x(x_) {
992 x.computePartials(value_type(1.), partials);
993 }
994 template <typename ArgT>
996 void operator () (ArgT arg) const {
997 const int Arg = ArgT::value;
998 t += partials[Arg] * x.template getTangent<Arg>(i);
999 }
1000 };
1001
1002 template <typename ExprT>
1003 struct SlowLocalAccumOp : FastLocalAccumOp<ExprT> {
1005 SlowLocalAccumOp(const ExprT& x_) :
1006 FastLocalAccumOp<ExprT>(x_) {}
1007 template <typename ArgT>
1009 void operator () (ArgT arg) const {
1010 const int Arg = ArgT::value;
1011 if (this->x.template isActive<Arg>())
1012 this->t += this->partials[Arg] * this->x.template getTangent<Arg>(this->i);
1013 }
1014 };
1015
1016 }; // class GeneralFad
1017
1018
1019 template <typename T, typename Storage>
1020 std::ostream& operator << (std::ostream& os,
1021 const GeneralFad<T,Storage>& x) {
1022 os << x.val() << " [";
1023
1024 for (int i=0; i< x.size(); i++) {
1025 os << " " << x.dx(i);
1026 }
1027
1028 os << " ]";
1029 return os;
1030 }
1031
1032 } // namespace ELRFad
1033
1034} // namespace Sacado
1035
1036#endif // SACADO_ELRFAD_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
const int N
Wrapper for a generic expression template.
Forward-mode AD class templated on the storage for the derivative array.
SACADO_INLINE_FUNCTION GeneralFad(const GeneralFad &x)
Copy constructor.
SACADO_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
SACADO_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
SACADO_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
SACADO_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
SACADO_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.
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 GeneralFad(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
SACADO_INLINE_FUNCTION SACADO_ENABLE_VALUE_FUNC(GeneralFad &) operator
Assignment operator with constant right-hand-side.
ScalarType< value_type >::type scalar_type
Typename of scalar's (which may be different from T)
SACADO_INLINE_FUNCTION void cache() const
Cache values.
SACADO_INLINE_FUNCTION ~GeneralFad()
Destructor.
SACADO_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
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 GeneralFad(const Storage &s)
Constructor with supplied storage s.
SACADO_INLINE_FUNCTION GeneralFad()
Default 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 bool isPassive() const
Returns true if derivative array is empty.
SACADO_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
static double x_
std::ostream & operator<<(std::ostream &os, const GeneralFad< T, Storage > &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.
SACADO_INLINE_FUNCTION SlowLocalAccumOp(const ExprT &x_)
SACADO_INLINE_FUNCTION void operator()(ArgT arg) const
Base template specification for testing equivalence.