Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
ad_drivers.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2020 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 
22 #include <deal.II/base/types.h>
23 
26 
28 #include <deal.II/lac/vector.h>
29 
30 #ifdef DEAL_II_WITH_ADOLC
31 
33 # include <adolc/internal/usrparms.h>
35 
36 #endif // DEAL_II_WITH_ADOLC
37 
38 #include <vector>
39 
40 
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 
147  template <typename ADNumberType, typename ScalarType, typename T = void>
149  {
150  // This dummy class definition safely supports compilation
151  // against tapeless numbers or taped number types that have
152  // not yet been implemented.
153 
158 
163  bool
164  is_recording() const;
165 
171  active_tape_index() const;
172 
177  bool
178  is_registered_tape(
179  const typename Types<ADNumberType>::tape_index tape_index) const;
180 
185  bool
186  keep_independent_values() const;
187 
213  void
214  set_tape_buffer_sizes(
215  const typename Types<ADNumberType>::tape_buffer_sizes obufsize = 64 *
216  1024 *
217  1024,
218  const typename Types<ADNumberType>::tape_buffer_sizes lbufsize = 64 *
219  1024 *
220  1024,
221  const typename Types<ADNumberType>::tape_buffer_sizes vbufsize = 64 *
222  1024 *
223  1024,
224  const typename Types<ADNumberType>::tape_buffer_sizes tbufsize = 64 *
225  1024 *
226  1024);
227 
236  void
237  start_taping(const typename Types<ADNumberType>::tape_index tape_index,
238  const bool keep_independent_values);
239 
248  void
249  stop_taping(
250  const typename Types<ADNumberType>::tape_index active_tape_index,
251  const bool write_tapes_to_file);
252 
256  std::vector<typename Types<ADNumberType>::tape_index>
257  get_registered_tape_indices() const;
258 
273  void
274  activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
275 
313  bool
314  requires_retaping(
315  const typename Types<ADNumberType>::tape_index tape_index) const;
316 
350  bool
351  last_action_requires_retaping() const;
352 
356  void
357  remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
358 
365  void
366  reset(const bool clear_registered_tapes);
367 
376  void
377  print(std::ostream &stream) const;
378 
386  void
387  print_tape_stats(
388  const typename Types<ADNumberType>::tape_index tape_index,
389  std::ostream & stream) const;
390 
392 
397 
408  ScalarType
409  value(const typename Types<ADNumberType>::tape_index active_tape_index,
410  const std::vector<ScalarType> &independent_variables) const;
411 
425  void
426  gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
427  const std::vector<ScalarType> &independent_variables,
428  Vector<ScalarType> & gradient) const;
429 
443  void
444  hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
445  const std::vector<ScalarType> &independent_variables,
446  FullMatrix<ScalarType> & hessian) const;
447 
449 
454 
467  void
468  values(const typename Types<ADNumberType>::tape_index active_tape_index,
469  const unsigned int n_dependent_variables,
470  const std::vector<ScalarType> &independent_variables,
471  Vector<ScalarType> & values) const;
472 
491  void
492  jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
493  const unsigned int n_dependent_variables,
494  const std::vector<ScalarType> &independent_variables,
495  FullMatrix<ScalarType> & jacobian) const;
496 
498  };
499 
500 
501 
517  template <typename ADNumberType, typename ScalarType, typename T = void>
519  {
520  // This dummy class definition safely supports compilation
521  // against taped numbers or tapeless number types that have
522  // not yet been implemented.
523 
528 
546  static void
547  initialize_global_environment(const unsigned int n_independent_variables);
548 
550 
555 
560  void
561  allow_dependent_variable_marking();
562 
567  void
568  prevent_dependent_variable_marking();
569 
574  bool
575  is_dependent_variable_marking_allowed() const;
576 
578 
583 
592  ScalarType
593  value(const std::vector<ADNumberType> &dependent_variables) const;
594 
608  void
609  gradient(const std::vector<ADNumberType> &independent_variables,
610  const std::vector<ADNumberType> &dependent_variables,
611  Vector<ScalarType> & gradient) const;
612 
626  void
627  hessian(const std::vector<ADNumberType> &independent_variables,
628  const std::vector<ADNumberType> &dependent_variables,
629  FullMatrix<ScalarType> & hessian) const;
630 
632 
637 
647  void
648  values(const std::vector<ADNumberType> &dependent_variables,
649  Vector<ScalarType> & values) const;
650 
668  void
669  jacobian(const std::vector<ADNumberType> &independent_variables,
670  const std::vector<ADNumberType> &dependent_variables,
671  FullMatrix<ScalarType> & jacobian) const;
672 
674  };
675 
676  } // namespace AD
677 } // namespace Differentiation
678 
679 
680 
681 /* --------------------- inline and template functions --------------------- */
682 
683 
684 #ifndef DOXYGEN
685 
686 namespace Differentiation
687 {
688  namespace AD
689  {
690 # ifdef DEAL_II_WITH_ADOLC
691 
695  template <typename ADNumberType>
696  struct Types<
697  ADNumberType,
698  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
699  NumberTypes::adolc_taped>::type>
700  {
705  using tape_index = unsigned short;
706 
710  using tape_buffer_sizes = unsigned int;
711  }; // struct Types
712 
713 
717  template <typename ADNumberType>
718  struct Numbers<
719  ADNumberType,
720  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
721  NumberTypes::adolc_taped>::type>
722  {
732  static const typename Types<ADNumberType>::tape_index invalid_tape_index =
733  0;
734 
738  // Note: This value is a limitation of ADOL-C, and not something that we
739  // have control over. See test adolc/helper_tape_index_01.cc for
740  // verification that we cannot use or exceed this value. This value is
741  // defined as TBUFNUM; see
742  // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/include/adolc/internal/usrparms.h#L34
743 # ifdef __clang__
744  static const typename Types<ADNumberType>::tape_index max_tape_index =
745  TBUFNUM;
746 # else
747  // For some reason, the test adolc/helper_tape_index_01 indicates that
748  // ADOL-C does not reliably perform correct computations for the full
749  // range of tape indices when GCC is the compiler. So we limit this number
750  // according to the results of the test.
751  static const typename Types<ADNumberType>::tape_index max_tape_index =
752  TBUFNUM - 2;
753 # endif
754  }; // struct Numbers
755 
756 
764  template <typename ADNumberType>
765  struct TapedDrivers<
766  ADNumberType,
767  double,
768  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
769  NumberTypes::adolc_taped>::type>
770  {
771  using scalar_type = double;
772 
776  TapedDrivers();
777 
778 
783 
784  bool
785  is_recording() const;
786 
788  active_tape_index() const;
789 
790  bool
791  keep_independent_values() const;
792 
793  bool
794  is_registered_tape(
795  const typename Types<ADNumberType>::tape_index tape_index) const;
796 
797  void
798  set_tape_buffer_sizes(
799  const typename Types<ADNumberType>::tape_buffer_sizes in_obufsize,
800  const typename Types<ADNumberType>::tape_buffer_sizes in_lbufsize,
801  const typename Types<ADNumberType>::tape_buffer_sizes in_vbufsize,
802  const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize);
803 
804  void
805  start_taping(const typename Types<ADNumberType>::tape_index tape_index,
806  const bool keep_independent_values);
807 
808  void
809  stop_taping(
810  const typename Types<ADNumberType>::tape_index active_tape_index,
811  const bool write_tapes_to_file);
812 
813  std::vector<typename Types<ADNumberType>::tape_index>
814  get_registered_tape_indices() const;
815 
816  void
817  activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
818 
819  bool
820  requires_retaping(
821  const typename Types<ADNumberType>::tape_index tape_index) const;
822 
823  bool
824  last_action_requires_retaping() const;
825 
826  void
827  remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
828 
829  void
830  reset(const bool clear_registered_tapes);
831 
832  void
833  print(std::ostream &stream) const;
834 
835  void
836  print_tape_stats(
837  const typename Types<ADNumberType>::tape_index tape_index,
838  std::ostream & stream) const;
839 
841 
846 
847  scalar_type
848  value(const typename Types<ADNumberType>::tape_index active_tape_index,
849  const std::vector<scalar_type> &independent_variables) const;
850 
851  void
852  gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
853  const std::vector<scalar_type> &independent_variables,
854  Vector<scalar_type> & gradient) const;
855 
856  void
857  hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
858  const std::vector<scalar_type> &independent_variables,
859  FullMatrix<scalar_type> & hessian) const;
860 
862 
867 
868  void
869  values(const typename Types<ADNumberType>::tape_index active_tape_index,
870  const unsigned int n_dependent_variables,
871  const std::vector<scalar_type> &independent_variables,
872  Vector<scalar_type> & values) const;
873 
874  void
875  jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
876  const unsigned int n_dependent_variables,
877  const std::vector<scalar_type> &independent_variables,
878  FullMatrix<scalar_type> & jacobian) const;
879 
881 
882  protected:
888  typename Types<ADNumberType>::tape_index active_tape;
889 
895  bool keep_values;
896 
901  bool is_recording_flag;
902 
930  mutable std::map<typename Types<ADNumberType>::tape_index, int> status;
931 
938  bool use_stored_taped_buffer_sizes;
939 
943  typename Types<ADNumberType>::tape_buffer_sizes obufsize;
944 
948  typename Types<ADNumberType>::tape_buffer_sizes lbufsize;
949 
953  typename Types<ADNumberType>::tape_buffer_sizes vbufsize;
954 
958  typename Types<ADNumberType>::tape_buffer_sizes tbufsize;
959  };
960 
961 # else
962 
971  template <typename ADNumberType>
972  struct TapedDrivers<
973  ADNumberType,
974  double,
975  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
976  NumberTypes::adolc_taped>::type>
977  {
978  using scalar_type = double;
979 
984 
985  bool
986  is_recording() const;
987 
989  active_tape_index() const;
990 
991  bool
992  keep_independent_values() const;
993 
994  bool
995  is_registered_tape(
996  const typename Types<ADNumberType>::tape_index tape_index) const;
997 
998  void
999  set_tape_buffer_sizes(
1004 
1005  void
1006  start_taping(const typename Types<ADNumberType>::tape_index, const bool);
1007 
1008  void
1009  stop_taping(const typename Types<ADNumberType>::tape_index, const bool);
1010 
1011  std::vector<typename Types<ADNumberType>::tape_index>
1012  get_registered_tape_indices() const;
1013 
1014  void
1015  activate_tape(const typename Types<ADNumberType>::tape_index);
1016 
1017  bool
1018  requires_retaping(const typename Types<ADNumberType>::tape_index) const;
1019 
1020  bool
1021  last_action_requires_retaping() const;
1022 
1023  void
1024  remove_tape(const typename Types<ADNumberType>::tape_index);
1025 
1026  void
1027  reset(const bool);
1028 
1029  void
1030  print(std::ostream &stream) const;
1031 
1032  void
1033  print_tape_stats(const typename Types<ADNumberType>::tape_index,
1034  std::ostream &) const;
1035 
1037 
1042 
1043  scalar_type
1044  value(const typename Types<ADNumberType>::tape_index,
1045  const std::vector<scalar_type> &) const;
1046 
1047  void
1048  gradient(const typename Types<ADNumberType>::tape_index,
1049  const std::vector<scalar_type> &,
1050  Vector<scalar_type> &) const;
1051 
1052  void
1053  hessian(const typename Types<ADNumberType>::tape_index,
1054  const std::vector<scalar_type> &,
1055  FullMatrix<scalar_type> &) const;
1056 
1058 
1063 
1064  void
1065  values(const typename Types<ADNumberType>::tape_index,
1066  const unsigned int,
1067  const std::vector<scalar_type> &,
1068  Vector<scalar_type> &) const;
1069 
1070  void
1071  jacobian(const typename Types<ADNumberType>::tape_index,
1072  const unsigned int,
1073  const std::vector<scalar_type> &,
1074  FullMatrix<scalar_type> &) const;
1075 
1077  };
1078 
1079 # endif // DEAL_II_WITH_ADOLC
1080 
1089  template <typename ADNumberType>
1090  struct TapedDrivers<
1091  ADNumberType,
1092  float,
1093  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
1094  NumberTypes::adolc_taped>::type>
1095  {
1096  using scalar_type = float;
1097 
1102 
1103  bool
1104  is_recording() const;
1105 
1107  active_tape_index() const;
1108 
1109  bool
1110  keep_independent_values() const;
1111 
1112  bool
1113  is_registered_tape(
1114  const typename Types<ADNumberType>::tape_index tape_index) const;
1115 
1116  void
1117  set_tape_buffer_sizes(
1118  const typename Types<ADNumberType>::tape_buffer_sizes obufsize,
1119  const typename Types<ADNumberType>::tape_buffer_sizes lbufsize,
1120  const typename Types<ADNumberType>::tape_buffer_sizes vbufsize,
1121  const typename Types<ADNumberType>::tape_buffer_sizes tbufsize);
1122 
1123  void
1124  start_taping(const typename Types<ADNumberType>::tape_index tape_index,
1125  const bool keep_independent_values);
1126 
1127  void
1128  stop_taping(
1129  const typename Types<ADNumberType>::tape_index active_tape_index,
1130  const bool write_tapes_to_file);
1131 
1132  std::vector<typename Types<ADNumberType>::tape_index>
1133  get_registered_tape_indices() const;
1134 
1135  void
1136  activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
1137 
1138  bool
1139  requires_retaping(
1140  const typename Types<ADNumberType>::tape_index tape_index) const;
1141 
1142  bool
1143  last_action_requires_retaping() const;
1144 
1145  void
1146  remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
1147 
1148  void
1149  reset(const bool clear_registered_tapes);
1150 
1151  void
1152  print(std::ostream &stream) const;
1153 
1154  void
1155  print_tape_stats(
1156  const typename Types<ADNumberType>::tape_index tape_index,
1157  std::ostream & stream) const;
1158 
1160 
1165 
1166  scalar_type
1167  value(const typename Types<ADNumberType>::tape_index active_tape_index,
1168  const std::vector<scalar_type> &independent_variables) const;
1169 
1170  void
1171  gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
1172  const std::vector<scalar_type> &independent_variables,
1173  Vector<scalar_type> & gradient) const;
1174 
1175  void
1176  hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
1177  const std::vector<scalar_type> &independent_variables,
1178  FullMatrix<scalar_type> & hessian) const;
1179 
1181 
1186 
1187  void
1188  values(const typename Types<ADNumberType>::tape_index active_tape_index,
1189  const unsigned int n_dependent_variables,
1190  const std::vector<scalar_type> &independent_variables,
1191  Vector<scalar_type> & values) const;
1192 
1193  void
1194  jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
1195  const unsigned int n_dependent_variables,
1196  const std::vector<scalar_type> &independent_variables,
1197  FullMatrix<scalar_type> & jacobian) const;
1198 
1200 
1201  private:
1205  std::vector<double>
1206  vector_float_to_double(const std::vector<float> &in) const;
1207 
1212  };
1213 
1214 
1215  // ------------- TapelessDrivers -------------
1216 
1217 
1223  template <typename ADNumberType, typename ScalarType>
1224  struct TapelessDrivers<
1225  ADNumberType,
1226  ScalarType,
1227  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
1228  NumberTypes::sacado_rad ||
1229  ADNumberTraits<ADNumberType>::type_code ==
1230  NumberTypes::sacado_rad_dfad>::type>
1231  {
1235  TapelessDrivers();
1236 
1241 
1242  static void
1243  initialize_global_environment(const unsigned int n_independent_variables);
1244 
1246 
1251 
1252  void
1253  allow_dependent_variable_marking();
1254 
1255  void
1256  prevent_dependent_variable_marking();
1257 
1258  bool
1259  is_dependent_variable_marking_allowed() const;
1260 
1262 
1267 
1268  ScalarType
1269  value(const std::vector<ADNumberType> &dependent_variables) const;
1270 
1271  void
1272  gradient(const std::vector<ADNumberType> &independent_variables,
1273  const std::vector<ADNumberType> &dependent_variables,
1274  Vector<ScalarType> & gradient) const;
1275 
1276  void
1277  hessian(const std::vector<ADNumberType> &independent_variables,
1278  const std::vector<ADNumberType> &dependent_variables,
1279  FullMatrix<ScalarType> & hessian) const;
1280 
1282 
1287 
1288  void
1289  values(const std::vector<ADNumberType> &dependent_variables,
1290  Vector<ScalarType> & values) const;
1291 
1292  void
1293  jacobian(const std::vector<ADNumberType> &independent_variables,
1294  const std::vector<ADNumberType> &dependent_variables,
1295  FullMatrix<ScalarType> & jacobian) const;
1296 
1298 
1299  private:
1304  bool dependent_variable_marking_safe;
1305  };
1306 
1307 
1313  template <typename ADNumberType, typename ScalarType>
1314  struct TapelessDrivers<
1315  ADNumberType,
1316  ScalarType,
1317  typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
1318  NumberTypes::adolc_tapeless ||
1319  ADNumberTraits<ADNumberType>::type_code ==
1320  NumberTypes::sacado_dfad ||
1321  ADNumberTraits<ADNumberType>::type_code ==
1322  NumberTypes::sacado_dfad_dfad>::type>
1323  {
1327  TapelessDrivers();
1328 
1333 
1334  static void
1335  initialize_global_environment(const unsigned int n_independent_variables);
1336 
1338 
1343 
1344  void
1345  allow_dependent_variable_marking();
1346 
1347  void
1348  prevent_dependent_variable_marking();
1349 
1350  bool
1351  is_dependent_variable_marking_allowed() const;
1352 
1354 
1359 
1360  ScalarType
1361  value(const std::vector<ADNumberType> &dependent_variables) const;
1362 
1363  void
1364  gradient(const std::vector<ADNumberType> &independent_variables,
1365  const std::vector<ADNumberType> &dependent_variables,
1366  Vector<ScalarType> & gradient) const;
1367 
1368  void
1369  hessian(const std::vector<ADNumberType> &independent_variables,
1370  const std::vector<ADNumberType> &dependent_variables,
1371  FullMatrix<ScalarType> & hessian) const;
1372 
1374 
1379 
1380  void
1381  values(const std::vector<ADNumberType> &dependent_variables,
1382  Vector<ScalarType> & values) const;
1383 
1384  void
1385  jacobian(const std::vector<ADNumberType> &independent_variables,
1386  const std::vector<ADNumberType> &dependent_variables,
1387  FullMatrix<ScalarType> & jacobian) const;
1388 
1390 
1391  private:
1396  bool dependent_variable_marking_safe;
1397  };
1398 
1399  } // namespace AD
1400 } // namespace Differentiation
1401 
1402 
1403 #endif // DOXYGEN
1404 
1406 
1407 
1408 #endif
static ::ExceptionBase & ExcRequiresADNumberSpecialization()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
static ::ExceptionBase & ExcSupportedDerivativeLevels(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcRequiresADOLC()
unsigned int tape_buffer_sizes
Definition: ad_drivers.h:104
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:408
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:445
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394