Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_ELRCacheFad_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_ELRCACHEFAD_GENERALFAD_HPP
53#define SACADO_ELRCACHEFAD_GENERALFAD_HPP
54
58
59namespace Sacado {
60
62 namespace ELRCacheFad {
63
65
70 template <typename T, typename Storage>
71 class GeneralFad : public Storage {
72
73 public:
74
77
80
85
88 GeneralFad() : Storage(T(0.)) {}
89
91
94 template <typename S>
97 Storage(x) {}
98
100
104 GeneralFad(const int sz, const T & x, const DerivInit zero_out = InitDerivArray) :
105 Storage(sz, x, zero_out) {}
106
108
114 GeneralFad(const int sz, const int i, const T & x) :
115 Storage(sz, x, InitDerivArray) {
116 this->fastAccessDx(i)=1.;
117 }
118
121 GeneralFad(const Storage& s) : Storage(s) {}
122
126 Storage(x) {}
127
129 template <typename S>
132 Storage(x.size(), T(0.), NoInitDerivArray) {
133 x.cache();
134
135 const int sz = x.size();
136
137 // Compute value
138 this->val() = x.val();
139
140 if (sz) {
141
142 if (Expr<S>::is_linear) {
143 if (x.hasFastAccess())
144 for(int i=0; i<sz; ++i)
145 this->fastAccessDx(i) = x.fastAccessDx(i);
146 else
147 for(int i=0; i<sz; ++i)
148 this->fastAccessDx(i) = x.dx(i);
149 }
150 else {
151
152 if (x.hasFastAccess()) {
153 // Compute partials
154 FastLocalAccumOp< Expr<S> > op(x);
155
156 // Compute each tangent direction
157 for(op.i=0; op.i<sz; ++op.i) {
158 op.t = T(0.);
159
160 // Automatically unrolled loop that computes
161 // for (int j=0; j<N; j++)
162 // op.t += op.partials[j] * x.getTangent<j>(i);
164
165 this->fastAccessDx(op.i) = op.t;
166 }
167 }
168 else {
169 // Compute partials
171
172 // Compute each tangent direction
173 for(op.i=0; op.i<sz; ++op.i) {
174 op.t = T(0.);
175
176 // Automatically unrolled loop that computes
177 // for (int j=0; j<N; j++)
178 // op.t += op.partials[j] * x.getTangent<j>(i);
180
181 this->fastAccessDx(op.i) = op.t;
182 }
183 }
184
185 }
186 }
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 x.cache();
293
294 const int xsz = x.size();
295
296 if (xsz != this->size())
297 this->resizeAndZero(xsz);
298
299 const int sz = this->size();
300
301 // For ViewStorage, the resize above may not in fact resize the
302 // derivative array, so it is possible that sz != xsz at this point.
303 // The only valid use case here is sz > xsz == 0, so we use sz in the
304 // assignment below
305
306 if (sz) {
307
308 if (Expr<S>::is_linear) {
309 if (x.hasFastAccess())
310 for(int i=0; i<sz; ++i)
311 this->fastAccessDx(i) = x.fastAccessDx(i);
312 else
313 for(int i=0; i<sz; ++i)
314 this->fastAccessDx(i) = x.dx(i);
315 }
316 else {
317
318 if (x.hasFastAccess()) {
319 // Compute partials
320 FastLocalAccumOp< Expr<S> > op(x);
321
322 // Compute each tangent direction
323 for(op.i=0; op.i<sz; ++op.i) {
324 op.t = T(0.);
325
326 // Automatically unrolled loop that computes
327 // for (int j=0; j<N; j++)
328 // op.t += op.partials[j] * x.getTangent<j>(i);
330
331 this->fastAccessDx(op.i) = op.t;
332 }
333 }
334 else {
335 // Compute partials
336 SlowLocalAccumOp< Expr<S> > op(x);
337
338 // Compute each tangent direction
339 for(op.i=0; op.i<sz; ++op.i) {
340 op.t = T(0.);
341
342 // Automatically unrolled loop that computes
343 // for (int j=0; j<N; j++)
344 // op.t += op.partials[j] * x.getTangent<j>(i);
346
347 this->fastAccessDx(op.i) = op.t;
348 }
349 }
350
351 }
352
353 }
354
355 this->val() = x.val();
356
357 return *this;
358 }
359
361
366
368 template <typename S>
370 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
371 this->val() += v;
372 return *this;
373 }
374
376 template <typename S>
378 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator -= (const S& v) {
379 this->val() -= v;
380 return *this;
381 }
382
384 template <typename S>
386 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator *= (const S& v) {
387 const int sz = this->size();
388 this->val() *= v;
389 for (int i=0; i<sz; ++i)
390 this->fastAccessDx(i) *= v;
391 return *this;
392 }
393
395 template <typename S>
397 SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
398 const int sz = this->size();
399 this->val() /= v;
400 for (int i=0; i<sz; ++i)
401 this->fastAccessDx(i) /= v;
402 return *this;
403 }
404
407 GeneralFad& operator += (const GeneralFad& x) {
408 const int xsz = x.size(), sz = this->size();
409
410#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
411 if ((xsz != sz) && (xsz != 0) && (sz != 0))
412 throw "Fad Error: Attempt to assign with incompatible sizes";
413#endif
414
415 if (xsz) {
416 if (sz) {
417 for (int i=0; i<sz; ++i)
418 this->fastAccessDx(i) += x.fastAccessDx(i);
419 }
420 else {
421 this->resizeAndZero(xsz);
422 for (int i=0; i<xsz; ++i)
423 this->fastAccessDx(i) = x.fastAccessDx(i);
424 }
425 }
426
427 this->val() += x.val();
428
429 return *this;
430 }
431
434 GeneralFad& operator -= (const GeneralFad& x) {
435 const int xsz = x.size(), sz = this->size();
436
437#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
438 if ((xsz != sz) && (xsz != 0) && (sz != 0))
439 throw "Fad Error: Attempt to assign with incompatible sizes";
440#endif
441
442 if (xsz) {
443 if (sz) {
444 for(int i=0; i<sz; ++i)
445 this->fastAccessDx(i) -= x.fastAccessDx(i);
446 }
447 else {
448 this->resizeAndZero(xsz);
449 for(int i=0; i<xsz; ++i)
450 this->fastAccessDx(i) = -x.fastAccessDx(i);
451 }
452 }
453
454 this->val() -= x.val();
455
456
457 return *this;
458 }
459
462 GeneralFad& operator *= (const GeneralFad& x) {
463 const int xsz = x.size(), sz = this->size();
464 T xval = x.val();
465 T v = this->val();
466
467#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
468 if ((xsz != sz) && (xsz != 0) && (sz != 0))
469 throw "Fad Error: Attempt to assign with incompatible sizes";
470#endif
471
472 if (xsz) {
473 if (sz) {
474 for(int i=0; i<sz; ++i)
475 this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
476 }
477 else {
478 this->resizeAndZero(xsz);
479 for(int i=0; i<xsz; ++i)
480 this->fastAccessDx(i) = v*x.fastAccessDx(i);
481 }
482 }
483 else {
484 if (sz) {
485 for (int i=0; i<sz; ++i)
486 this->fastAccessDx(i) *= xval;
487 }
488 }
489
490 this->val() *= xval;
491
492 return *this;
493 }
494
497 GeneralFad& operator /= (const GeneralFad& x) {
498 const int xsz = x.size(), sz = this->size();
499 T xval = x.val();
500 T v = this->val();
501
502#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
503 if ((xsz != sz) && (xsz != 0) && (sz != 0))
504 throw "Fad Error: Attempt to assign with incompatible sizes";
505#endif
506
507 if (xsz) {
508 if (sz) {
509 for(int i=0; i<sz; ++i)
510 this->fastAccessDx(i) =
511 ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
512 }
513 else {
514 this->resizeAndZero(xsz);
515 for(int i=0; i<xsz; ++i)
516 this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
517 }
518 }
519 else {
520 if (sz) {
521 for (int i=0; i<sz; ++i)
522 this->fastAccessDx(i) /= xval;
523 }
524 }
525
526 this->val() /= xval;
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
539#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
540 if ((xsz != sz) && (xsz != 0) && (sz != 0))
541 throw "Fad Error: Attempt to assign with incompatible sizes";
542#endif
543
544 if (Expr<S>::is_linear) {
545 if (xsz) {
546 if (sz) {
547 if (x.hasFastAccess())
548 for (int i=0; i<sz; ++i)
549 this->fastAccessDx(i) += x.fastAccessDx(i);
550 else
551 for (int i=0; i<sz; ++i)
552 this->fastAccessDx(i) += x.dx(i);
553 }
554 else {
555 this->resizeAndZero(xsz);
556 if (x.hasFastAccess())
557 for (int i=0; i<xsz; ++i)
558 this->fastAccessDx(i) = x.fastAccessDx(i);
559 else
560 for (int i=0; i<xsz; ++i)
561 this->fastAccessDx(i) = x.dx(i);
562 }
563 }
564 }
565 else {
566
567 if (xsz) {
568
569 if (sz != xsz)
570 this->resizeAndZero(xsz);
571
572 if (x.hasFastAccess()) {
573 // Compute partials
574 FastLocalAccumOp< Expr<S> > op(x);
575
576 // Compute each tangent direction
577 for(op.i=0; op.i<xsz; ++op.i) {
578 op.t = T(0.);
579
580 // Automatically unrolled loop that computes
581 // for (int j=0; j<N; j++)
582 // op.t += op.partials[j] * x.getTangent<j>(i);
584
585 this->fastAccessDx(op.i) += op.t;
586 }
587 }
588 else {
589 // Compute partials
590 SlowLocalAccumOp< Expr<S> > op(x);
591
592 // Compute each tangent direction
593 for(op.i=0; op.i<xsz; ++op.i) {
594 op.t = T(0.);
595
596 // Automatically unrolled loop that computes
597 // for (int j=0; j<N; j++)
598 // op.t += op.partials[j] * x.getTangent<j>(i);
600
601 this->fastAccessDx(op.i) += op.t;
602 }
603 }
604
605 }
606
607 }
608
609 this->val() += x.val();
610
611 return *this;
612 }
613
615 template <typename S>
617 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
618 x.cache();
619
620 const int xsz = x.size(), sz = this->size();
621
622#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
623 if ((xsz != sz) && (xsz != 0) && (sz != 0))
624 throw "Fad Error: Attempt to assign with incompatible sizes";
625#endif
626
627 if (Expr<S>::is_linear) {
628 if (xsz) {
629 if (sz) {
630 if (x.hasFastAccess())
631 for(int i=0; i<sz; ++i)
632 this->fastAccessDx(i) -= x.fastAccessDx(i);
633 else
634 for (int i=0; i<sz; ++i)
635 this->fastAccessDx(i) -= x.dx(i);
636 }
637 else {
638 this->resizeAndZero(xsz);
639 if (x.hasFastAccess())
640 for(int i=0; i<xsz; ++i)
641 this->fastAccessDx(i) = -x.fastAccessDx(i);
642 else
643 for (int i=0; i<xsz; ++i)
644 this->fastAccessDx(i) = -x.dx(i);
645 }
646 }
647 }
648 else {
649
650 if (xsz) {
651
652 if (sz != xsz)
653 this->resizeAndZero(xsz);
654
655 if (x.hasFastAccess()) {
656 // Compute partials
657 FastLocalAccumOp< Expr<S> > op(x);
658
659 // Compute each tangent direction
660 for(op.i=0; op.i<xsz; ++op.i) {
661 op.t = T(0.);
662
663 // Automatically unrolled loop that computes
664 // for (int j=0; j<N; j++)
665 // op.t += op.partials[j] * x.getTangent<j>(i);
667
668 this->fastAccessDx(op.i) -= op.t;
669 }
670 }
671 else {
672 // Compute partials
673 SlowLocalAccumOp< Expr<S> > op(x);
674
675 // Compute each tangent direction
676 for(op.i=0; op.i<xsz; ++op.i) {
677 op.t = T(0.);
678
679 // Automatically unrolled loop that computes
680 // for (int j=0; j<N; j++)
681 // op.t += op.partials[j] * x.getTangent<j>(i);
683
684 this->fastAccessDx(op.i) -= op.t;
685 }
686 }
687
688 }
689
690 }
691
692 this->val() -= x.val();
693
694 return *this;
695 }
696
698 template <typename S>
700 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
701 x.cache();
702
703 const int xsz = x.size(), sz = this->size();
704 T xval = x.val();
705 T v = this->val();
706
707#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
708 if ((xsz != sz) && (xsz != 0) && (sz != 0))
709 throw "Fad Error: Attempt to assign with incompatible sizes";
710#endif
711
712 if (Expr<S>::is_linear) {
713 if (xsz) {
714 if (sz) {
715 if (x.hasFastAccess())
716 for(int i=0; i<sz; ++i)
717 this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
718 else
719 for (int i=0; i<sz; ++i)
720 this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
721 }
722 else {
723 this->resizeAndZero(xsz);
724 if (x.hasFastAccess())
725 for(int i=0; i<xsz; ++i)
726 this->fastAccessDx(i) = v*x.fastAccessDx(i);
727 else
728 for (int i=0; i<xsz; ++i)
729 this->fastAccessDx(i) = v*x.dx(i);
730 }
731 }
732 else {
733 if (sz) {
734 for (int i=0; i<sz; ++i)
735 this->fastAccessDx(i) *= xval;
736 }
737 }
738 }
739 else {
740
741 if (xsz) {
742
743 if (sz) {
744
745 if (x.hasFastAccess()) {
746 // Compute partials
747 FastLocalAccumOp< Expr<S> > op(x);
748
749 // Compute each tangent direction
750 for(op.i=0; op.i<xsz; ++op.i) {
751 op.t = T(0.);
752
753 // Automatically unrolled loop that computes
754 // for (int j=0; j<N; j++)
755 // op.t += op.partials[j] * x.getTangent<j>(i);
757
758 this->fastAccessDx(op.i) =
759 v * op.t + this->fastAccessDx(op.i) * xval;
760 }
761 }
762 else {
763 // Compute partials
764 SlowLocalAccumOp< Expr<S> > op(x);
765
766 // Compute each tangent direction
767 for(op.i=0; op.i<xsz; ++op.i) {
768 op.t = T(0.);
769
770 // Automatically unrolled loop that computes
771 // for (int j=0; j<N; j++)
772 // op.t += op.partials[j] * x.getTangent<j>(i);
774
775 this->fastAccessDx(op.i) =
776 v * op.t + this->fastAccessDx(op.i) * xval;
777 }
778 }
779
780 }
781
782 else {
783
784 this->resizeAndZero(xsz);
785
786 if (x.hasFastAccess()) {
787 // Compute partials
788 FastLocalAccumOp< Expr<S> > op(x);
789
790 // Compute each tangent direction
791 for(op.i=0; op.i<xsz; ++op.i) {
792 op.t = T(0.);
793
794 // Automatically unrolled loop that computes
795 // for (int j=0; j<N; j++)
796 // op.t += op.partials[j] * x.getTangent<j>(i);
798
799 this->fastAccessDx(op.i) = v * op.t;
800 }
801 }
802 else {
803 // Compute partials
804 SlowLocalAccumOp< Expr<S> > op(x);
805
806 // Compute each tangent direction
807 for(op.i=0; op.i<xsz; ++op.i) {
808 op.t = T(0.);
809
810 // Automatically unrolled loop that computes
811 // for (int j=0; j<N; j++)
812 // op.t += op.partials[j] * x.getTangent<j>(i);
814
815 this->fastAccessDx(op.i) = v * op.t;
816 }
817 }
818
819 }
820
821 }
822
823 else {
824
825 if (sz) {
826 for (int i=0; i<sz; ++i)
827 this->fastAccessDx(i) *= xval;
828 }
829
830 }
831
832 }
833
834 this->val() *= xval;
835
836 return *this;
837 }
838
840 template <typename S>
842 SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
843 x.cache();
844
845 const int xsz = x.size(), sz = this->size();
846 T xval = x.val();
847 T v = this->val();
848
849#if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
850 if ((xsz != sz) && (xsz != 0) && (sz != 0))
851 throw "Fad Error: Attempt to assign with incompatible sizes";
852#endif
853
854 if (Expr<S>::is_linear) {
855 if (xsz) {
856 if (sz) {
857 if (x.hasFastAccess())
858 for(int i=0; i<sz; ++i)
859 this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
860 else
861 for (int i=0; i<sz; ++i)
862 this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
863 }
864 else {
865 this->resizeAndZero(xsz);
866 if (x.hasFastAccess())
867 for(int i=0; i<xsz; ++i)
868 this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
869 else
870 for (int i=0; i<xsz; ++i)
871 this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
872 }
873 }
874 else {
875 if (sz) {
876 for (int i=0; i<sz; ++i)
877 this->fastAccessDx(i) /= xval;
878 }
879 }
880 }
881 else {
882
883 if (xsz) {
884
885 T xval2 = xval*xval;
886
887 if (sz) {
888
889 if (x.hasFastAccess()) {
890 // Compute partials
891 FastLocalAccumOp< Expr<S> > op(x);
892
893 // Compute each tangent direction
894 for(op.i=0; op.i<xsz; ++op.i) {
895 op.t = T(0.);
896
897 // Automatically unrolled loop that computes
898 // for (int j=0; j<N; j++)
899 // op.t += op.partials[j] * x.getTangent<j>(i);
901
902 this->fastAccessDx(op.i) =
903 (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
904 }
905 }
906 else {
907 // Compute partials
908 SlowLocalAccumOp< Expr<S> > op(x);
909
910 // Compute each tangent direction
911 for(op.i=0; op.i<xsz; ++op.i) {
912 op.t = T(0.);
913
914 // Automatically unrolled loop that computes
915 // for (int j=0; j<N; j++)
916 // op.t += op.partials[j] * x.getTangent<j>(i);
918
919 this->fastAccessDx(op.i) =
920 (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
921 }
922 }
923
924 }
925
926 else {
927
928 this->resizeAndZero(xsz);
929
930 if (x.hasFastAccess()) {
931 // Compute partials
932 FastLocalAccumOp< Expr<S> > op(x);
933
934 // Compute each tangent direction
935 for(op.i=0; op.i<xsz; ++op.i) {
936 op.t = T(0.);
937
938 // Automatically unrolled loop that computes
939 // for (int j=0; j<N; j++)
940 // op.t += op.partials[j] * x.getTangent<j>(i);
942
943 this->fastAccessDx(op.i) = (-v * op.t) / xval2;
944 }
945 }
946 else {
947 // Compute partials
948 SlowLocalAccumOp< Expr<S> > op(x);
949
950 // Compute each tangent direction
951 for(op.i=0; op.i<xsz; ++op.i) {
952 op.t = T(0.);
953
954 // Automatically unrolled loop that computes
955 // for (int j=0; j<N; j++)
956 // op.t += op.partials[j] * x.getTangent<j>(i);
958
959 this->fastAccessDx(op.i) = (-v * op.t) / xval2;
960 }
961 }
962
963 }
964
965 }
966
967 else {
968
969 if (sz) {
970 for (int i=0; i<sz; ++i)
971 this->fastAccessDx(i) /= xval;
972 }
973
974 }
975
976 }
977
978 this->val() /= xval;
979
980 return *this;
981 }
982
984
985 protected:
986
987 // Functor for mpl::for_each to compute the local accumulation
988 // of a tangent derivative
989 //
990 // We use getTangent<>() to get dx components from expression
991 // arguments instead of getting the argument directly or extracting
992 // the dx array due to striding in ViewFad (or could use striding
993 // directly here if we need to get dx array).
994 template <typename ExprT>
995 struct FastLocalAccumOp {
996 typedef typename ExprT::value_type value_type;
997 static const int N = ExprT::num_args;
998 const ExprT& x;
999 mutable value_type t;
1000 value_type partials[N];
1001 int i;
1003 FastLocalAccumOp(const ExprT& x_) : x(x_) {
1004 x.computePartials(value_type(1.), partials);
1005 }
1006 template <typename ArgT>
1008 void operator () (ArgT arg) const {
1009 const int Arg = ArgT::value;
1010 t += partials[Arg] * x.template getTangent<Arg>(i);
1011 }
1012 };
1013
1014 template <typename ExprT>
1015 struct SlowLocalAccumOp : FastLocalAccumOp<ExprT> {
1017 SlowLocalAccumOp(const ExprT& x_) :
1018 FastLocalAccumOp<ExprT>(x_) {}
1019 template <typename ArgT>
1021 void operator () (ArgT arg) const {
1022 const int Arg = ArgT::value;
1023 if (this->x.template isActive<Arg>())
1024 this->t += this->partials[Arg] * this->x.template getTangent<Arg>(this->i);
1025 }
1026 };
1027
1028 }; // class GeneralFad
1029
1030
1031 template <typename T, typename Storage>
1032 std::ostream& operator << (std::ostream& os,
1033 const GeneralFad<T,Storage>& x) {
1034 os << x.val() << " [";
1035
1036 for (int i=0; i< x.size(); i++) {
1037 os << " " << x.dx(i);
1038 }
1039
1040 os << " ]";
1041 return os;
1042 }
1043
1044 } // namespace ELRCacheFad
1045
1046} // namespace Sacado
1047
1048#endif // SACADO_ELRCACHEFAD_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
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 int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
ScalarType< value_type >::type scalar_type
Typename of scalar's (which may be different from T)
RemoveConst< T >::type value_type
Typename of values.
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 bool isPassive() const
Returns true if derivative array is empty.
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 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 GeneralFad &x)
Copy constructor.
SACADO_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
SACADO_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
SACADO_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
SACADO_INLINE_FUNCTION ~GeneralFad()
Destructor.
SACADO_INLINE_FUNCTION void cache() const
Cache values.
SACADO_INLINE_FUNCTION GeneralFad(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
SACADO_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
SACADO_INLINE_FUNCTION SACADO_ENABLE_VALUE_FUNC(GeneralFad &) operator
Assignment operator with constant right-hand-side.
SACADO_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
SACADO_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
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.
Base template specification for testing equivalence.