Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10:00+00:00
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
ad_drivers.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_differentiation_ad_ad_drivers_h
16#define dealii_differentiation_ad_ad_drivers_h
17
18#include <deal.II/base/config.h>
19
21#include <deal.II/base/types.h>
22
25
27#include <deal.II/lac/vector.h>
28
29#ifdef DEAL_II_WITH_ADOLC
30
31# include <adolc/adouble.h> // Taped double
32# include <adolc/adtl.h> // Tapeless double
33# include <adolc/internal/usrparms.h>
34
35#endif // DEAL_II_WITH_ADOLC
36
37#include <limits>
38#include <vector>
39
40
42
43
44namespace Differentiation
45{
46 namespace AD
47 {
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
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).
122 0u;
123
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
163 bool
164 is_recording() const;
165
171 active_tape_index() const;
172
177 bool
179 const typename Types<ADNumberType>::tape_index tape_index) const;
180
185 bool
187
213 void
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
251 const bool write_tapes_to_file);
252
256 std::vector<typename Types<ADNumberType>::tape_index>
258
273 void
274 activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
275
313 bool
315 const typename Types<ADNumberType>::tape_index tape_index) const;
316
350 bool
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
388 const typename Types<ADNumberType>::tape_index tape_index,
389 std::ostream &stream) const;
390
408 ScalarType
410 const std::vector<ScalarType> &independent_variables) const;
411
425 void
427 const std::vector<ScalarType> &independent_variables,
429
443 void
445 const std::vector<ScalarType> &independent_variables,
447
467 void
469 const unsigned int n_dependent_variables,
470 const std::vector<ScalarType> &independent_variables,
472
491 void
493 const unsigned int n_dependent_variables,
494 const std::vector<ScalarType> &independent_variables,
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
546 static void
547 initialize_global_environment(const unsigned int n_independent_variables);
548
560 void
562
567 void
569
574 bool
576
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,
612
626 void
627 hessian(const std::vector<ADNumberType> &independent_variables,
628 const std::vector<ADNumberType> &dependent_variables,
630
647 void
648 values(const std::vector<ADNumberType> &dependent_variables,
650
668 void
669 jacobian(const std::vector<ADNumberType> &independent_variables,
670 const std::vector<ADNumberType> &dependent_variables,
672
674 };
675
676 } // namespace AD
677} // namespace Differentiation
678
679
680
681/* --------------------- inline and template functions --------------------- */
682
683
684#ifndef DOXYGEN
685
686namespace Differentiation
687{
688 namespace AD
689 {
690# ifdef DEAL_II_WITH_ADOLC
691
695 template <typename ADNumberType>
696 struct Types<ADNumberType,
697 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
698 NumberTypes::adolc_taped>>
699 {
704 using tape_index = unsigned short;
705
709 using tape_buffer_sizes = unsigned int;
710 }; // struct Types
711
712
716 template <typename ADNumberType>
717 struct Numbers<ADNumberType,
718 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
719 NumberTypes::adolc_taped>>
720 {
731 0;
732
736 // Note: This value is a limitation of ADOL-C, and not something that we
737 // have control over. See test adolc/helper_tape_index_01.cc for
738 // verification that we cannot use or exceed this value. This value is
739 // defined as TBUFNUM; see
740 // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/include/adolc/internal/usrparms.h#L34
741# ifdef __clang__
742 static const typename Types<ADNumberType>::tape_index max_tape_index =
743 TBUFNUM;
744# else
745 // For some reason, the test adolc/helper_tape_index_01 indicates that
746 // ADOL-C does not reliably perform correct computations for the full
747 // range of tape indices when GCC is the compiler. So we limit this number
748 // according to the results of the test.
749 static const typename Types<ADNumberType>::tape_index max_tape_index =
750 TBUFNUM - 2;
751# endif
752 }; // struct Numbers
753
754
762 template <typename ADNumberType>
763 struct TapedDrivers<
764 ADNumberType,
765 double,
766 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
767 NumberTypes::adolc_taped>>
768 {
769 using scalar_type = double;
770
774 TapedDrivers();
775
776
782 bool
783 is_recording() const;
784
786 active_tape_index() const;
787
788 bool
790
791 bool
793 const typename Types<ADNumberType>::tape_index tape_index) const;
794
795 void
797 const typename Types<ADNumberType>::tape_buffer_sizes in_obufsize,
798 const typename Types<ADNumberType>::tape_buffer_sizes in_lbufsize,
799 const typename Types<ADNumberType>::tape_buffer_sizes in_vbufsize,
800 const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize);
801
802 void
803 start_taping(const typename Types<ADNumberType>::tape_index tape_index,
804 const bool keep_independent_values);
805
806 void
809 const bool write_tapes_to_file);
810
811 std::vector<typename Types<ADNumberType>::tape_index>
813
814 void
815 activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
816
817 bool
819 const typename Types<ADNumberType>::tape_index tape_index) const;
820
821 bool
823
824 void
825 remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
826
827 void
828 reset(const bool clear_registered_tapes);
829
830 void
831 print(std::ostream &stream) const;
832
833 void
835 const typename Types<ADNumberType>::tape_index tape_index,
836 std::ostream &stream) const;
837
845 scalar_type
847 const std::vector<scalar_type> &independent_variables) const;
848
849 void
851 const std::vector<scalar_type> &independent_variables,
853
854 void
856 const std::vector<scalar_type> &independent_variables,
858
866 void
868 const unsigned int n_dependent_variables,
869 const std::vector<scalar_type> &independent_variables,
870 Vector<scalar_type> &values) const;
871
872 void
874 const unsigned int n_dependent_variables,
875 const std::vector<scalar_type> &independent_variables,
877
880 protected:
886 typename Types<ADNumberType>::tape_index active_tape;
887
893 bool keep_values;
894
899 bool is_recording_flag;
900
928 mutable std::map<typename Types<ADNumberType>::tape_index, int> status;
929
936 bool use_stored_taped_buffer_sizes;
937
942
947
952
957 };
958
959# else
960
969 template <typename ADNumberType>
970 struct TapedDrivers<
971 ADNumberType,
972 double,
973 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
974 NumberTypes::adolc_taped>>
975 {
976 using scalar_type = double;
977
983 bool
984 is_recording() const;
985
987 active_tape_index() const;
988
989 bool
991
992 bool
994 const typename Types<ADNumberType>::tape_index tape_index) const;
995
996 void
1002
1003 void
1004 start_taping(const typename Types<ADNumberType>::tape_index, const bool);
1005
1006 void
1007 stop_taping(const typename Types<ADNumberType>::tape_index, const bool);
1008
1009 std::vector<typename Types<ADNumberType>::tape_index>
1011
1012 void
1014
1015 bool
1017
1018 bool
1020
1021 void
1023
1024 void
1025 reset(const bool);
1026
1027 void
1028 print(std::ostream &stream) const;
1029
1030 void
1032 std::ostream &) const;
1033
1041 scalar_type
1043 const std::vector<scalar_type> &) const;
1044
1045 void
1047 const std::vector<scalar_type> &,
1048 Vector<scalar_type> &) const;
1049
1050 void
1052 const std::vector<scalar_type> &,
1053 FullMatrix<scalar_type> &) const;
1054
1062 void
1064 const unsigned int,
1065 const std::vector<scalar_type> &,
1066 Vector<scalar_type> &) const;
1067
1068 void
1070 const unsigned int,
1071 const std::vector<scalar_type> &,
1072 FullMatrix<scalar_type> &) const;
1073
1075 };
1076
1077# endif // DEAL_II_WITH_ADOLC
1078
1087 template <typename ADNumberType>
1088 struct TapedDrivers<
1089 ADNumberType,
1090 float,
1091 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1092 NumberTypes::adolc_taped>>
1093 {
1094 using scalar_type = float;
1095
1101 bool
1102 is_recording() const;
1103
1105 active_tape_index() const;
1106
1107 bool
1109
1110 bool
1112 const typename Types<ADNumberType>::tape_index tape_index) const;
1113
1114 void
1116 const typename Types<ADNumberType>::tape_buffer_sizes obufsize,
1117 const typename Types<ADNumberType>::tape_buffer_sizes lbufsize,
1118 const typename Types<ADNumberType>::tape_buffer_sizes vbufsize,
1119 const typename Types<ADNumberType>::tape_buffer_sizes tbufsize);
1120
1121 void
1122 start_taping(const typename Types<ADNumberType>::tape_index tape_index,
1123 const bool keep_independent_values);
1124
1125 void
1128 const bool write_tapes_to_file);
1129
1130 std::vector<typename Types<ADNumberType>::tape_index>
1132
1133 void
1134 activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
1135
1136 bool
1138 const typename Types<ADNumberType>::tape_index tape_index) const;
1139
1140 bool
1142
1143 void
1144 remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
1145
1146 void
1147 reset(const bool clear_registered_tapes);
1148
1149 void
1150 print(std::ostream &stream) const;
1151
1152 void
1154 const typename Types<ADNumberType>::tape_index tape_index,
1155 std::ostream &stream) const;
1156
1164 scalar_type
1166 const std::vector<scalar_type> &independent_variables) const;
1167
1168 void
1170 const std::vector<scalar_type> &independent_variables,
1172
1173 void
1175 const std::vector<scalar_type> &independent_variables,
1177
1185 void
1187 const unsigned int n_dependent_variables,
1188 const std::vector<scalar_type> &independent_variables,
1189 Vector<scalar_type> &values) const;
1190
1191 void
1193 const unsigned int n_dependent_variables,
1194 const std::vector<scalar_type> &independent_variables,
1196
1199 private:
1203 std::vector<double>
1204 vector_float_to_double(const std::vector<float> &in) const;
1205
1209 TapedDrivers<ADNumberType, double> taped_driver;
1210 };
1211
1212
1213 // ------------- TapelessDrivers -------------
1214
1215
1221 template <typename ADNumberType, typename ScalarType>
1222 struct TapelessDrivers<
1223 ADNumberType,
1224 ScalarType,
1225 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1226 NumberTypes::sacado_rad ||
1227 ADNumberTraits<ADNumberType>::type_code ==
1228 NumberTypes::sacado_rad_dfad>>
1229 {
1233 TapelessDrivers();
1234
1240 static void
1241 initialize_global_environment(const unsigned int n_independent_variables);
1242
1250 void
1252
1253 void
1255
1256 bool
1258
1266 ScalarType
1267 value(const std::vector<ADNumberType> &dependent_variables) const;
1268
1269 void
1270 gradient(const std::vector<ADNumberType> &independent_variables,
1271 const std::vector<ADNumberType> &dependent_variables,
1273
1274 void
1275 hessian(const std::vector<ADNumberType> &independent_variables,
1276 const std::vector<ADNumberType> &dependent_variables,
1278
1286 void
1287 values(const std::vector<ADNumberType> &dependent_variables,
1288 Vector<ScalarType> &values) const;
1289
1290 void
1291 jacobian(const std::vector<ADNumberType> &independent_variables,
1292 const std::vector<ADNumberType> &dependent_variables,
1294
1297 private:
1302 bool dependent_variable_marking_safe;
1303 };
1304
1305
1311 template <typename ADNumberType, typename ScalarType>
1312 struct TapelessDrivers<
1313 ADNumberType,
1314 ScalarType,
1315 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1316 NumberTypes::adolc_tapeless ||
1317 ADNumberTraits<ADNumberType>::type_code ==
1318 NumberTypes::sacado_dfad ||
1319 ADNumberTraits<ADNumberType>::type_code ==
1320 NumberTypes::sacado_dfad_dfad>>
1321 {
1325 TapelessDrivers();
1326
1332 static void
1333 initialize_global_environment(const unsigned int n_independent_variables);
1334
1342 void
1344
1345 void
1347
1348 bool
1350
1358 ScalarType
1359 value(const std::vector<ADNumberType> &dependent_variables) const;
1360
1361 void
1362 gradient(const std::vector<ADNumberType> &independent_variables,
1363 const std::vector<ADNumberType> &dependent_variables,
1365
1366 void
1367 hessian(const std::vector<ADNumberType> &independent_variables,
1368 const std::vector<ADNumberType> &dependent_variables,
1370
1378 void
1379 values(const std::vector<ADNumberType> &dependent_variables,
1380 Vector<ScalarType> &values) const;
1381
1382 void
1383 jacobian(const std::vector<ADNumberType> &independent_variables,
1384 const std::vector<ADNumberType> &dependent_variables,
1386
1389 private:
1394 bool dependent_variable_marking_safe;
1395 };
1396
1397 } // namespace AD
1398} // namespace Differentiation
1399
1400
1401#endif // DOXYGEN
1402
1404
1405
1406#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcSupportedDerivativeLevels(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcRequiresADNumberSpecialization()
static ::ExceptionBase & ExcRequiresADOLC()
#define DeclException2(Exception2, type1, type2, outsequence)
#define DeclExceptionMsg(Exception, defaulttext)
STL namespace.
static const Types< ADNumberType >::tape_index max_tape_index
Definition ad_drivers.h:127
static const Types< ADNumberType >::tape_index invalid_tape_index
Definition ad_drivers.h:121
void gradient(const typename Types< ADNumberType >::tape_index active_tape_index, const std::vector< ScalarType > &independent_variables, Vector< ScalarType > &gradient) const
bool is_registered_tape(const typename Types< ADNumberType >::tape_index tape_index) const
Definition ad_drivers.cc:71
void stop_taping(const typename Types< ADNumberType >::tape_index active_tape_index, const bool write_tapes_to_file)
void remove_tape(const typename Types< ADNumberType >::tape_index tape_index)
void start_taping(const typename Types< ADNumberType >::tape_index tape_index, const bool keep_independent_values)
void hessian(const typename Types< ADNumberType >::tape_index active_tape_index, const std::vector< ScalarType > &independent_variables, FullMatrix< ScalarType > &hessian) const
void set_tape_buffer_sizes(const typename Types< ADNumberType >::tape_buffer_sizes obufsize=64 *1024 *1024, const typename Types< ADNumberType >::tape_buffer_sizes lbufsize=64 *1024 *1024, const typename Types< ADNumberType >::tape_buffer_sizes vbufsize=64 *1024 *1024, const typename Types< ADNumberType >::tape_buffer_sizes tbufsize=64 *1024 *1024)
Definition ad_drivers.cc:90
void activate_tape(const typename Types< ADNumberType >::tape_index tape_index)
bool requires_retaping(const typename Types< ADNumberType >::tape_index tape_index) const
void print(std::ostream &stream) const
void print_tape_stats(const typename Types< ADNumberType >::tape_index tape_index, std::ostream &stream) const
ScalarType value(const typename Types< ADNumberType >::tape_index active_tape_index, const std::vector< ScalarType > &independent_variables) const
Types< ADNumberType >::tape_index active_tape_index() const
Definition ad_drivers.cc:62
void reset(const bool clear_registered_tapes)
std::vector< typename Types< ADNumberType >::tape_index > get_registered_tape_indices() const
void values(const typename Types< ADNumberType >::tape_index active_tape_index, const unsigned int n_dependent_variables, const std::vector< ScalarType > &independent_variables, Vector< ScalarType > &values) const
void jacobian(const typename Types< ADNumberType >::tape_index active_tape_index, const unsigned int n_dependent_variables, const std::vector< ScalarType > &independent_variables, FullMatrix< ScalarType > &jacobian) const
void gradient(const std::vector< ADNumberType > &independent_variables, const std::vector< ADNumberType > &dependent_variables, Vector< ScalarType > &gradient) const
void hessian(const std::vector< ADNumberType > &independent_variables, const std::vector< ADNumberType > &dependent_variables, FullMatrix< ScalarType > &hessian) const
void values(const std::vector< ADNumberType > &dependent_variables, Vector< ScalarType > &values) const
void jacobian(const std::vector< ADNumberType > &independent_variables, const std::vector< ADNumberType > &dependent_variables, FullMatrix< ScalarType > &jacobian) const
static void initialize_global_environment(const unsigned int n_independent_variables)
ScalarType value(const std::vector< ADNumberType > &dependent_variables) const
unsigned int tape_buffer_sizes
Definition ad_drivers.h:104