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