Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
56 #include "Sacado_mpl_range_c.hpp"
57 #include "Sacado_mpl_for_each.hpp"
58 
59 namespace Sacado {
60 
62  namespace ELRCacheFad {
63 
65 
70  template <typename T, typename Storage>
71  class GeneralFad : public Storage {
72 
73  public:
74 
76  typedef typename RemoveConst<T>::type value_type;
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 
125  GeneralFad(const GeneralFad& x) :
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 
281  GeneralFad&
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
KOKKOS_INLINE_FUNCTION SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr< S > &x) const
Returns whether two Fad objects have the same values.
KOKKOS_INLINE_FUNCTION ~GeneralFad()
Destructor.
void f()
expr expr dx(i)
KOKKOS_INLINE_FUNCTION void diff(const int ith, const int n)
Set GeneralFad object as the ith independent variable.
#define SACADO_ENABLE_VALUE_CTOR_DECL
KOKKOS_INLINE_FUNCTION GeneralFad(const int sz, const T &x, const DerivInit zero_out=InitDerivArray)
Constructor with size sz and value x.
KOKKOS_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
KOKKOS_INLINE_FUNCTION GeneralFad(const Storage &s)
Constructor with supplied storage s.
KOKKOS_INLINE_FUNCTION GeneralFad(const GeneralFad &x)
Copy constructor.
RemoveConst< T >::type value_type
Typename of values.
#define SACADO_ENABLE_EXPR_CTOR_DECL
KOKKOS_INLINE_FUNCTION void operator()(ArgT arg) const
expr val()
#define KOKKOS_INLINE_FUNCTION
#define T
Definition: Sacado_rad.hpp:573
KOKKOS_INLINE_FUNCTION SACADO_ENABLE_VALUE_FUNC(GeneralFad &) operator
Assignment operator with constant right-hand-side.
KOKKOS_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
KOKKOS_INLINE_FUNCTION void cache() const
Cache values.
Base template specification for testing equivalence.
Wrapper for a generic expression template.
KOKKOS_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
Do not initialize the derivative array.
const int N
KOKKOS_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
KOKKOS_INLINE_FUNCTION GeneralFad()
Default constructor.
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors...
KOKKOS_INLINE_FUNCTION GeneralFad(const int sz, const int i, const T &x)
Constructor with size sz, index i, and value x.
std::ostream & operator<<(std::ostream &os, const GeneralFad< T, Storage > &x)
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
Forward-mode AD class templated on the storage for the derivative array.
Initialize the derivative array.
KOKKOS_INLINE_FUNCTION GeneralFad(const S &x, SACADO_ENABLE_VALUE_CTOR_DECL)
Constructor with supplied value x.
KOKKOS_INLINE_FUNCTION GeneralFad(const Expr< S > &x, SACADO_ENABLE_EXPR_CTOR_DECL)
Copy constructor from any Expression object.
KOKKOS_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.
KOKKOS_INLINE_FUNCTION int availableSize() const
Returns number of derivative components that can be stored without reallocation.