Reference documentation for deal.II version Git 54f37b890c 2020-02-28 19:52:28 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
ad_drivers.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2018 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_ad_drivers_h
17 #define dealii_differentiation_ad_ad_drivers_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/types.h>
23 
24 #include <deal.II/differentiation/ad/ad_number_traits.h>
25 #include <deal.II/differentiation/ad/ad_number_types.h>
26 
27 #include <deal.II/lac/full_matrix.h>
28 #include <deal.II/lac/vector.h>
29 
30 #ifdef DEAL_II_WITH_ADOLC
31 
32 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
33 # include <adolc/internal/usrparms.h>
34 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
35 
36 #endif // DEAL_II_WITH_ADOLC
37 
38 #include <vector>
39 
40 
41 DEAL_II_NAMESPACE_OPEN
42 
43 
44 namespace Differentiation
45 {
46  namespace AD
47  {
52 
59  "This function is called in a class that is expected to be specialized "
60  "for auto-differentiable numbers.");
61 
67  "This function is only available if deal.II is compiled with ADOL-C.");
68 
79  std::size_t,
80  std::size_t,
81  << "The number of derivative levels that this auto-differentiable number type supports is "
82  << arg1
83  << ", but to perform the intended operation the number must support at least "
84  << arg2 << " levels.");
85 
87 
88 
93  template <typename ADNumberType, typename T = void>
94  struct Types
95  {
99  using tape_index = unsigned int;
100 
104  using tape_buffer_sizes = unsigned int;
105  }; // struct Types
106 
107 
112  template <typename ADNumberType, typename T = void>
113  struct Numbers
114  {
119  // Note: We stipulate in the documentation for the helper classes that the
120  // valid tape index range is (invalid_tape_index,max_tape_index).
121  static const typename Types<ADNumberType>::tape_index invalid_tape_index =
122  0u;
123 
127  static const typename Types<ADNumberType>::tape_index max_tape_index =
128  std::numeric_limits<typename Types<ADNumberType>::tape_index>::max();
129  }; // struct Numbers
130 
131 
149  template <typename ADNumberType, typename ScalarType, typename T = void>
151  {
152  // This dummy class definition safely supports compilation
153  // against tapeless numbers or taped number types that have
154  // not yet been implemented.
155 
160 
165  bool
166  is_recording() const;
167 
173  active_tape_index() const;
174 
179  bool
180  is_registered_tape(
181  const typename Types<ADNumberType>::tape_index tape_index) const;
182 
187  bool
188  keep_independent_values() const;
189 
215  void
216  set_tape_buffer_sizes(
217  const typename Types<ADNumberType>::tape_buffer_sizes obufsize = 64 *
218  1024 *
219  1024,
220  const typename Types<ADNumberType>::tape_buffer_sizes lbufsize = 64 *
221  1024 *
222  1024,
223  const typename Types<ADNumberType>::tape_buffer_sizes vbufsize = 64 *
224  1024 *
225  1024,
226  const typename Types<ADNumberType>::tape_buffer_sizes tbufsize = 64 *
227  1024 *
228  1024);
229 
238  void
239  start_taping(const typename Types<ADNumberType>::tape_index tape_index,
240  const bool keep_independent_values);
241 
250  void
251  stop_taping(
252  const typename Types<ADNumberType>::tape_index active_tape_index,
253  const bool write_tapes_to_file);
254 
258  std::vector<typename Types<ADNumberType>::tape_index>
259  get_registered_tape_indices() const;
260 
275  void
276  activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
277 
315  bool
316  requires_retaping(
317  const typename Types<ADNumberType>::tape_index tape_index) const;
318 
352  bool
353  last_action_requires_retaping() const;
354 
358  void
359  remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
360 
367  void
368  reset(const bool clear_registered_tapes);
369 
378  void
379  print(std::ostream &stream) const;
380 
388  void
389  print_tape_stats(
390  const typename Types<ADNumberType>::tape_index tape_index,
391  std::ostream & stream) const;
392 
394 
399 
410  ScalarType
411  value(const typename Types<ADNumberType>::tape_index active_tape_index,
412  const std::vector<ScalarType> &independent_variables) const;
413 
427  void
428  gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
429  const std::vector<ScalarType> &independent_variables,
430  Vector<ScalarType> & gradient) const;
431 
445  void
446  hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
447  const std::vector<ScalarType> &independent_variables,
448  FullMatrix<ScalarType> & hessian) const;
449 
451 
456 
469  void
470  values(const typename Types<ADNumberType>::tape_index active_tape_index,
471  const unsigned int n_dependent_variables,
472  const std::vector<ScalarType> &independent_variables,
473  Vector<ScalarType> & values) const;
474 
493  void
494  jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
495  const unsigned int n_dependent_variables,
496  const std::vector<ScalarType> &independent_variables,
497  FullMatrix<ScalarType> & jacobian) const;
498 
500  };
501 
502 
503 
521  template <typename ADNumberType, typename ScalarType, typename T = void>
523  {
524  // This dummy class definition safely supports compilation
525  // against taped numbers or tapeless number types that have
526  // not yet been implemented.
527 
532 
550  static void
551  initialize_global_environment(const unsigned int n_independent_variables);
552 
554 
559 
564  void
565  allow_dependent_variable_marking();
566 
571  void
572  prevent_dependent_variable_marking();
573 
578  bool
579  is_dependent_variable_marking_allowed() const;
580 
582 
587 
596  ScalarType
597  value(const std::vector<ADNumberType> &dependent_variables) const;
598 
612  void
613  gradient(const std::vector<ADNumberType> &independent_variables,
614  const std::vector<ADNumberType> &dependent_variables,
615  Vector<ScalarType> & gradient) const;
616 
630  void
631  hessian(const std::vector<ADNumberType> &independent_variables,
632  const std::vector<ADNumberType> &dependent_variables,
633  FullMatrix<ScalarType> & hessian) const;
634 
636 
641 
651  void
652  values(const std::vector<ADNumberType> &dependent_variables,
653  Vector<ScalarType> & values) const;
654 
672  void
673  jacobian(const std::vector<ADNumberType> &independent_variables,
674  const std::vector<ADNumberType> &dependent_variables,
675  FullMatrix<ScalarType> & jacobian) const;
676 
678  };
679 
680  } // namespace AD
681 } // namespace Differentiation
682 
683 
684 
685 /* --------------------- inline and template functions --------------------- */
686 
687 
688 #ifndef DOXYGEN
689 
690 namespace Differentiation
691 {
692  namespace AD
693  {
694 # ifdef DEAL_II_WITH_ADOLC
695 
699  template <typename ADNumberType>
700  struct Types<
701  ADNumberType,
702  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
703  NumberTypes::adolc_taped>::type>
704  {
709  using tape_index = unsigned short;
710 
714  using tape_buffer_sizes = unsigned int;
715  }; // struct Types
716 
717 
721  template <typename ADNumberType>
722  struct Numbers<
723  ADNumberType,
724  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
725  NumberTypes::adolc_taped>::type>
726  {
736  static const typename Types<ADNumberType>::tape_index invalid_tape_index =
737  0;
738 
742  // Note: This value is a limitation of ADOL-C, and not something that we
743  // have control over. See test adolc/helper_tape_index_01.cc for
744  // verification that we cannot use or exceed this value. This value is
745  // defined as TBUFNUM; see
746  // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/include/adolc/internal/usrparms.h#L34
747 # ifdef __clang__
748  static const typename Types<ADNumberType>::tape_index max_tape_index =
749  TBUFNUM;
750 # else
751  // For some reason, the test adolc/helper_tape_index_01 indicates that
752  // ADOL-C does not reliably perform correct computations for the full
753  // range of tape indices when GCC is the compiler. So we limit this number
754  // according to the results of the test.
755  static const typename Types<ADNumberType>::tape_index max_tape_index =
756  TBUFNUM - 2;
757 # endif
758  }; // struct Numbers
759 
760 
768  template <typename ADNumberType>
769  struct TapedDrivers<
770  ADNumberType,
771  double,
772  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
773  NumberTypes::adolc_taped>::type>
774  {
775  using scalar_type = double;
776 
780  TapedDrivers();
781 
782 
787 
788  bool
789  is_recording() const;
790 
792  active_tape_index() const;
793 
794  bool
795  keep_independent_values() const;
796 
797  bool
798  is_registered_tape(
799  const typename Types<ADNumberType>::tape_index tape_index) const;
800 
801  void
802  set_tape_buffer_sizes(
803  const typename Types<ADNumberType>::tape_buffer_sizes in_obufsize,
804  const typename Types<ADNumberType>::tape_buffer_sizes in_lbufsize,
805  const typename Types<ADNumberType>::tape_buffer_sizes in_vbufsize,
806  const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize);
807 
808  void
809  start_taping(const typename Types<ADNumberType>::tape_index tape_index,
810  const bool keep_independent_values);
811 
812  void
813  stop_taping(
814  const typename Types<ADNumberType>::tape_index active_tape_index,
815  const bool write_tapes_to_file);
816 
817  std::vector<typename Types<ADNumberType>::tape_index>
818  get_registered_tape_indices() const;
819 
820  void
821  activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
822 
823  bool
824  requires_retaping(
825  const typename Types<ADNumberType>::tape_index tape_index) const;
826 
827  bool
828  last_action_requires_retaping() const;
829 
830  void
831  remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
832 
833  void
834  reset(const bool clear_registered_tapes);
835 
836  void
837  print(std::ostream &stream) const;
838 
839  void
840  print_tape_stats(
841  const typename Types<ADNumberType>::tape_index tape_index,
842  std::ostream & stream) const;
843 
845 
850 
851  scalar_type
852  value(const typename Types<ADNumberType>::tape_index active_tape_index,
853  const std::vector<scalar_type> &independent_variables) const;
854 
855  void
856  gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
857  const std::vector<scalar_type> &independent_variables,
858  Vector<scalar_type> & gradient) const;
859 
860  void
861  hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
862  const std::vector<scalar_type> &independent_variables,
863  FullMatrix<scalar_type> & hessian) const;
864 
866 
871 
872  void
873  values(const typename Types<ADNumberType>::tape_index active_tape_index,
874  const unsigned int n_dependent_variables,
875  const std::vector<scalar_type> &independent_variables,
876  Vector<scalar_type> & values) const;
877 
878  void
879  jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
880  const unsigned int n_dependent_variables,
881  const std::vector<scalar_type> &independent_variables,
882  FullMatrix<scalar_type> & jacobian) const;
883 
885 
886  protected:
892  typename Types<ADNumberType>::tape_index active_tape;
893 
899  bool keep_values;
900 
905  bool is_recording_flag;
906 
934  mutable std::map<typename Types<ADNumberType>::tape_index, int> status;
935 
942  bool use_stored_taped_buffer_sizes;
943 
947  typename Types<ADNumberType>::tape_buffer_sizes obufsize;
948 
952  typename Types<ADNumberType>::tape_buffer_sizes lbufsize;
953 
957  typename Types<ADNumberType>::tape_buffer_sizes vbufsize;
958 
962  typename Types<ADNumberType>::tape_buffer_sizes tbufsize;
963  };
964 
965 # else
966 
975  template <typename ADNumberType>
976  struct TapedDrivers<
977  ADNumberType,
978  double,
979  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
980  NumberTypes::adolc_taped>::type>
981  {
982  using scalar_type = double;
983 
988 
989  bool
990  is_recording() const;
991 
993  active_tape_index() const;
994 
995  bool
996  keep_independent_values() const;
997 
998  bool
999  is_registered_tape(
1000  const typename Types<ADNumberType>::tape_index tape_index) const;
1001 
1002  void
1003  set_tape_buffer_sizes(
1008 
1009  void
1010  start_taping(const typename Types<ADNumberType>::tape_index, const bool);
1011 
1012  void
1013  stop_taping(const typename Types<ADNumberType>::tape_index, const bool);
1014 
1015  std::vector<typename Types<ADNumberType>::tape_index>
1016  get_registered_tape_indices() const;
1017 
1018  void
1019  activate_tape(const typename Types<ADNumberType>::tape_index);
1020 
1021  bool
1022  requires_retaping(const typename Types<ADNumberType>::tape_index) const;
1023 
1024  bool
1025  last_action_requires_retaping() const;
1026 
1027  void
1028  remove_tape(const typename Types<ADNumberType>::tape_index);
1029 
1030  void
1031  reset(const bool);
1032 
1033  void
1034  print(std::ostream &stream) const;
1035 
1036  void
1037  print_tape_stats(const typename Types<ADNumberType>::tape_index,
1038  std::ostream &) const;
1039 
1041 
1046 
1047  scalar_type
1048  value(const typename Types<ADNumberType>::tape_index,
1049  const std::vector<scalar_type> &) const;
1050 
1051  void
1052  gradient(const typename Types<ADNumberType>::tape_index,
1053  const std::vector<scalar_type> &,
1054  Vector<scalar_type> &) const;
1055 
1056  void
1057  hessian(const typename Types<ADNumberType>::tape_index,
1058  const std::vector<scalar_type> &,
1059  FullMatrix<scalar_type> &) const;
1060 
1062 
1067 
1068  void
1069  values(const typename Types<ADNumberType>::tape_index,
1070  const unsigned int,
1071  const std::vector<scalar_type> &,
1072  Vector<scalar_type> &) const;
1073 
1074  void
1075  jacobian(const typename Types<ADNumberType>::tape_index,
1076  const unsigned int,
1077  const std::vector<scalar_type> &,
1078  FullMatrix<scalar_type> &) const;
1079 
1081  };
1082 
1083 # endif // DEAL_II_WITH_ADOLC
1084 
1093  template <typename ADNumberType>
1094  struct TapedDrivers<
1095  ADNumberType,
1096  float,
1097  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
1098  NumberTypes::adolc_taped>::type>
1099  {
1100  using scalar_type = float;
1101 
1106 
1107  bool
1108  is_recording() const;
1109 
1111  active_tape_index() const;
1112 
1113  bool
1114  keep_independent_values() const;
1115 
1116  bool
1117  is_registered_tape(
1118  const typename Types<ADNumberType>::tape_index tape_index) const;
1119 
1120  void
1121  set_tape_buffer_sizes(
1122  const typename Types<ADNumberType>::tape_buffer_sizes obufsize,
1123  const typename Types<ADNumberType>::tape_buffer_sizes lbufsize,
1124  const typename Types<ADNumberType>::tape_buffer_sizes vbufsize,
1125  const typename Types<ADNumberType>::tape_buffer_sizes tbufsize);
1126 
1127  void
1128  start_taping(const typename Types<ADNumberType>::tape_index tape_index,
1129  const bool keep_independent_values);
1130 
1131  void
1132  stop_taping(
1133  const typename Types<ADNumberType>::tape_index active_tape_index,
1134  const bool write_tapes_to_file);
1135 
1136  std::vector<typename Types<ADNumberType>::tape_index>
1137  get_registered_tape_indices() const;
1138 
1139  void
1140  activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
1141 
1142  bool
1143  requires_retaping(
1144  const typename Types<ADNumberType>::tape_index tape_index) const;
1145 
1146  bool
1147  last_action_requires_retaping() const;
1148 
1149  void
1150  remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
1151 
1152  void
1153  reset(const bool clear_registered_tapes);
1154 
1155  void
1156  print(std::ostream &stream) const;
1157 
1158  void
1159  print_tape_stats(
1160  const typename Types<ADNumberType>::tape_index tape_index,
1161  std::ostream & stream) const;
1162 
1164 
1169 
1170  scalar_type
1171  value(const typename Types<ADNumberType>::tape_index active_tape_index,
1172  const std::vector<scalar_type> &independent_variables) const;
1173 
1174  void
1175  gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
1176  const std::vector<scalar_type> &independent_variables,
1177  Vector<scalar_type> & gradient) const;
1178 
1179  void
1180  hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
1181  const std::vector<scalar_type> &independent_variables,
1182  FullMatrix<scalar_type> & hessian) const;
1183 
1185 
1190 
1191  void
1192  values(const typename Types<ADNumberType>::tape_index active_tape_index,
1193  const unsigned int n_dependent_variables,
1194  const std::vector<scalar_type> &independent_variables,
1195  Vector<scalar_type> & values) const;
1196 
1197  void
1198  jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
1199  const unsigned int n_dependent_variables,
1200  const std::vector<scalar_type> &independent_variables,
1201  FullMatrix<scalar_type> & jacobian) const;
1202 
1204 
1205  private:
1209  std::vector<double>
1210  vector_float_to_double(const std::vector<float> &in) const;
1211 
1216  };
1217 
1218 
1219  // ------------- TapelessDrivers -------------
1220 
1221 
1227  template <typename ADNumberType, typename ScalarType>
1228  struct TapelessDrivers<
1229  ADNumberType,
1230  ScalarType,
1231  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
1232  NumberTypes::sacado_rad ||
1233  ADNumberTraits<ADNumberType>::type_code ==
1234  NumberTypes::sacado_rad_dfad>::type>
1235  {
1239  TapelessDrivers();
1240 
1245 
1246  static void
1247  initialize_global_environment(const unsigned int n_independent_variables);
1248 
1250 
1255 
1256  void
1257  allow_dependent_variable_marking();
1258 
1259  void
1260  prevent_dependent_variable_marking();
1261 
1262  bool
1263  is_dependent_variable_marking_allowed() const;
1264 
1266 
1271 
1272  ScalarType
1273  value(const std::vector<ADNumberType> &dependent_variables) const;
1274 
1275  void
1276  gradient(const std::vector<ADNumberType> &independent_variables,
1277  const std::vector<ADNumberType> &dependent_variables,
1278  Vector<ScalarType> & gradient) const;
1279 
1280  void
1281  hessian(const std::vector<ADNumberType> &independent_variables,
1282  const std::vector<ADNumberType> &dependent_variables,
1283  FullMatrix<ScalarType> & hessian) const;
1284 
1286 
1291 
1292  void
1293  values(const std::vector<ADNumberType> &dependent_variables,
1294  Vector<ScalarType> & values) const;
1295 
1296  void
1297  jacobian(const std::vector<ADNumberType> &independent_variables,
1298  const std::vector<ADNumberType> &dependent_variables,
1299  FullMatrix<ScalarType> & jacobian) const;
1300 
1302 
1303  private:
1308  bool dependent_variable_marking_safe;
1309  };
1310 
1311 
1317  template <typename ADNumberType, typename ScalarType>
1318  struct TapelessDrivers<
1319  ADNumberType,
1320  ScalarType,
1321  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
1322  NumberTypes::adolc_tapeless ||
1323  ADNumberTraits<ADNumberType>::type_code ==
1324  NumberTypes::sacado_dfad ||
1325  ADNumberTraits<ADNumberType>::type_code ==
1326  NumberTypes::sacado_dfad_dfad>::type>
1327  {
1331  TapelessDrivers();
1332 
1337 
1338  static void
1339  initialize_global_environment(const unsigned int n_independent_variables);
1340 
1342 
1347 
1348  void
1349  allow_dependent_variable_marking();
1350 
1351  void
1352  prevent_dependent_variable_marking();
1353 
1354  bool
1355  is_dependent_variable_marking_allowed() const;
1356 
1358 
1363 
1364  ScalarType
1365  value(const std::vector<ADNumberType> &dependent_variables) const;
1366 
1367  void
1368  gradient(const std::vector<ADNumberType> &independent_variables,
1369  const std::vector<ADNumberType> &dependent_variables,
1370  Vector<ScalarType> & gradient) const;
1371 
1372  void
1373  hessian(const std::vector<ADNumberType> &independent_variables,
1374  const std::vector<ADNumberType> &dependent_variables,
1375  FullMatrix<ScalarType> & hessian) const;
1376 
1378 
1383 
1384  void
1385  values(const std::vector<ADNumberType> &dependent_variables,
1386  Vector<ScalarType> & values) const;
1387 
1388  void
1389  jacobian(const std::vector<ADNumberType> &independent_variables,
1390  const std::vector<ADNumberType> &dependent_variables,
1391  FullMatrix<ScalarType> & jacobian) const;
1392 
1394 
1395  private:
1400  bool dependent_variable_marking_safe;
1401  };
1402 
1403  } // namespace AD
1404 } // namespace Differentiation
1405 
1406 
1407 #endif // DOXYGEN
1408 
1409 DEAL_II_NAMESPACE_CLOSE
1410 
1411 
1412 #endif
static ::ExceptionBase & ExcRequiresADNumberSpecialization()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
static ::ExceptionBase & ExcSupportedDerivativeLevels(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcRequiresADOLC()
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496