Reference documentation for deal.II version 9.0.0
sacado_number_types.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2017 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_differentiation_ad_sacado_number_types_h
17 #define dealii_differentiation_ad_sacado_number_types_h
18 
19 #include <deal.II/base/config.h>
20 #include <type_traits>
21 
22 
23 DEAL_II_NAMESPACE_OPEN
24 
25 
26 namespace Differentiation
27 {
28  namespace AD
29  {
30 
38  template <typename NumberType, typename = void>
40  : std::false_type
41  {};
42 
43 
51  template <typename NumberType, typename = void>
53  : std::false_type
54  {};
55 
56 
64  template <typename NumberType, typename = void>
66  : std::false_type
67  {};
68 
69  } // namespace AD
70 } // namespace Differentiation
71 
72 
73 DEAL_II_NAMESPACE_CLOSE
74 
75 
76 
77 #ifdef DEAL_II_WITH_TRILINOS
78 
79 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
80 #include <Sacado.hpp>
81 // It appears that some versions of Trilinos do not directly or indirectly
82 // include all the headers for all forward and reverse Sacado AD types.
83 // So we directly include these both here as a precaution.
84 // Standard forward AD classes (templated)
85 #include <Sacado_Fad_DFad.hpp>
86 // Reverse AD classes (templated)
87 #include <Sacado_trad.hpp>
88 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
89 
90 #include <deal.II/base/numbers.h>
91 #include <deal.II/base/exceptions.h>
92 
93 #include <deal.II/differentiation/ad/ad_number_types.h>
94 #include <deal.II/differentiation/ad/ad_number_traits.h>
95 
96 #include <complex>
97 #include <type_traits>
98 
99 DEAL_II_NAMESPACE_OPEN
100 
101 
102 namespace Differentiation
103 {
104  namespace AD
105  {
106 
107 
108  namespace internal
109  {
110 
119  template <typename SacadoNumber, typename = void>
121 
122  }
123 
124 
125 
126  } // namespace AD
127 } // namespace Differentiation
128 
129 
130 /* --------------------------- inline and template functions and specializations ------------------------- */
131 
132 
133 #ifndef DOXYGEN
134 
135 namespace Differentiation
136 {
137  namespace AD
138  {
139 
140  namespace internal
141  {
142 
143  // The documentation on Sacado numbers is pretty sparse and/or hard to
144  // navigate. As a point of reference, see
145  // https://trilinos.org/docs/dev/packages/sacado/doc/html/classSacado_1_1Fad_1_1SimpleFad.html
146  // for semi-applicable documentation for the Sacado::Fad::Dfad class.
147  // and the examples in
148  // https://github.com/trilinos/Trilinos/tree/master/packages/sacado/example
149  //
150  // If one dares to venture there, the relevant files for the classes supported here are:
151  //
152  // Forward-mode auto-differentiable types:
153  // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_DFad.hpp
154  // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/sacado_dfad_GeneralFad.hpp
155  //
156  // Reverse-mode auto-differentiable types:
157  // https://github.com/trilinos/Trilinos/blob/master/packages/sacado/src/Sacado_trad.hpp
158 
159 
163  template <typename SacadoNumber>
164  struct SacadoNumberInfo<SacadoNumber, typename std::enable_if<
165  std::is_same<SacadoNumber,Sacado::Fad::DFad<typename SacadoNumber::value_type> >::value
166  >::type>
167  {
168  typedef SacadoNumber ad_type;
169  typedef typename ad_type::scalar_type scalar_type;
170  typedef typename ad_type::value_type value_type;
171  typedef typename ad_type::value_type derivative_type;
172 
173  static const unsigned int n_supported_derivative_levels
174  = 1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
175  };
176 
177 
181  template <typename SacadoNumber>
182  struct SacadoNumberInfo<SacadoNumber, typename std::enable_if<
183  std::is_same<SacadoNumber,Sacado::Rad::ADvar<typename SacadoNumber::value_type> >::value
184  >::type>
185  {
186  typedef SacadoNumber ad_type;
187  typedef typename ad_type::ADVari::scalar_type scalar_type;
188  typedef typename ad_type::ADVari::value_type value_type;
189  typedef typename ad_type::ADVari::value_type derivative_type;
190 
191  static const unsigned int n_supported_derivative_levels
192  = 1 + SacadoNumberInfo<derivative_type>::n_supported_derivative_levels;
193  };
194 
195 
202  template <typename Number>
203  struct SacadoNumberInfo<Number, typename std::enable_if<
204  std::is_arithmetic<typename std::decay<Number>::type>::value
205  >::type>
206  {
207  static const unsigned int n_supported_derivative_levels = 0;
208  };
209 
210 
215  template<typename ScalarType>
216  struct ADNumberInfoFromEnum<ScalarType, Differentiation::AD::NumberTypes::sacado_dfad,
217  typename std::enable_if<std::is_floating_point<ScalarType>::value>::type
218  >
219  {
220  static const bool is_taped = false;
221  typedef Sacado::Fad::DFad<ScalarType> real_type;
222  typedef typename SacadoNumberInfo<real_type>::derivative_type derivative_type;
223  static const unsigned int n_supported_derivative_levels
224  = SacadoNumberInfo<real_type>::n_supported_derivative_levels;
225  };
226 
227 
232  template<typename ScalarType>
233  struct ADNumberInfoFromEnum<ScalarType, Differentiation::AD::NumberTypes::sacado_dfad_dfad,
234  typename std::enable_if<std::is_floating_point<ScalarType>::value>::type
235  >
236  {
237  static const bool is_taped = false;
238  typedef Sacado::Fad::DFad< Sacado::Fad::DFad<ScalarType> > real_type;
239  typedef typename SacadoNumberInfo<real_type>::derivative_type derivative_type;
240  static const unsigned int n_supported_derivative_levels
241  = SacadoNumberInfo<real_type>::n_supported_derivative_levels;
242  };
243 
244 
249  template<typename ScalarType>
250  struct ADNumberInfoFromEnum<ScalarType, Differentiation::AD::NumberTypes::sacado_rad,
251  typename std::enable_if<std::is_floating_point<ScalarType>::value>::type
252  >
253  {
254  static const bool is_taped = false;
255  typedef Sacado::Rad::ADvar<ScalarType> real_type;
256  typedef typename SacadoNumberInfo<real_type>::derivative_type derivative_type;
257  static const unsigned int n_supported_derivative_levels
258  = SacadoNumberInfo<real_type>::n_supported_derivative_levels;
259  };
260 
261 
266  template<typename ScalarType>
267  struct ADNumberInfoFromEnum<ScalarType, Differentiation::AD::NumberTypes::sacado_rad_dfad,
268  typename std::enable_if<std::is_floating_point<ScalarType>::value>::type
269  >
270  {
271  static const bool is_taped = false;
272  typedef Sacado::Rad::ADvar< Sacado::Fad::DFad<ScalarType> > real_type;
273  typedef typename SacadoNumberInfo<real_type>::derivative_type derivative_type;
274  static const unsigned int n_supported_derivative_levels
275  = SacadoNumberInfo<real_type>::n_supported_derivative_levels;
276  };
277 
278 
283  template<typename NumberType>
284  struct Marking< Sacado::Fad::DFad<NumberType> >
285  {
286  typedef typename SacadoNumberInfo< Sacado::Fad::DFad<NumberType> >::ad_type ad_type;
287  typedef typename SacadoNumberInfo< Sacado::Fad::DFad<NumberType> >::derivative_type derivative_type;
288  typedef typename SacadoNumberInfo< Sacado::Fad::DFad<NumberType> >::scalar_type scalar_type;
289 
290  /*
291  * Initialize the state of an independent variable.
292  */
293  static void
294  independent_variable(const scalar_type &in,
295  const unsigned int index,
296  const unsigned int n_independent_variables,
297  ad_type &out)
298  {
299  // It is required that we first initialise the outer number before
300  // any of the nested ones.
301  out = ad_type(n_independent_variables, index, in);
302 
303  // Initialize potential nested directional derivatives
304  Marking<derivative_type>::independent_variable(
305  in, index, n_independent_variables, out.val());
306  }
307 
308  /*
309  * Initialize the state of a dependent variable.
310  */
311  static void
312  dependent_variable(ad_type &out,
313  const ad_type &func)
314  {
315  out = func;
316  }
317  };
318 
319 
324  template<typename NumberType>
325  struct Marking< Sacado::Rad::ADvar<NumberType> >
326  {
327  typedef typename SacadoNumberInfo< Sacado::Rad::ADvar<NumberType> >::ad_type ad_type;
328  typedef typename SacadoNumberInfo< Sacado::Rad::ADvar<NumberType> >::derivative_type derivative_type;
329  typedef typename SacadoNumberInfo< Sacado::Rad::ADvar<NumberType> >::scalar_type scalar_type;
330 
331  /*
332  * Initialize the state of an independent variable.
333  */
334  static void
335  independent_variable(const scalar_type &in,
336  const unsigned int index,
337  const unsigned int n_independent_variables,
338  ad_type &out)
339  {
340  // For Sacado::Rad::ADvar numbers, we have to initialize the
341  // ADNumber with an already fully-configured value. This means
342  // that if this nests another ADNumber then the nested number
343  // must already be setup and ready for use.
344 
345  // Initialize potential nested directional derivatives
346  derivative_type derivative_initializer;
347  Marking<derivative_type>::independent_variable(
348  in, index, n_independent_variables,
349  derivative_initializer);
350 
351  // Initialize the outer ad_type
352  out = derivative_initializer;
353  }
354 
355  /*
356  * Initialize the state of a dependent variable.
357  */
358  static void
359  dependent_variable(ad_type &out,
360  const ad_type &func)
361  {
362  out = func;
363  }
364  };
365 
366 
374  template<typename NumberType>
375  struct ExtractData< Sacado::Fad::DFad<NumberType> >
376  {
377  typedef typename SacadoNumberInfo< Sacado::Fad::DFad<NumberType> >::derivative_type derivative_type;
378  typedef typename SacadoNumberInfo< Sacado::Fad::DFad<NumberType> >::scalar_type scalar_type;
379  typedef typename SacadoNumberInfo< Sacado::Fad::DFad<NumberType> >::value_type value_type;
380 
384  static scalar_type
385  value (const Sacado::Fad::DFad<NumberType> &x)
386  {
387  return ExtractData<value_type>::value(x.val());
388  }
389 
390 
394  static unsigned int
395  n_directional_derivatives (const Sacado::Fad::DFad<NumberType> &x)
396  {
397  return x.size();
398  }
399 
400 
404  static derivative_type
405  directional_derivative (const Sacado::Fad::DFad<NumberType> &x,
406  const unsigned int direction)
407  {
408  if (x.hasFastAccess())
409  return x.fastAccessDx(direction);
410  else
411  return x.dx(direction);
412  }
413  };
414 
415 
423  template<typename NumberType>
424  struct ExtractData< Sacado::Rad::ADvar<NumberType> >
425  {
426  typedef typename SacadoNumberInfo< Sacado::Rad::ADvar<NumberType> >::derivative_type derivative_type;
427  typedef typename SacadoNumberInfo< Sacado::Rad::ADvar<NumberType> >::scalar_type scalar_type;
428  typedef typename SacadoNumberInfo< Sacado::Rad::ADvar<NumberType> >::value_type value_type;
429 
433  static scalar_type
434  value (const Sacado::Rad::ADvar<NumberType> &x)
435  {
436  return ExtractData<value_type>::value(x.val());
437  }
438 
439 
443  static unsigned int
444  n_directional_derivatives (const Sacado::Rad::ADvar<NumberType> &)
445  {
446  // There are as many directional derivatives as there are
447  // independent variables, but each independent variable can
448  // only return one directional derivative.
449  return 1;
450  }
451 
452 
460  static derivative_type
461  directional_derivative (const Sacado::Rad::ADvar<NumberType> &x,
462  const unsigned int )
463  {
464  return x.adj();
465  }
466  };
467 
468  } // namespace internal
469 
470 
471  /* -------------- NumberTypes::sacado_dfad -------------- */
472 
473 
480  template<typename ADNumberType>
481  struct ADNumberTraits<ADNumberType, typename std::enable_if<
482  std::is_same<ADNumberType, Sacado::Fad::DFad<typename ADNumberType::scalar_type> >::value
483  >::type>
484  : NumberTraits<typename ADNumberType::scalar_type, NumberTypes::sacado_dfad>
485  {};
486 
487 
494  template<typename ADNumberType>
495  struct ADNumberTraits<ADNumberType, typename std::enable_if<
496  std::is_same<ADNumberType, std::complex<Sacado::Fad::DFad<typename ADNumberType::value_type::scalar_type> > >::value
497  >::type>
498  : NumberTraits<std::complex<typename ADNumberType::value_type::scalar_type>, NumberTypes::sacado_dfad>
499  {};
500 
501 
506  template<>
507  struct NumberTraits<Sacado::Fad::DFad<float>,NumberTypes::sacado_dfad>
508  : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<float> >::scalar_type,NumberTypes::sacado_dfad>
509  {};
510 
511 
516  template<>
517  struct NumberTraits<std::complex< Sacado::Fad::DFad<float> >,NumberTypes::sacado_dfad>
518  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Fad::DFad<float> > >::scalar_type,NumberTypes::sacado_dfad>
519  {};
520 
521 
526  template<>
527  struct NumberTraits<Sacado::Fad::DFad<double>,NumberTypes::sacado_dfad>
528  : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<double> >::scalar_type,NumberTypes::sacado_dfad>
529  {};
530 
531 
536  template<>
537  struct NumberTraits<std::complex< Sacado::Fad::DFad<double> >,NumberTypes::sacado_dfad>
538  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Fad::DFad<double> > >::scalar_type,NumberTypes::sacado_dfad>
539  {};
540 
541 
542  /* -------------- NumberTypes::sacado_rad -------------- */
543 
544 
551  template<typename ADNumberType>
552  struct ADNumberTraits<ADNumberType, typename std::enable_if<
553  std::is_same<ADNumberType, Sacado::Rad::ADvar<typename ADNumberType::ADVari::scalar_type> >::value
554  >::type>
555  : NumberTraits<typename ADNumberType::ADVari::scalar_type, NumberTypes::sacado_rad>
556  {};
557 
558 
563  template<>
564  struct NumberTraits<Sacado::Rad::ADvar<float>,NumberTypes::sacado_rad>
565  : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<float> >::scalar_type,NumberTypes::sacado_rad>
566  {};
567 
568 
573  template<>
574  struct NumberTraits<Sacado::Rad::ADvar<double>,NumberTypes::sacado_rad>
575  : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<double> >::scalar_type,NumberTypes::sacado_rad>
576  {};
577 
578 
579 #ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
580 
581 
588  template<typename ADNumberType>
589  struct ADNumberTraits<ADNumberType, typename std::enable_if<
590  std::is_same<ADNumberType, std::complex<Sacado::Rad::ADvar<typename ADNumberType::value_type::ADVari::scalar_type> > >::value
591  >::type>
592  : NumberTraits<std::complex<typename ADNumberType::value_type::ADVari::scalar_type>, NumberTypes::sacado_rad>
593  {};
594 
595 
600  template<>
601  struct NumberTraits<std::complex< Sacado::Rad::ADvar<float> >,NumberTypes::sacado_rad>
602  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Rad::ADvar<float> > >::scalar_type,NumberTypes::sacado_rad>
603  {};
604 
605 
610  template<>
611  struct NumberTraits<std::complex< Sacado::Rad::ADvar<double> >,NumberTypes::sacado_rad>
612  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Rad::ADvar<double> > >::scalar_type,NumberTypes::sacado_rad>
613  {};
614 
615 
616 #endif
617 
618 
619  /* -------------- NumberTypes::sacado_dfad_dfad -------------- */
620 
628  template<typename ADNumberType>
629  struct ADNumberTraits<ADNumberType, typename std::enable_if<
630  std::is_same<ADNumberType, Sacado::Fad::DFad<Sacado::Fad::DFad<typename ADNumberType::scalar_type> > >::value
631  >::type>
632  : NumberTraits<typename ADNumberType::scalar_type, NumberTypes::sacado_dfad_dfad>
633  {};
634 
635 
643  template<typename ADNumberType>
644  struct ADNumberTraits<ADNumberType, typename std::enable_if<
645  std::is_same<ADNumberType, std::complex<Sacado::Fad::DFad<Sacado::Fad::DFad<typename ADNumberType::value_type::scalar_type> > > >::value
646  >::type>
647  : NumberTraits<std::complex<typename ADNumberType::value_type::scalar_type>, NumberTypes::sacado_dfad_dfad>
648  {};
649 
650 
655  template<>
656  struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<float> >,NumberTypes::sacado_dfad_dfad>
657  : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<float> > >::scalar_type,NumberTypes::sacado_dfad_dfad>
658  {};
659 
660 
665  template<>
666  struct NumberTraits<std::complex< Sacado::Fad::DFad<Sacado::Fad::DFad<float> > >,NumberTypes::sacado_dfad_dfad>
667  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Fad::DFad<Sacado::Fad::DFad<float> > > >::scalar_type,NumberTypes::sacado_dfad_dfad>
668  {};
669 
670 
675  template<>
676  struct NumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<double> >,NumberTypes::sacado_dfad_dfad>
677  : NumberTraits<typename ADNumberTraits<Sacado::Fad::DFad<Sacado::Fad::DFad<double> > >::scalar_type,NumberTypes::sacado_dfad_dfad>
678  {};
679 
680 
685  template<>
686  struct NumberTraits<std::complex< Sacado::Fad::DFad<Sacado::Fad::DFad<double> > >,NumberTypes::sacado_dfad_dfad>
687  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Fad::DFad<Sacado::Fad::DFad<double> > > >::scalar_type,NumberTypes::sacado_dfad_dfad>
688  {};
689 
690 
691  /* -------------- NumberTypes::sacado_rad_dfad -------------- */
692 
700  template<typename ADNumberType>
701  struct ADNumberTraits<ADNumberType, typename std::enable_if<
702  std::is_same<ADNumberType, Sacado::Rad::ADvar<Sacado::Fad::DFad<typename ADNumberType::ADVari::scalar_type> > >::value
703  >::type>
704  : NumberTraits<typename ADNumberType::ADVari::scalar_type, NumberTypes::sacado_rad_dfad>
705  {};
706 
707 
712  template<>
713  struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<float> >,NumberTypes::sacado_rad_dfad>
714  : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<float> > >::scalar_type,NumberTypes::sacado_rad_dfad>
715  {};
716 
717 
722  template<>
723  struct NumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<double> >,NumberTypes::sacado_rad_dfad>
724  : NumberTraits<typename ADNumberTraits<Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > >::scalar_type,NumberTypes::sacado_rad_dfad>
725  {};
726 
727 
728 #ifdef DEAL_II_TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD
729 
730 
738  template<typename ADNumberType>
739  struct ADNumberTraits<ADNumberType, typename std::enable_if<
740  std::is_same<ADNumberType, std::complex<Sacado::Rad::ADvar<Sacado::Fad::DFad<typename ADNumberType::value_type::ADVari::scalar_type> > > >::value
741  >::type>
742  : NumberTraits<std::complex<typename ADNumberType::value_type::ADVari::scalar_type>, NumberTypes::sacado_rad_dfad>
743  {};
744 
745 
750  template<>
751  struct NumberTraits<std::complex< Sacado::Rad::ADvar<Sacado::Fad::DFad<float> > >,NumberTypes::sacado_rad_dfad>
752  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Rad::ADvar<Sacado::Fad::DFad<float> > > >::scalar_type,NumberTypes::sacado_rad_dfad>
753  {};
754 
755 
760  template<>
761  struct NumberTraits<std::complex< Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > >,NumberTypes::sacado_rad_dfad>
762  : NumberTraits<typename ADNumberTraits<std::complex< Sacado::Rad::ADvar<Sacado::Fad::DFad<double> > > >::scalar_type,NumberTypes::sacado_rad_dfad>
763  {};
764 
765 
766 #endif
767 
768 
769  /* -------------- Additional type traits -------------- */
770 
771 
772  template <typename NumberType>
773  struct is_sacado_dfad_number<NumberType, typename std::enable_if<
774  ADNumberTraits<typename std::decay<NumberType>::type>::type_code == NumberTypes::sacado_dfad ||
775  ADNumberTraits<typename std::decay<NumberType>::type>::type_code == NumberTypes::sacado_dfad_dfad
776  >::type>
777  : std::true_type
778  {};
779 
780 
781  template <typename NumberType>
782  struct is_sacado_rad_number<NumberType, typename std::enable_if<
783  ADNumberTraits<typename std::decay<NumberType>::type>::type_code == NumberTypes::sacado_rad ||
784  ADNumberTraits<typename std::decay<NumberType>::type>::type_code == NumberTypes::sacado_rad_dfad
785  >::type>
786  : std::true_type
787  {};
788 
789 
790  template <typename NumberType>
791  struct is_sacado_number<NumberType, typename std::enable_if<
792  is_sacado_dfad_number<NumberType>::value ||
793  is_sacado_rad_number<NumberType>::value
794  >::type>
795  : std::true_type
796  {};
797 
798  }
799 }
800 
801 #endif // DOXYGEN
802 
803 
804 
805 DEAL_II_NAMESPACE_CLOSE
806 
807 
808 #endif // DEAL_II_WITH_TRILINOS
809 
810 #endif
STL namespace.