Sacado Package Browser (Single Doxygen Collection)  Version of the Day
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) {
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.)) {
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 
219  template <typename S>
221  SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr<S>& x) const {
222  typedef IsEqual<value_type> IE;
223  if (x.size() != this->size()) return false;
224  bool eq = IE::eval(x.val(), this->val());
225  for (int i=0; i<this->size(); i++)
226  eq = eq && IE::eval(x.dx(i), this->dx(i));
227  return eq;
228  }
229 
231 
236 
242  int availableSize() const { return this->length(); }
243 
246  bool hasFastAccess() const { return this->size()!=0;}
247 
250  bool isPassive() const { return this->size()==0;}
251 
254  void setIsConstant(bool is_const) {
255  if (is_const && this->size()!=0)
256  this->resize(0);
257  }
258 
260 
265 
267  template <typename S>
269  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator=(const S& v) {
270  this->val() = v;
271  if (this->size()) this->resize(0);
272  return *this;
273  }
274 
277  GeneralFad&
278  operator=(const GeneralFad& x) {
279  // Copy value & dx_
280  Storage::operator=(x);
281  return *this;
282  }
283 
285  template <typename S>
287  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator=(const Expr<S>& x) {
288  const int xsz = x.size();
289  if (xsz != this->size())
290  this->resizeAndZero(xsz);
291 
292  const int sz = this->size();
293 
294  // For ViewStorage, the resize above may not in fact resize the
295  // derivative array, so it is possible that sz != xsz at this point.
296  // The only valid use case here is sz > xsz == 0, so we use sz in the
297  // assignment below
298 
299  if (sz) {
300 
301  if (Expr<S>::is_linear) {
302  if (x.hasFastAccess())
303  for(int i=0; i<sz; ++i)
304  this->fastAccessDx(i) = x.fastAccessDx(i);
305  else
306  for(int i=0; i<sz; ++i)
307  this->fastAccessDx(i) = x.dx(i);
308  }
309  else {
310 
311  if (x.hasFastAccess()) {
312  // Compute partials
313  FastLocalAccumOp< Expr<S> > op(x);
314 
315  // Compute each tangent direction
316  for(op.i=0; op.i<sz; ++op.i) {
317  op.t = T(0.);
318 
319  // Automatically unrolled loop that computes
320  // for (int j=0; j<N; j++)
321  // op.t += op.partials[j] * x.getTangent<j>(i);
323 
324  this->fastAccessDx(op.i) = op.t;
325  }
326  }
327  else {
328  // Compute partials
329  SlowLocalAccumOp< Expr<S> > op(x);
330 
331  // Compute each tangent direction
332  for(op.i=0; op.i<sz; ++op.i) {
333  op.t = T(0.);
334 
335  // Automatically unrolled loop that computes
336  // for (int j=0; j<N; j++)
337  // op.t += op.partials[j] * x.getTangent<j>(i);
339 
340  this->fastAccessDx(op.i) = op.t;
341  }
342  }
343  }
344  }
345 
346  // Compute value
347  this->val() = x.val();
348 
349  return *this;
350  }
351 
353 
358 
360  template <typename S>
362  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator += (const S& v) {
363  this->val() += v;
364  return *this;
365  }
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  const int sz = this->size();
380  this->val() *= v;
381  for (int i=0; i<sz; ++i)
382  this->fastAccessDx(i) *= v;
383  return *this;
384  }
385 
387  template <typename S>
389  SACADO_ENABLE_VALUE_FUNC(GeneralFad&) operator /= (const S& v) {
390  const int sz = this->size();
391  this->val() /= v;
392  for (int i=0; i<sz; ++i)
393  this->fastAccessDx(i) /= v;
394  return *this;
395  }
396 
399  GeneralFad& operator += (const GeneralFad& x) {
400  const int xsz = x.size(), sz = this->size();
401 
402 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
403  if ((xsz != sz) && (xsz != 0) && (sz != 0))
404  throw "Fad Error: Attempt to assign with incompatible sizes";
405 #endif
406 
407  if (xsz) {
408  if (sz) {
409  for (int i=0; i<sz; ++i)
410  this->fastAccessDx(i) += x.fastAccessDx(i);
411  }
412  else {
413  this->resizeAndZero(xsz);
414  for (int i=0; i<xsz; ++i)
415  this->fastAccessDx(i) = x.fastAccessDx(i);
416  }
417  }
418 
419  this->val() += x.val();
420 
421  return *this;
422  }
423 
426  GeneralFad& operator -= (const GeneralFad& x) {
427  const int xsz = x.size(), sz = this->size();
428 
429 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
430  if ((xsz != sz) && (xsz != 0) && (sz != 0))
431  throw "Fad Error: Attempt to assign with incompatible sizes";
432 #endif
433 
434  if (xsz) {
435  if (sz) {
436  for(int i=0; i<sz; ++i)
437  this->fastAccessDx(i) -= x.fastAccessDx(i);
438  }
439  else {
440  this->resizeAndZero(xsz);
441  for(int i=0; i<xsz; ++i)
442  this->fastAccessDx(i) = -x.fastAccessDx(i);
443  }
444  }
445 
446  this->val() -= x.val();
447 
448 
449  return *this;
450  }
451 
454  GeneralFad& operator *= (const GeneralFad& x) {
455  const int xsz = x.size(), sz = this->size();
456  T xval = x.val();
457  T v = this->val();
458 
459 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
460  if ((xsz != sz) && (xsz != 0) && (sz != 0))
461  throw "Fad Error: Attempt to assign with incompatible sizes";
462 #endif
463 
464  if (xsz) {
465  if (sz) {
466  for(int i=0; i<sz; ++i)
467  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
468  }
469  else {
470  this->resizeAndZero(xsz);
471  for(int i=0; i<xsz; ++i)
472  this->fastAccessDx(i) = v*x.fastAccessDx(i);
473  }
474  }
475  else {
476  if (sz) {
477  for (int i=0; i<sz; ++i)
478  this->fastAccessDx(i) *= xval;
479  }
480  }
481 
482  this->val() *= xval;
483 
484  return *this;
485  }
486 
489  GeneralFad& operator /= (const GeneralFad& x) {
490  const int xsz = x.size(), sz = this->size();
491  T xval = x.val();
492  T v = this->val();
493 
494 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
495  if ((xsz != sz) && (xsz != 0) && (sz != 0))
496  throw "Fad Error: Attempt to assign with incompatible sizes";
497 #endif
498 
499  if (xsz) {
500  if (sz) {
501  for(int i=0; i<sz; ++i)
502  this->fastAccessDx(i) =
503  ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
504  }
505  else {
506  this->resizeAndZero(xsz);
507  for(int i=0; i<xsz; ++i)
508  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
509  }
510  }
511  else {
512  if (sz) {
513  for (int i=0; i<sz; ++i)
514  this->fastAccessDx(i) /= xval;
515  }
516  }
517 
518  this->val() /= xval;
519 
520  return *this;
521  }
522 
524  template <typename S>
526  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator += (const Expr<S>& x) {
527  const int xsz = x.size(), sz = this->size();
528 
529 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
530  if ((xsz != sz) && (xsz != 0) && (sz != 0))
531  throw "Fad Error: Attempt to assign with incompatible sizes";
532 #endif
533 
534  if (Expr<S>::is_linear) {
535  if (xsz) {
536  if (sz) {
537  if (x.hasFastAccess())
538  for (int i=0; i<sz; ++i)
539  this->fastAccessDx(i) += x.fastAccessDx(i);
540  else
541  for (int i=0; i<sz; ++i)
542  this->fastAccessDx(i) += x.dx(i);
543  }
544  else {
545  this->resizeAndZero(xsz);
546  if (x.hasFastAccess())
547  for (int i=0; i<xsz; ++i)
548  this->fastAccessDx(i) = x.fastAccessDx(i);
549  else
550  for (int i=0; i<xsz; ++i)
551  this->fastAccessDx(i) = x.dx(i);
552  }
553  }
554  }
555  else {
556 
557  if (xsz) {
558 
559  if (sz != xsz)
560  this->resizeAndZero(xsz);
561 
562  if (x.hasFastAccess()) {
563  // Compute partials
564  FastLocalAccumOp< Expr<S> > op(x);
565 
566  // Compute each tangent direction
567  for(op.i=0; op.i<xsz; ++op.i) {
568  op.t = T(0.);
569 
570  // Automatically unrolled loop that computes
571  // for (int j=0; j<N; j++)
572  // op.t += op.partials[j] * x.getTangent<j>(i);
574 
575  this->fastAccessDx(op.i) += op.t;
576  }
577  }
578  else {
579  // Compute partials
580  SlowLocalAccumOp< Expr<S> > op(x);
581 
582  // Compute each tangent direction
583  for(op.i=0; op.i<xsz; ++op.i) {
584  op.t = T(0.);
585 
586  // Automatically unrolled loop that computes
587  // for (int j=0; j<N; j++)
588  // op.t += op.partials[j] * x.getTangent<j>(i);
590 
591  this->fastAccessDx(op.i) += op.t;
592  }
593  }
594 
595  }
596 
597  }
598 
599  // Compute value
600  this->val() += x.val();
601 
602  return *this;
603  }
604 
606  template <typename S>
608  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator -= (const Expr<S>& x) {
609  const int xsz = x.size(), sz = this->size();
610 
611 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
612  if ((xsz != sz) && (xsz != 0) && (sz != 0))
613  throw "Fad Error: Attempt to assign with incompatible sizes";
614 #endif
615 
616  if (Expr<S>::is_linear) {
617  if (xsz) {
618  if (sz) {
619  if (x.hasFastAccess())
620  for(int i=0; i<sz; ++i)
621  this->fastAccessDx(i) -= x.fastAccessDx(i);
622  else
623  for (int i=0; i<sz; ++i)
624  this->fastAccessDx(i) -= x.dx(i);
625  }
626  else {
627  this->resizeAndZero(xsz);
628  if (x.hasFastAccess())
629  for(int i=0; i<xsz; ++i)
630  this->fastAccessDx(i) = -x.fastAccessDx(i);
631  else
632  for (int i=0; i<xsz; ++i)
633  this->fastAccessDx(i) = -x.dx(i);
634  }
635  }
636  }
637  else {
638 
639  if (xsz) {
640 
641  if (sz != xsz)
642  this->resizeAndZero(xsz);
643 
644  if (x.hasFastAccess()) {
645  // Compute partials
646  FastLocalAccumOp< Expr<S> > op(x);
647 
648  // Compute each tangent direction
649  for(op.i=0; op.i<xsz; ++op.i) {
650  op.t = T(0.);
651 
652  // Automatically unrolled loop that computes
653  // for (int j=0; j<N; j++)
654  // op.t += op.partials[j] * x.getTangent<j>(i);
656 
657  this->fastAccessDx(op.i) -= op.t;
658  }
659  }
660  else {
661  // Compute partials
662  SlowLocalAccumOp< Expr<S> > op(x);
663 
664  // Compute each tangent direction
665  for(op.i=0; op.i<xsz; ++op.i) {
666  op.t = T(0.);
667 
668  // Automatically unrolled loop that computes
669  // for (int j=0; j<N; j++)
670  // op.t += op.partials[j] * x.getTangent<j>(i);
672 
673  this->fastAccessDx(op.i) -= op.t;
674  }
675  }
676  }
677 
678  }
679 
680  this->val() -= x.val();
681 
682  return *this;
683  }
684 
686  template <typename S>
688  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator *= (const Expr<S>& x) {
689  const int xsz = x.size(), sz = this->size();
690  T xval = x.val();
691  T v = this->val();
692 
693 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
694  if ((xsz != sz) && (xsz != 0) && (sz != 0))
695  throw "Fad Error: Attempt to assign with incompatible sizes";
696 #endif
697 
698  if (Expr<S>::is_linear) {
699  if (xsz) {
700  if (sz) {
701  if (x.hasFastAccess())
702  for(int i=0; i<sz; ++i)
703  this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
704  else
705  for (int i=0; i<sz; ++i)
706  this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
707  }
708  else {
709  this->resizeAndZero(xsz);
710  if (x.hasFastAccess())
711  for(int i=0; i<xsz; ++i)
712  this->fastAccessDx(i) = v*x.fastAccessDx(i);
713  else
714  for (int i=0; i<xsz; ++i)
715  this->fastAccessDx(i) = v*x.dx(i);
716  }
717  }
718  else {
719  if (sz) {
720  for (int i=0; i<sz; ++i)
721  this->fastAccessDx(i) *= xval;
722  }
723  }
724  }
725  else {
726 
727  if (xsz) {
728 
729  if (sz) {
730 
731  if (x.hasFastAccess()) {
732  // Compute partials
733  FastLocalAccumOp< Expr<S> > op(x);
734 
735  // Compute each tangent direction
736  for(op.i=0; op.i<xsz; ++op.i) {
737  op.t = T(0.);
738 
739  // Automatically unrolled loop that computes
740  // for (int j=0; j<N; j++)
741  // op.t += op.partials[j] * x.getTangent<j>(i);
743 
744  this->fastAccessDx(op.i) =
745  v * op.t + this->fastAccessDx(op.i) * xval;
746  }
747  }
748  else {
749  // Compute partials
750  SlowLocalAccumOp< Expr<S> > op(x);
751 
752  // Compute each tangent direction
753  for(op.i=0; op.i<xsz; ++op.i) {
754  op.t = T(0.);
755 
756  // Automatically unrolled loop that computes
757  // for (int j=0; j<N; j++)
758  // op.t += op.partials[j] * x.getTangent<j>(i);
760 
761  this->fastAccessDx(op.i) =
762  v * op.t + this->fastAccessDx(op.i) * xval;
763  }
764  }
765 
766  }
767 
768  else {
769 
770  this->resizeAndZero(xsz);
771 
772  if (x.hasFastAccess()) {
773  // Compute partials
774  FastLocalAccumOp< Expr<S> > op(x);
775 
776  // Compute each tangent direction
777  for(op.i=0; op.i<xsz; ++op.i) {
778  op.t = T(0.);
779 
780  // Automatically unrolled loop that computes
781  // for (int j=0; j<N; j++)
782  // op.t += op.partials[j] * x.getTangent<j>(i);
784 
785  this->fastAccessDx(op.i) = v * op.t;
786  }
787  }
788  else {
789  // Compute partials
790  SlowLocalAccumOp< Expr<S> > op(x);
791 
792  // Compute each tangent direction
793  for(op.i=0; op.i<xsz; ++op.i) {
794  op.t = T(0.);
795 
796  // Automatically unrolled loop that computes
797  // for (int j=0; j<N; j++)
798  // op.t += op.partials[j] * x.getTangent<j>(i);
800 
801  this->fastAccessDx(op.i) = v * op.t;
802  }
803  }
804 
805  }
806 
807  }
808 
809  else {
810 
811  if (sz) {
812  for (int i=0; i<sz; ++i)
813  this->fastAccessDx(i) *= xval;
814  }
815 
816  }
817 
818  }
819 
820  this->val() *= xval;
821 
822  return *this;
823  }
824 
826  template <typename S>
828  SACADO_ENABLE_EXPR_FUNC(GeneralFad&) operator /= (const Expr<S>& x) {
829  const int xsz = x.size(), sz = this->size();
830  T xval = x.val();
831  T v = this->val();
832 
833 #if defined(SACADO_DEBUG) && !defined(__CUDA_ARCH__ )
834  if ((xsz != sz) && (xsz != 0) && (sz != 0))
835  throw "Fad Error: Attempt to assign with incompatible sizes";
836 #endif
837 
838  if (Expr<S>::is_linear) {
839  if (xsz) {
840  if (sz) {
841  if (x.hasFastAccess())
842  for(int i=0; i<sz; ++i)
843  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
844  else
845  for (int i=0; i<sz; ++i)
846  this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
847  }
848  else {
849  this->resizeAndZero(xsz);
850  if (x.hasFastAccess())
851  for(int i=0; i<xsz; ++i)
852  this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
853  else
854  for (int i=0; i<xsz; ++i)
855  this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
856  }
857  }
858  else {
859  if (sz) {
860  for (int i=0; i<sz; ++i)
861  this->fastAccessDx(i) /= xval;
862  }
863  }
864  }
865  else {
866 
867  if (xsz) {
868 
869  T xval2 = xval*xval;
870 
871  if (sz) {
872 
873  if (x.hasFastAccess()) {
874  // Compute partials
875  FastLocalAccumOp< Expr<S> > op(x);
876 
877  // Compute each tangent direction
878  for(op.i=0; op.i<xsz; ++op.i) {
879  op.t = T(0.);
880 
881  // Automatically unrolled loop that computes
882  // for (int j=0; j<N; j++)
883  // op.t += op.partials[j] * x.getTangent<j>(i);
885 
886  this->fastAccessDx(op.i) =
887  (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
888  }
889  }
890  else {
891  // Compute partials
892  SlowLocalAccumOp< Expr<S> > op(x);
893 
894  // Compute each tangent direction
895  for(op.i=0; op.i<xsz; ++op.i) {
896  op.t = T(0.);
897 
898  // Automatically unrolled loop that computes
899  // for (int j=0; j<N; j++)
900  // op.t += op.partials[j] * x.getTangent<j>(i);
902 
903  this->fastAccessDx(op.i) =
904  (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
905  }
906  }
907 
908  }
909 
910  else {
911 
912  this->resizeAndZero(xsz);
913 
914  if (x.hasFastAccess()) {
915  // Compute partials
916  FastLocalAccumOp< Expr<S> > op(x);
917 
918  // Compute each tangent direction
919  for(op.i=0; op.i<xsz; ++op.i) {
920  op.t = T(0.);
921 
922  // Automatically unrolled loop that computes
923  // for (int j=0; j<N; j++)
924  // op.t += op.partials[j] * x.getTangent<j>(i);
926 
927  this->fastAccessDx(op.i) = (-v * op.t) / xval2;
928  }
929  }
930  else {
931  // Compute partials
932  SlowLocalAccumOp< 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 
947  }
948 
949  }
950 
951  else {
952 
953  if (sz) {
954  for (int i=0; i<sz; ++i)
955  this->fastAccessDx(i) /= xval;
956  }
957 
958  }
959 
960  }
961 
962  this->val() /= xval;
963 
964  return *this;
965  }
966 
968 
969  protected:
970 
971  // Functor for mpl::for_each to compute the local accumulation
972  // of a tangent derivative
973  //
974  // We use getTangent<>() to get dx components from expression
975  // arguments instead of getting the argument directly or extracting
976  // the dx array due to striding in ViewFad (or could use striding
977  // directly here if we need to get dx array).
978  template <typename ExprT>
979  struct FastLocalAccumOp {
980  typedef typename ExprT::value_type value_type;
981  static const int N = ExprT::num_args;
982  const ExprT& x;
983  mutable value_type t;
984  value_type partials[N];
985  int i;
987  FastLocalAccumOp(const ExprT& x_) : x(x_) {
988  x.computePartials(value_type(1.), partials);
989  }
990  template <typename ArgT>
992  void operator () (ArgT arg) const {
993  const int Arg = ArgT::value;
994  t += partials[Arg] * x.template getTangent<Arg>(i);
995  }
996  };
997 
998  template <typename ExprT>
999  struct SlowLocalAccumOp : FastLocalAccumOp<ExprT> {
1001  SlowLocalAccumOp(const ExprT& x_) :
1002  FastLocalAccumOp<ExprT>(x_) {}
1003  template <typename ArgT>
1005  void operator () (ArgT arg) const {
1006  const int Arg = ArgT::value;
1007  if (this->x.template isActive<Arg>())
1008  this->t += this->partials[Arg] * this->x.template getTangent<Arg>(this->i);
1009  }
1010  };
1011 
1012  }; // class GeneralFad
1013 
1014 
1015  template <typename T, typename Storage>
1016  std::ostream& operator << (std::ostream& os,
1017  const GeneralFad<T,Storage>& x) {
1018  os << x.val() << " [";
1019 
1020  for (int i=0; i< x.size(); i++) {
1021  os << " " << x.dx(i);
1022  }
1023 
1024  os << " ]";
1025  return os;
1026  }
1027 
1028  } // namespace ELRFad
1029 
1030 } // namespace Sacado
1031 
1032 #endif // SACADO_ELRFAD_GENERALFAD_HPP
RemoveConst< T >::type value_type
Typename of values.
void f()
expr val()
KOKKOS_INLINE_FUNCTION GeneralFad()
Default constructor.
KOKKOS_INLINE_FUNCTION SACADO_ENABLE_EXPR_FUNC(bool) isEqualTo(const Expr< S > &x) const
Returns whether two Fad objects have the same values.
#define SACADO_ENABLE_VALUE_CTOR_DECL
KOKKOS_INLINE_FUNCTION GeneralFad(const Storage &s)
Constructor with supplied storage s.
KOKKOS_INLINE_FUNCTION bool isPassive() const
Returns true if derivative array is empty.
#define SACADO_ENABLE_EXPR_CTOR_DECL
KOKKOS_INLINE_FUNCTION void setUpdateValue(bool update_val)
Set whether this Fad object should update values.
#define KOKKOS_INLINE_FUNCTION
#define T
Definition: Sacado_rad.hpp:573
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
KOKKOS_INLINE_FUNCTION SlowLocalAccumOp(const ExprT &x_)
KOKKOS_INLINE_FUNCTION void operator()(ArgT arg) const
Base template specification for testing equivalence.
KOKKOS_INLINE_FUNCTION bool updateValue() const
Return whether this Fad object has an updated value.
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.
DerivInit
Enum use to signal whether the derivative array should be initialized in AD object constructors...
Forward-mode AD class templated on the storage for the derivative array.
Wrapper for a generic expression template.
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.
expr expr dx(i)
ScalarType< value_type >::type scalar_type
Typename of scalar&#39;s (which may be different from T)
int n
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 int availableSize() const
Returns number of derivative components that can be stored without reallocation.
KOKKOS_INLINE_FUNCTION void setIsConstant(bool is_const)
Set whether variable is constant.
KOKKOS_INLINE_FUNCTION bool hasFastAccess() const
Returns true if derivative array is not empty.