Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
ad_drivers.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2018 - 2023 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#include <deal.II/base/config.h>
17
18#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
19
20# include <deal.II/base/exceptions.h>
21# include <deal.II/base/types.h>
22# include <deal.II/base/utilities.h>
23
24# include <deal.II/differentiation/ad/ad_drivers.h>
25# include <deal.II/differentiation/ad/ad_number_traits.h>
26# include <deal.II/differentiation/ad/ad_number_types.h>
27# include <deal.II/differentiation/ad/adolc_number_types.h>
28# include <deal.II/differentiation/ad/sacado_number_types.h>
29
30# include <deal.II/lac/full_matrix.h>
31# include <deal.II/lac/vector.h>
32
33# ifdef DEAL_II_WITH_ADOLC
34# include <adolc/adolc_fatalerror.h>
35# include <adolc/drivers/drivers.h>
36# include <adolc/taping.h>
37# endif // DEAL_II_WITH_ADOLC
38
39# include <vector>
40
41
42DEAL_II_NAMESPACE_OPEN
43
44
45namespace Differentiation
46{
47 namespace AD
48 {
49 // ------------- TapedDrivers -------------
50
51
52 template <typename ADNumberType, typename ScalarType, typename T>
53 bool
54 TapedDrivers<ADNumberType, ScalarType, T>::is_recording() const
55 {
56 AssertThrow(false, ExcRequiresADNumberSpecialization());
57 return false;
58 }
59
60
61 template <typename ADNumberType, typename ScalarType, typename T>
62 typename Types<ADNumberType>::tape_index
63 TapedDrivers<ADNumberType, ScalarType, T>::active_tape_index() const
64 {
65 AssertThrow(false, ExcRequiresADNumberSpecialization());
66 return Numbers<ADNumberType>::invalid_tape_index;
67 }
68
69
70 template <typename ADNumberType, typename ScalarType, typename T>
71 bool
72 TapedDrivers<ADNumberType, ScalarType, T>::is_registered_tape(
73 const typename Types<ADNumberType>::tape_index) const
74 {
75 AssertThrow(false, ExcRequiresADNumberSpecialization());
76 return false;
77 }
78
79
80 template <typename ADNumberType, typename ScalarType, typename T>
81 bool
82 TapedDrivers<ADNumberType, ScalarType, T>::keep_independent_values() const
83 {
84 AssertThrow(false, ExcRequiresADNumberSpecialization());
85 return false;
86 }
87
88
89 template <typename ADNumberType, typename ScalarType, typename T>
90 void
91 TapedDrivers<ADNumberType, ScalarType, T>::set_tape_buffer_sizes(
92 const typename Types<ADNumberType>::tape_buffer_sizes,
93 const typename Types<ADNumberType>::tape_buffer_sizes,
94 const typename Types<ADNumberType>::tape_buffer_sizes,
95 const typename Types<ADNumberType>::tape_buffer_sizes)
96 {
97 AssertThrow(false, ExcRequiresADNumberSpecialization());
98 }
99
100
101 template <typename ADNumberType, typename ScalarType, typename T>
102 void
103 TapedDrivers<ADNumberType, ScalarType, T>::start_taping(
104 const typename Types<ADNumberType>::tape_index,
105 const bool)
106 {
107 AssertThrow(false, ExcRequiresADNumberSpecialization());
108 }
109
110
111 template <typename ADNumberType, typename ScalarType, typename T>
112 void
113 TapedDrivers<ADNumberType, ScalarType, T>::stop_taping(
114 const typename Types<ADNumberType>::tape_index,
115 const bool)
116 {
117 AssertThrow(false, ExcRequiresADNumberSpecialization());
118 }
119
120
121 template <typename ADNumberType, typename ScalarType, typename T>
122 std::vector<typename Types<ADNumberType>::tape_index>
123 TapedDrivers<ADNumberType, ScalarType, T>::get_registered_tape_indices()
124 const
125 {
126 AssertThrow(false, ExcRequiresADNumberSpecialization());
127 return std::vector<typename Types<ADNumberType>::tape_index>();
128 }
129
130
131 template <typename ADNumberType, typename ScalarType, typename T>
132 void
133 TapedDrivers<ADNumberType, ScalarType, T>::activate_tape(
134 const typename Types<ADNumberType>::tape_index)
135 {
136 AssertThrow(false, ExcRequiresADNumberSpecialization());
137 }
138
139
140 template <typename ADNumberType, typename ScalarType, typename T>
141 bool
142 TapedDrivers<ADNumberType, ScalarType, T>::requires_retaping(
143 const typename Types<ADNumberType>::tape_index) const
144 {
145 AssertThrow(false, ExcRequiresADNumberSpecialization());
146 return false;
147 }
148
149
150 template <typename ADNumberType, typename ScalarType, typename T>
151 bool
152 TapedDrivers<ADNumberType, ScalarType, T>::last_action_requires_retaping()
153 const
154 {
155 AssertThrow(false, ExcRequiresADNumberSpecialization());
156 return false;
157 }
158
159
160 template <typename ADNumberType, typename ScalarType, typename T>
161 void
162 TapedDrivers<ADNumberType, ScalarType, T>::remove_tape(
163 const typename Types<ADNumberType>::tape_index)
164 {
165 AssertThrow(false, ExcRequiresADNumberSpecialization());
166 }
167
168
169 template <typename ADNumberType, typename ScalarType, typename T>
170 void
171 TapedDrivers<ADNumberType, ScalarType, T>::reset(const bool)
172 {
173 AssertThrow(false, ExcRequiresADNumberSpecialization());
174 }
175
176
177 template <typename ADNumberType, typename ScalarType, typename T>
178 void
179 TapedDrivers<ADNumberType, ScalarType, T>::print(std::ostream &) const
180 {
181 AssertThrow(false, ExcRequiresADNumberSpecialization());
182 }
183
184
185 template <typename ADNumberType, typename ScalarType, typename T>
186 void
187 TapedDrivers<ADNumberType, ScalarType, T>::print_tape_stats(
188 const typename Types<ADNumberType>::tape_index,
189 std::ostream &) const
190 {
191 AssertThrow(false, ExcRequiresADNumberSpecialization());
192 }
193
194
195 template <typename ADNumberType, typename ScalarType, typename T>
196 ScalarType
197 TapedDrivers<ADNumberType, ScalarType, T>::value(
198 const typename Types<ADNumberType>::tape_index,
199 const std::vector<ScalarType> &) const
200 {
201 AssertThrow(false, ExcRequiresADNumberSpecialization());
202 return ScalarType(0.0);
203 }
204
205
206 template <typename ADNumberType, typename ScalarType, typename T>
207 void
208 TapedDrivers<ADNumberType, ScalarType, T>::gradient(
209 const typename Types<ADNumberType>::tape_index,
210 const std::vector<ScalarType> &,
211 Vector<ScalarType> &) const
212 {
213 AssertThrow(false, ExcRequiresADNumberSpecialization());
214 }
215
216
217 template <typename ADNumberType, typename ScalarType, typename T>
218 void
219 TapedDrivers<ADNumberType, ScalarType, T>::hessian(
220 const typename Types<ADNumberType>::tape_index,
221 const std::vector<ScalarType> &,
222 FullMatrix<ScalarType> &) const
223 {
224 AssertThrow(false, ExcRequiresADNumberSpecialization());
225 }
226
227
228 template <typename ADNumberType, typename ScalarType, typename T>
229 void
230 TapedDrivers<ADNumberType, ScalarType, T>::values(
231 const typename Types<ADNumberType>::tape_index,
232 const unsigned int,
233 const std::vector<ScalarType> &,
234 Vector<ScalarType> &) const
235 {
236 AssertThrow(false, ExcRequiresADNumberSpecialization());
237 }
238
239
240 template <typename ADNumberType, typename ScalarType, typename T>
241 void
242 TapedDrivers<ADNumberType, ScalarType, T>::jacobian(
243 const typename Types<ADNumberType>::tape_index,
244 const unsigned int,
245 const std::vector<ScalarType> &,
246 FullMatrix<ScalarType> &) const
247 {
248 AssertThrow(false, ExcRequiresADNumberSpecialization());
249 }
250
251
252
253# ifdef DEAL_II_WITH_ADOLC
254
255# ifndef DOXYGEN
256 // Specialization for taped ADOL-C auto-differentiable numbers.
257
258 template <typename ADNumberType>
259 TapedDrivers<ADNumberType,
260 double,
261 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
262 NumberTypes::adolc_taped>>::TapedDrivers()
263 : active_tape(Numbers<ADNumberType>::invalid_tape_index)
264 , keep_values(true)
265 , is_recording_flag(false)
266 , use_stored_taped_buffer_sizes(false)
267 , obufsize(0u)
268 , lbufsize(0u)
269 , vbufsize(0u)
270 , tbufsize(0u)
271 {}
272# endif
273
274
275 template <typename ADNumberType>
276 bool
277 TapedDrivers<ADNumberType,
278 double,
279 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
280 NumberTypes::adolc_taped>>::is_recording()
281 const
282 {
283 return is_recording_flag;
284 }
285
286
287 template <typename ADNumberType>
288 typename Types<ADNumberType>::tape_index
289 TapedDrivers<
290 ADNumberType,
291 double,
292 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
293 NumberTypes::adolc_taped>>::active_tape_index() const
294 {
295 return active_tape;
296 }
297
298
299 template <typename ADNumberType>
300 bool
301 TapedDrivers<ADNumberType,
302 double,
303 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
304 NumberTypes::adolc_taped>>::
305 keep_independent_values() const
306 {
307 return keep_values;
308 }
309
310
311 template <typename ADNumberType>
312 bool
313 TapedDrivers<ADNumberType,
314 double,
315 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
316 NumberTypes::adolc_taped>>::
317 is_registered_tape(
318 const typename Types<ADNumberType>::tape_index tape_index) const
319 {
320 // Sigh... This is a mess :-/
321 // The most succinct way to get this piece of information, would be to
322 // use the getTapeInfos() function, but would come at the expense of
323 // creating an inactive tape data within ADOL-C's global store. For
324 // hints as to why this is the way it is, see the getTapeInfos()
325 // function in
326 // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/src/tape_handling.cpp
327 // An alternative solution would be to manually access the tape data;
328 // see the removeTape() function in
329 // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/src/tape_handling.cpp
330 // with the consideration of the #defines in
331 // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/src/taping_p.h
332 // as to how this would be performed.
333 // Doing things "manually" (the second way) without creating the
334 // additional data object would be a lot more work...
335 //
336 // Both of the above solutions would be possible IF ADOL-C exposed
337 // this data object or a method to access it to the outside world.
338 // But they don't :-(
339 // Instead, what we'll have to do is get the statistics for this tape
340 // and make our own determination as to whether or not this tape exists.
341 // This effectively executes the first solution, with even more
342 // overhead! If the tape is in existence, we can assume that it should a
343 // non-zero number of dependent and independent variables. Those are
344 // stored in the zeroth and first entries of the statistics vector.
345 //
346 // But, oh wait... Surprise! This will trigger an error if the tape
347 // doesn't exist at all! So lets first check their tape info cache to
348 // see if the tape REALLY exists (i.e. has been touched, even if nothing
349 // has been written to it) before trying to access it. It'll only take
350 // an O(n) search, but at this point who really cares about efficiency?
351 //
352 // It has been suggested in
353 // https://gitlab.com/adol-c/adol-c/issues/11
354 // that a simply try-catch block around ::tapestats is the solution that
355 // we want here. Unfortunately this results in unwanted pollution of
356 // the terminal, of the form
357 // ADOL-C error: reading integer tape number 4!
358 // >>> File or directory not found! <<<
359 // , every time a query is made about a non-existent tape.
360 // So either way we have to guard that check with something more
361 // conservative so that we don't output useless messages for our users.
362 const std::vector<typename Types<ADNumberType>::tape_index>
363 registered_tape_indices = get_registered_tape_indices();
364 const auto it = std::find(registered_tape_indices.begin(),
365 registered_tape_indices.end(),
366 tape_index);
367 if (it == registered_tape_indices.end())
368 return false;
369
370 // See https://gitlab.com/adol-c/adol-c/issues/11#note_108341333
371 try
372 {
373 std::vector<std::size_t> counts(STAT_SIZE);
374 ::tapestats(tape_index, counts.data());
375 return true;
376 }
377 catch (const ::FatalError &exc)
378 {
379 return false;
380 }
381 }
382
383
384 template <typename ADNumberType>
385 void
386 TapedDrivers<ADNumberType,
387 double,
388 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
389 NumberTypes::adolc_taped>>::
390 set_tape_buffer_sizes(
391 const typename Types<ADNumberType>::tape_buffer_sizes in_obufsize,
392 const typename Types<ADNumberType>::tape_buffer_sizes in_lbufsize,
393 const typename Types<ADNumberType>::tape_buffer_sizes in_vbufsize,
394 const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize)
395 {
396 // When valid for the chosen AD number type, these values will be used
397 // the next time start_recording_operations() is called.
398 obufsize = in_obufsize;
399 lbufsize = in_lbufsize;
400 vbufsize = in_vbufsize;
401 tbufsize = in_tbufsize;
402 use_stored_taped_buffer_sizes = true;
403 }
404
405
406 template <typename ADNumberType>
407 void
408 TapedDrivers<ADNumberType,
409 double,
410 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
411 NumberTypes::adolc_taped>>::
412 start_taping(const typename Types<ADNumberType>::tape_index tape_index,
413 const bool keep_independent_values)
414 {
415 if (use_stored_taped_buffer_sizes)
416 trace_on(tape_index,
417 keep_independent_values,
418 obufsize,
419 lbufsize,
420 vbufsize,
421 tbufsize);
422 else
423 trace_on(tape_index, keep_independent_values);
424
425 // Set some other flags to their indicated / required values
426 keep_values = keep_independent_values;
427 is_recording_flag = true;
428 }
429
430
431 template <typename ADNumberType>
432 void
433 TapedDrivers<ADNumberType,
434 double,
435 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
436 NumberTypes::adolc_taped>>::
437 stop_taping(
438 const typename Types<ADNumberType>::tape_index active_tape_index,
439 const bool write_tapes_to_file)
440 {
441 if (write_tapes_to_file)
442 trace_off(active_tape_index); // Slow
443 else
444 trace_off(); // Fast(er)
445
446 // Now that we've turned tracing off, we've definitely
447 // stopped all tape recording.
448 is_recording_flag = false;
449
450 // If the keep_values flag is set, then we expect the user to use this
451 // tape immediately after recording it. There is therefore no need to
452 // invalidate it. However, there is now also no way to double-check
453 // that the newly recorded tape is indeed the active tape.
454 if (keep_independent_values() == false)
455 active_tape = Numbers<ADNumberType>::invalid_tape_index;
456 }
457
458
459 template <typename ADNumberType>
460 std::vector<typename Types<ADNumberType>::tape_index>
461 TapedDrivers<ADNumberType,
462 double,
463 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
464 NumberTypes::adolc_taped>>::
465 get_registered_tape_indices() const
466 {
467 // We've chosen to use unsigned shorts for the tape
468 // index type (a safety precaution) so we need to
469 // perform a conversion between ADOL-C's native tape
470 // index type and that chosen by us.
471 std::vector<short> registered_tape_indices_s;
472 cachedTraceTags(registered_tape_indices_s);
473
474 return std::vector<typename Types<ADNumberType>::tape_index>(
475 registered_tape_indices_s.begin(), registered_tape_indices_s.end());
476 }
477
478
479 template <typename ADNumberType>
480 void
481 TapedDrivers<ADNumberType,
482 double,
483 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
484 NumberTypes::adolc_taped>>::
485 activate_tape(const typename Types<ADNumberType>::tape_index tape_index)
486 {
487 active_tape = tape_index;
488 }
489
490
491 template <typename ADNumberType>
492 bool
493 TapedDrivers<ADNumberType,
494 double,
495 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
496 NumberTypes::adolc_taped>>::
497 requires_retaping(
498 const typename Types<ADNumberType>::tape_index tape_index) const
499 {
500 Assert(
501 is_registered_tape(tape_index) == true,
502 ExcMessage(
503 "Cannot ask for the status of a tape that has not yet been recorded and used."));
504
505 const auto it_status_tape = status.find(tape_index);
506
507 // This tape's status has not been found in the map. This could be
508 // because a non-existent tape_index has been used as an argument, or
509 // because the tape exists but has not been used (in a way that
510 // initiated a status update). For example, on can create a tape in one
511 // section of code and then query the status of all existing tapes in a
512 // completely different section of code that knows nothing about the
513 // first tape. There is no prerequisite that the first tape is ever
514 // used, and it can therefore not have a status. So, in this case
515 // there's not much we can do other than to report that the tape does
516 // not require retaping.
517 if (it_status_tape == status.end())
518 return false;
519
520 const auto status_tape = it_status_tape->second;
521
522 // See ADOL-C manual section 1.7 and comments in last paragraph of
523 // section 3.1
524 Assert(
525 status_tape < 4 && status_tape >= -2,
526 ExcMessage(
527 "The tape status is not within the range specified within the ADOL-C documentation."));
528 return (status_tape < 0);
529 }
530
531
532 template <typename ADNumberType>
533 bool
534 TapedDrivers<ADNumberType,
535 double,
536 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
537 NumberTypes::adolc_taped>>::
538 last_action_requires_retaping() const
539 {
540 return requires_retaping(active_tape);
541 }
542
543
544 template <typename ADNumberType>
545 void
546 TapedDrivers<ADNumberType,
547 double,
548 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
549 NumberTypes::adolc_taped>>::
550 remove_tape(const typename Types<ADNumberType>::tape_index tape_index)
551 {
552 Assert(is_registered_tape(tape_index),
553 ExcMessage(
554 "This tape does not exist, and therefore cannot be cleared."));
555 removeTape(tape_index, TapeRemovalType::ADOLC_REMOVE_COMPLETELY);
556 }
557
558
559 template <typename ADNumberType>
560 void
561 TapedDrivers<ADNumberType,
562 double,
563 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
564 NumberTypes::adolc_taped>>::
565 reset(const bool clear_registered_tapes)
566 {
567 active_tape = Numbers<ADNumberType>::invalid_tape_index;
568 is_recording_flag = false;
569 status.clear();
570 if (clear_registered_tapes)
571 {
572 const std::vector<typename Types<ADNumberType>::tape_index>
573 registered_tape_indices = get_registered_tape_indices();
574 for (const auto &tape_index : registered_tape_indices)
575 remove_tape(tape_index);
576 }
577 }
578
579
580 template <typename ADNumberType>
581 void
582 TapedDrivers<ADNumberType,
583 double,
584 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
585 NumberTypes::adolc_taped>>::print(std::ostream
586 &stream)
587 const
588 {
589 const std::vector<typename Types<ADNumberType>::tape_index>
590 registered_tape_indices = get_registered_tape_indices();
591 stream << "Registered tapes and their status: ";
592 auto it_registered_tape = registered_tape_indices.begin();
593 for (unsigned int i = 0; i < registered_tape_indices.size();
594 ++i, ++it_registered_tape)
595 {
596 const auto tape_index = *it_registered_tape;
597 const auto it_status_tape = status.find(tape_index);
598 Assert(it_status_tape != status.end(),
599 ExcMessage(
600 "This tape's status has not been found in the map."));
601 const auto status_tape = it_status_tape->second;
602
603 stream << tape_index << "->" << status_tape
604 << (i < (registered_tape_indices.size() - 1) ? "," : "");
605 }
606 stream << '\n';
607
608 stream << "Keep values? " << keep_independent_values() << '\n';
609 stream << "Use stored tape buffer sizes? "
610 << use_stored_taped_buffer_sizes << '\n';
611 }
612
613
614 template <typename ADNumberType>
615 void
616 TapedDrivers<ADNumberType,
617 double,
618 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
619 NumberTypes::adolc_taped>>::
620 print_tape_stats(
621 const typename Types<ADNumberType>::tape_index tape_index,
622 std::ostream & stream) const
623 {
624 // See ADOL-C manual section 2.1
625 // and adolc/taping.h
626 std::vector<std::size_t> counts(STAT_SIZE);
627 ::tapestats(tape_index, counts.data());
628 Assert(counts.size() >= 18, ExcInternalError());
629 stream
630 << "Tape index: " << tape_index << '\n'
631 << "Number of independent variables: " << counts[0] << '\n'
632 << "Number of dependent variables: " << counts[1] << '\n'
633 << "Max number of live, active variables: " << counts[2] << '\n'
634 << "Size of taylor stack (number of overwrites): " << counts[3] << '\n'
635 << "Operations buffer size: " << counts[4] << '\n'
636 << "Total number of recorded operations: " << counts[5] << '\n'
637 << "Operations file written or not: " << counts[6] << '\n'
638 << "Overall number of locations: " << counts[7] << '\n'
639 << "Locations file written or not: " << counts[8] << '\n'
640 << "Overall number of values: " << counts[9] << '\n'
641 << "Values file written or not: " << counts[10] << '\n'
642 << "Locations buffer size: " << counts[11] << '\n'
643 << "Values buffer size: " << counts[12] << '\n'
644 << "Taylor buffer size: " << counts[13] << '\n'
645 << "Number of eq_*_prod for sparsity pattern: " << counts[14] << '\n'
646 << "Use of 'min_op', deferred to 'abs_op' for piecewise calculations: "
647 << counts[15] << '\n'
648 << "Number of 'abs' calls that can switch branch: " << counts[16]
649 << '\n'
650 << "Number of parameters (doubles) interchangeable without retaping: "
651 << counts[17] << '\n'
652 << std::flush;
653 }
654
655
656 template <typename ADNumberType>
657 typename TapedDrivers<
658 ADNumberType,
659 double,
660 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
661 NumberTypes::adolc_taped>>::scalar_type
662 TapedDrivers<ADNumberType,
663 double,
664 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
665 NumberTypes::adolc_taped>>::
666 value(const typename Types<ADNumberType>::tape_index active_tape_index,
667 const std::vector<scalar_type> &independent_variables) const
668 {
669 Assert(is_registered_tape(active_tape_index),
670 ExcMessage("This tape has not yet been recorded."));
671
672 scalar_type value = 0.0;
673
674 status[active_tape_index] =
675 ::function(active_tape_index,
676 1, // Only one dependent variable
677 independent_variables.size(),
678 const_cast<double *>(independent_variables.data()),
679 &value);
680
681 return value;
682 }
683
684
685 template <typename ADNumberType>
686 void
687 TapedDrivers<ADNumberType,
688 double,
689 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
690 NumberTypes::adolc_taped>>::
691 gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
692 const std::vector<scalar_type> &independent_variables,
693 Vector<scalar_type> & gradient) const
694 {
695 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
696 1,
697 ExcSupportedDerivativeLevels(
698 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
699 1));
700 Assert(gradient.size() == independent_variables.size(),
701 ExcDimensionMismatch(gradient.size(),
702 independent_variables.size()));
703 Assert(is_registered_tape(active_tape_index),
704 ExcMessage("This tape has not yet been recorded."));
705
706 // Note: ADOL-C's ::gradient function expects a *double as the last
707 // parameter. Here we take advantage of the fact that the data in the
708 // Vector class is aligned (e.g. stored as an Array)
709 status[active_tape_index] =
710 ::gradient(active_tape_index,
711 independent_variables.size(),
712 const_cast<scalar_type *>(independent_variables.data()),
713 gradient.data());
714 }
715
716
717 template <typename ADNumberType>
718 void
719 TapedDrivers<ADNumberType,
720 double,
721 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
722 NumberTypes::adolc_taped>>::
723 hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
724 const std::vector<scalar_type> &independent_variables,
725 FullMatrix<scalar_type> & hessian) const
726 {
727 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
728 2,
729 ExcSupportedDerivativeLevels(
730 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
731 2));
732 Assert(hessian.m() == independent_variables.size(),
733 ExcDimensionMismatch(hessian.m(), independent_variables.size()));
734 Assert(hessian.n() == independent_variables.size(),
735 ExcDimensionMismatch(hessian.n(), independent_variables.size()));
736 Assert(is_registered_tape(active_tape_index),
737 ExcMessage("This tape has not yet been recorded."));
738
739 const unsigned int n_independent_variables = independent_variables.size();
740 std::vector<scalar_type *> H(n_independent_variables);
741 for (unsigned int i = 0; i < n_independent_variables; ++i)
742 H[i] = &hessian[i][0];
743
744 status[active_tape_index] =
745 ::hessian(active_tape_index,
746 n_independent_variables,
747 const_cast<scalar_type *>(independent_variables.data()),
748 H.data());
749
750 // ADOL-C builds only the lower-triangular part of the
751 // symmetric Hessian, so we should copy the relevant
752 // entries into the upper triangular part.
753 for (unsigned int i = 0; i < n_independent_variables; ++i)
754 for (unsigned int j = 0; j < i; ++j)
755 hessian[j][i] = hessian[i][j]; // Symmetry
756 }
757
758
759 template <typename ADNumberType>
760 void
761 TapedDrivers<ADNumberType,
762 double,
763 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
764 NumberTypes::adolc_taped>>::
765 values(const typename Types<ADNumberType>::tape_index active_tape_index,
766 const unsigned int n_dependent_variables,
767 const std::vector<scalar_type> &independent_variables,
768 Vector<scalar_type> & values) const
769 {
770 Assert(values.size() == n_dependent_variables,
771 ExcDimensionMismatch(values.size(), n_dependent_variables));
772 Assert(is_registered_tape(active_tape_index),
773 ExcMessage("This tape has not yet been recorded."));
774
775 // Note: ADOL-C's ::function function expects a *double as the last
776 // parameter. Here we take advantage of the fact that the data in the
777 // Vector class is aligned (e.g. stored as an Array)
778 status[active_tape_index] =
779 ::function(active_tape_index,
780 n_dependent_variables,
781 independent_variables.size(),
782 const_cast<scalar_type *>(independent_variables.data()),
783 values.data());
784 }
785
786
787 template <typename ADNumberType>
788 void
789 TapedDrivers<ADNumberType,
790 double,
791 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
792 NumberTypes::adolc_taped>>::
793 jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
794 const unsigned int n_dependent_variables,
795 const std::vector<scalar_type> &independent_variables,
796 FullMatrix<scalar_type> & jacobian) const
797 {
798 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
799 1,
800 ExcSupportedDerivativeLevels(
801 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
802 1));
803 Assert(jacobian.m() == n_dependent_variables,
804 ExcDimensionMismatch(jacobian.m(), n_dependent_variables));
805 Assert(jacobian.n() == independent_variables.size(),
806 ExcDimensionMismatch(jacobian.n(), independent_variables.size()));
807 Assert(is_registered_tape(active_tape_index),
808 ExcMessage("This tape has not yet been recorded."));
809
810 std::vector<scalar_type *> J(n_dependent_variables);
811 for (unsigned int i = 0; i < n_dependent_variables; ++i)
812 J[i] = &jacobian[i][0];
813
814 status[active_tape_index] = ::jacobian(active_tape_index,
815 n_dependent_variables,
816 independent_variables.size(),
817 independent_variables.data(),
818 J.data());
819 }
820
821# else
822
823 // Specialization for taped ADOL-C auto-differentiable numbers.
824
825 template <typename ADNumberType>
826 bool
827 TapedDrivers<ADNumberType,
828 double,
829 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
830 NumberTypes::adolc_taped>>::is_recording()
831 const
832 {
833 AssertThrow(false, ExcRequiresADOLC());
834 return false;
835 }
836
837
838 template <typename ADNumberType>
839 typename Types<ADNumberType>::tape_index
840 TapedDrivers<
841 ADNumberType,
842 double,
843 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
844 NumberTypes::adolc_taped>>::active_tape_index() const
845 {
846 AssertThrow(false, ExcRequiresADOLC());
847 return Numbers<ADNumberType>::invalid_tape_index;
848 }
849
850
851 template <typename ADNumberType>
852 bool
853 TapedDrivers<ADNumberType,
854 double,
855 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
856 NumberTypes::adolc_taped>>::
857 keep_independent_values() const
858 {
859 AssertThrow(false, ExcRequiresADOLC());
860 return false;
861 }
862
863
864 template <typename ADNumberType>
865 bool
866 TapedDrivers<ADNumberType,
867 double,
868 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
869 NumberTypes::adolc_taped>>::
870 is_registered_tape(const typename Types<ADNumberType>::tape_index) const
871 {
872 AssertThrow(false, ExcRequiresADOLC());
873 return false;
874 }
875
876
877 template <typename ADNumberType>
878 void
879 TapedDrivers<ADNumberType,
880 double,
881 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
882 NumberTypes::adolc_taped>>::
883 set_tape_buffer_sizes(
884 const typename Types<ADNumberType>::tape_buffer_sizes,
885 const typename Types<ADNumberType>::tape_buffer_sizes,
886 const typename Types<ADNumberType>::tape_buffer_sizes,
887 const typename Types<ADNumberType>::tape_buffer_sizes)
888 {
889 AssertThrow(false, ExcRequiresADOLC());
890 }
891
892
893 template <typename ADNumberType>
894 void
895 TapedDrivers<ADNumberType,
896 double,
897 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
898 NumberTypes::adolc_taped>>::
899 start_taping(const typename Types<ADNumberType>::tape_index, const bool)
900 {
901 AssertThrow(false, ExcRequiresADOLC());
902 }
903
904
905 template <typename ADNumberType>
906 void
907 TapedDrivers<ADNumberType,
908 double,
909 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
910 NumberTypes::adolc_taped>>::
911 stop_taping(const typename Types<ADNumberType>::tape_index, const bool)
912 {
913 AssertThrow(false, ExcRequiresADOLC());
914 }
915
916
917 template <typename ADNumberType>
918 std::vector<typename Types<ADNumberType>::tape_index>
919 TapedDrivers<ADNumberType,
920 double,
921 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
922 NumberTypes::adolc_taped>>::
923 get_registered_tape_indices() const
924 {
925 AssertThrow(false, ExcRequiresADOLC());
926 return std::vector<typename Types<ADNumberType>::tape_index>();
927 }
928
929
930 template <typename ADNumberType>
931 void
932 TapedDrivers<ADNumberType,
933 double,
934 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
935 NumberTypes::adolc_taped>>::
936 activate_tape(const typename Types<ADNumberType>::tape_index)
937 {
938 AssertThrow(false, ExcRequiresADOLC());
939 }
940
941
942 template <typename ADNumberType>
943 bool
944 TapedDrivers<ADNumberType,
945 double,
946 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
947 NumberTypes::adolc_taped>>::
948 requires_retaping(const typename Types<ADNumberType>::tape_index) const
949 {
950 AssertThrow(false, ExcRequiresADOLC());
951 return false;
952 }
953
954
955 template <typename ADNumberType>
956 bool
957 TapedDrivers<ADNumberType,
958 double,
959 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
960 NumberTypes::adolc_taped>>::
961 last_action_requires_retaping() const
962 {
963 AssertThrow(false, ExcRequiresADOLC());
964 return false;
965 }
966
967
968 template <typename ADNumberType>
969 void
970 TapedDrivers<ADNumberType,
971 double,
972 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
973 NumberTypes::adolc_taped>>::
974 remove_tape(const typename Types<ADNumberType>::tape_index)
975 {
976 AssertThrow(false, ExcRequiresADOLC());
977 }
978
979
980 template <typename ADNumberType>
981 void
982 TapedDrivers<ADNumberType,
983 double,
984 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
985 NumberTypes::adolc_taped>>::reset(const bool)
986 {
987 AssertThrow(false, ExcRequiresADOLC());
988 }
989
990
991 template <typename ADNumberType>
992 void
993 TapedDrivers<ADNumberType,
994 double,
995 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
996 NumberTypes::adolc_taped>>::print(std::ostream
997 &) const
998 {
999 AssertThrow(false, ExcRequiresADOLC());
1000 }
1001
1002
1003 template <typename ADNumberType>
1004 void
1005 TapedDrivers<ADNumberType,
1006 double,
1007 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1008 NumberTypes::adolc_taped>>::
1009 print_tape_stats(const typename Types<ADNumberType>::tape_index,
1010 std::ostream &) const
1011 {
1012 AssertThrow(false, ExcRequiresADOLC());
1013 }
1014
1015
1016 template <typename ADNumberType>
1017 typename TapedDrivers<
1018 ADNumberType,
1019 double,
1020 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1021 NumberTypes::adolc_taped>>::scalar_type
1022 TapedDrivers<ADNumberType,
1023 double,
1024 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1025 NumberTypes::adolc_taped>>::
1026 value(const typename Types<ADNumberType>::tape_index,
1027 const std::vector<scalar_type> &) const
1028 {
1029 AssertThrow(false, ExcRequiresADOLC());
1030 return 0.0;
1031 }
1032
1033
1034 template <typename ADNumberType>
1035 void
1036 TapedDrivers<ADNumberType,
1037 double,
1038 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1039 NumberTypes::adolc_taped>>::
1040 gradient(const typename Types<ADNumberType>::tape_index,
1041 const std::vector<scalar_type> &,
1042 Vector<scalar_type> &) const
1043 {
1044 AssertThrow(false, ExcRequiresADOLC());
1045 }
1046
1047
1048 template <typename ADNumberType>
1049 void
1050 TapedDrivers<ADNumberType,
1051 double,
1052 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1053 NumberTypes::adolc_taped>>::
1054 hessian(const typename Types<ADNumberType>::tape_index,
1055 const std::vector<scalar_type> &,
1056 FullMatrix<scalar_type> &) const
1057 {
1058 AssertThrow(false, ExcRequiresADOLC());
1059 }
1060
1061
1062 template <typename ADNumberType>
1063 void
1064 TapedDrivers<ADNumberType,
1065 double,
1066 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1067 NumberTypes::adolc_taped>>::
1068 values(const typename Types<ADNumberType>::tape_index,
1069 const unsigned int,
1070 const std::vector<scalar_type> &,
1071 Vector<scalar_type> &) const
1072 {
1073 AssertThrow(false, ExcRequiresADOLC());
1074 }
1075
1076
1077 template <typename ADNumberType>
1078 void
1079 TapedDrivers<ADNumberType,
1080 double,
1081 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1082 NumberTypes::adolc_taped>>::
1083 jacobian(const typename Types<ADNumberType>::tape_index,
1084 const unsigned int,
1085 const std::vector<scalar_type> &,
1086 FullMatrix<scalar_type> &) const
1087 {
1088 AssertThrow(false, ExcRequiresADOLC());
1089 }
1090
1091# endif // DEAL_II_WITH_ADOLC
1092
1093
1094 // Specialization for ADOL-C taped numbers. It is expected that the
1095 // scalar return type for this class is a float.
1096
1097 template <typename ADNumberType>
1098 bool
1099 TapedDrivers<ADNumberType,
1100 float,
1101 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1102 NumberTypes::adolc_taped>>::is_recording()
1103 const
1104 {
1105 // ADOL-C only supports 'double', not 'float', so we can forward to
1106 // the 'double' implementation of this function
1107 return taped_driver.is_recording();
1108 }
1109
1110
1111 template <typename ADNumberType>
1112 typename Types<ADNumberType>::tape_index
1113 TapedDrivers<
1114 ADNumberType,
1115 float,
1116 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1117 NumberTypes::adolc_taped>>::active_tape_index() const
1118 {
1119 // ADOL-C only supports 'double', not 'float', so we can forward to
1120 // the 'double' implementation of this function
1121 return taped_driver.active_tape_index();
1122 }
1123
1124
1125 template <typename ADNumberType>
1126 bool
1127 TapedDrivers<ADNumberType,
1128 float,
1129 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1130 NumberTypes::adolc_taped>>::
1131 keep_independent_values() const
1132 {
1133 return taped_driver.keep_independent_values();
1134 }
1135
1136
1137 template <typename ADNumberType>
1138 bool
1139 TapedDrivers<ADNumberType,
1140 float,
1141 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1142 NumberTypes::adolc_taped>>::
1143 is_registered_tape(
1144 const typename Types<ADNumberType>::tape_index tape_index) const
1145 {
1146 // ADOL-C only supports 'double', not 'float', so we can forward to
1147 // the 'double' implementation of this function
1148 return taped_driver.is_registered_tape(tape_index);
1149 }
1150
1151
1152 template <typename ADNumberType>
1153 void
1154 TapedDrivers<ADNumberType,
1155 float,
1156 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1157 NumberTypes::adolc_taped>>::
1158 set_tape_buffer_sizes(
1159 const typename Types<ADNumberType>::tape_buffer_sizes obufsize,
1160 const typename Types<ADNumberType>::tape_buffer_sizes lbufsize,
1161 const typename Types<ADNumberType>::tape_buffer_sizes vbufsize,
1162 const typename Types<ADNumberType>::tape_buffer_sizes tbufsize)
1163 {
1164 // ADOL-C only supports 'double', not 'float', so we can forward to
1165 // the 'double' implementation of this function
1166 taped_driver.set_tape_buffer_sizes(obufsize,
1167 lbufsize,
1168 vbufsize,
1169 tbufsize);
1170 }
1171
1172
1173 template <typename ADNumberType>
1174 void
1175 TapedDrivers<ADNumberType,
1176 float,
1177 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1178 NumberTypes::adolc_taped>>::
1179 start_taping(const typename Types<ADNumberType>::tape_index tape_index,
1180 const bool keep_independent_values)
1181 {
1182 // ADOL-C only supports 'double', not 'float', so we can forward to
1183 // the 'double' implementation of this function
1184 taped_driver.start_taping(tape_index, keep_independent_values);
1185 }
1186
1187
1188 template <typename ADNumberType>
1189 void
1190 TapedDrivers<ADNumberType,
1191 float,
1192 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1193 NumberTypes::adolc_taped>>::
1194 stop_taping(
1195 const typename Types<ADNumberType>::tape_index active_tape_index,
1196 const bool write_tapes_to_file)
1197 {
1198 // ADOL-C only supports 'double', not 'float', so we can forward to
1199 // the 'double' implementation of this function
1200 taped_driver.stop_taping(active_tape_index, write_tapes_to_file);
1201 }
1202
1203
1204 template <typename ADNumberType>
1205 std::vector<typename Types<ADNumberType>::tape_index>
1206 TapedDrivers<ADNumberType,
1207 float,
1208 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1209 NumberTypes::adolc_taped>>::
1210 get_registered_tape_indices() const
1211 {
1212 return taped_driver.get_registered_tape_indices();
1213 }
1214
1215
1216 template <typename ADNumberType>
1217 void
1218 TapedDrivers<ADNumberType,
1219 float,
1220 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1221 NumberTypes::adolc_taped>>::
1222 activate_tape(const typename Types<ADNumberType>::tape_index tape_index)
1223 {
1224 taped_driver.activate_tape(tape_index);
1225 }
1226
1227
1228 template <typename ADNumberType>
1229 bool
1230 TapedDrivers<ADNumberType,
1231 float,
1232 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1233 NumberTypes::adolc_taped>>::
1234 requires_retaping(
1235 const typename Types<ADNumberType>::tape_index tape_index) const
1236 {
1237 return taped_driver.requires_retaping(tape_index);
1238 }
1239
1240
1241 template <typename ADNumberType>
1242 bool
1243 TapedDrivers<ADNumberType,
1244 float,
1245 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1246 NumberTypes::adolc_taped>>::
1247 last_action_requires_retaping() const
1248 {
1249 return taped_driver.last_action_requires_retaping();
1250 }
1251
1252
1253 template <typename ADNumberType>
1254 void
1255 TapedDrivers<ADNumberType,
1256 float,
1257 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1258 NumberTypes::adolc_taped>>::
1259 remove_tape(const typename Types<ADNumberType>::tape_index tape_index)
1260 {
1261 taped_driver.remove_tape(tape_index);
1262 }
1263
1264
1265 template <typename ADNumberType>
1266 void
1267 TapedDrivers<ADNumberType,
1268 float,
1269 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1270 NumberTypes::adolc_taped>>::
1271 reset(const bool clear_registered_tapes)
1272 {
1273 taped_driver.reset(clear_registered_tapes);
1274 }
1275
1276
1277 template <typename ADNumberType>
1278 void
1279 TapedDrivers<ADNumberType,
1280 float,
1281 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1282 NumberTypes::adolc_taped>>::print(std::ostream
1283 &stream)
1284 const
1285 {
1286 taped_driver.print(stream);
1287 }
1288
1289
1290 template <typename ADNumberType>
1291 void
1292 TapedDrivers<ADNumberType,
1293 float,
1294 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1295 NumberTypes::adolc_taped>>::
1296 print_tape_stats(
1297 const typename Types<ADNumberType>::tape_index tape_index,
1298 std::ostream & stream) const
1299 {
1300 // ADOL-C only supports 'double', not 'float', so we can forward to
1301 // the 'double' implementation of this function
1302 taped_driver.print_tape_stats(tape_index, stream);
1303 }
1304
1305
1306 template <typename ADNumberType>
1307 typename TapedDrivers<
1308 ADNumberType,
1309 float,
1310 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1311 NumberTypes::adolc_taped>>::scalar_type
1312 TapedDrivers<ADNumberType,
1313 float,
1314 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1315 NumberTypes::adolc_taped>>::
1316 value(const typename Types<ADNumberType>::tape_index active_tape_index,
1317 const std::vector<scalar_type> &independent_variables) const
1318 {
1319 // ADOL-C only supports 'double', not 'float', so we can forward to
1320 // the 'double' implementation of this function
1321 return taped_driver.value(active_tape_index,
1322 vector_float_to_double(independent_variables));
1323 }
1324
1325
1326 template <typename ADNumberType>
1327 void
1328 TapedDrivers<ADNumberType,
1329 float,
1330 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1331 NumberTypes::adolc_taped>>::
1332 gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
1333 const std::vector<scalar_type> &independent_variables,
1334 Vector<scalar_type> & gradient) const
1335 {
1336 Vector<double> gradient_double(gradient.size());
1337 // ADOL-C only supports 'double', not 'float', so we can forward to
1338 // the 'double' implementation of this function
1339 taped_driver.gradient(active_tape_index,
1340 vector_float_to_double(independent_variables),
1341 gradient_double);
1342 gradient = gradient_double;
1343 }
1344
1345
1346 template <typename ADNumberType>
1347 void
1348 TapedDrivers<ADNumberType,
1349 float,
1350 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1351 NumberTypes::adolc_taped>>::
1352 hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
1353 const std::vector<scalar_type> &independent_variables,
1354 FullMatrix<scalar_type> & hessian) const
1355 {
1356 FullMatrix<double> hessian_double(hessian.m(), hessian.n());
1357 // ADOL-C only supports 'double', not 'float', so we can forward to
1358 // the 'double' implementation of this function
1359 taped_driver.hessian(active_tape_index,
1360 vector_float_to_double(independent_variables),
1361 hessian_double);
1362 hessian = hessian_double;
1363 }
1364
1365
1366 template <typename ADNumberType>
1367 void
1368 TapedDrivers<ADNumberType,
1369 float,
1370 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1371 NumberTypes::adolc_taped>>::
1372 values(const typename Types<ADNumberType>::tape_index active_tape_index,
1373 const unsigned int n_dependent_variables,
1374 const std::vector<scalar_type> &independent_variables,
1375 Vector<scalar_type> & values) const
1376 {
1377 Vector<double> values_double(values.size());
1378 // ADOL-C only supports 'double', not 'float', so we can forward to
1379 // the 'double' implementation of this function
1380 taped_driver.values(active_tape_index,
1381 n_dependent_variables,
1382 vector_float_to_double(independent_variables),
1383 values_double);
1384 values = values_double;
1385 }
1386
1387
1388 template <typename ADNumberType>
1389 void
1390 TapedDrivers<ADNumberType,
1391 float,
1392 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1393 NumberTypes::adolc_taped>>::
1394 jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
1395 const unsigned int n_dependent_variables,
1396 const std::vector<scalar_type> &independent_variables,
1397 FullMatrix<scalar_type> & jacobian) const
1398 {
1399 FullMatrix<double> jacobian_double(jacobian.m(), jacobian.n());
1400 // ADOL-C only supports 'double', not 'float', so we can forward to
1401 // the 'double' implementation of this function
1402 taped_driver.jacobian(active_tape_index,
1403 n_dependent_variables,
1404 vector_float_to_double(independent_variables),
1405 jacobian_double);
1406 jacobian = jacobian_double;
1407 }
1408
1409
1410# ifndef DOXYGEN
1411 template <typename ADNumberType>
1412 std::vector<double>
1413 TapedDrivers<ADNumberType,
1414 float,
1415 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1416 NumberTypes::adolc_taped>>::
1417 vector_float_to_double(const std::vector<float> &in) const
1418 {
1419 std::vector<double> out(in.size());
1420 std::copy(in.begin(), in.end(), out.begin());
1421 return out;
1422 }
1423# endif
1424
1425
1426 // ------------- TapelessDrivers -------------
1427
1428
1429 template <typename ADNumberType, typename ScalarType, typename T>
1430 void
1431 TapelessDrivers<ADNumberType, ScalarType, T>::initialize_global_environment(
1432 const unsigned int)
1433 {
1434 AssertThrow(false, ExcRequiresADNumberSpecialization());
1435 }
1436
1437 template <typename ADNumberType, typename ScalarType, typename T>
1438 void
1439 TapelessDrivers<ADNumberType, ScalarType, T>::
1440 allow_dependent_variable_marking()
1441 {
1442 AssertThrow(false, ExcRequiresADNumberSpecialization());
1443 }
1444
1445 template <typename ADNumberType, typename ScalarType, typename T>
1446 void
1447 TapelessDrivers<ADNumberType, ScalarType, T>::
1448 prevent_dependent_variable_marking()
1449 {
1450 AssertThrow(false, ExcRequiresADNumberSpecialization());
1451 }
1452
1453 template <typename ADNumberType, typename ScalarType, typename T>
1454 bool
1455 TapelessDrivers<ADNumberType, ScalarType, T>::
1456 is_dependent_variable_marking_allowed() const
1457 {
1458 AssertThrow(false, ExcRequiresADNumberSpecialization());
1459 return false;
1460 }
1461
1462
1463 template <typename ADNumberType, typename ScalarType, typename T>
1464 ScalarType
1465 TapelessDrivers<ADNumberType, ScalarType, T>::value(
1466 const std::vector<ADNumberType> &) const
1467 {
1468 AssertThrow(false, ExcRequiresADNumberSpecialization());
1469 return ScalarType(0.0);
1470 }
1471
1472
1473 template <typename ADNumberType, typename ScalarType, typename T>
1474 void
1475 TapelessDrivers<ADNumberType, ScalarType, T>::gradient(
1476 const std::vector<ADNumberType> &,
1477 const std::vector<ADNumberType> &,
1478 Vector<ScalarType> &) const
1479 {
1480 AssertThrow(false, ExcRequiresADNumberSpecialization());
1481 }
1482
1483
1484 template <typename ADNumberType, typename ScalarType, typename T>
1485 void
1486 TapelessDrivers<ADNumberType, ScalarType, T>::hessian(
1487 const std::vector<ADNumberType> &,
1488 const std::vector<ADNumberType> &,
1489 FullMatrix<ScalarType> &) const
1490 {
1491 AssertThrow(false, ExcRequiresADNumberSpecialization());
1492 }
1493
1494
1495 template <typename ADNumberType, typename ScalarType, typename T>
1496 void
1497 TapelessDrivers<ADNumberType, ScalarType, T>::values(
1498 const std::vector<ADNumberType> &,
1499 Vector<ScalarType> &) const
1500 {
1501 AssertThrow(false, ExcRequiresADNumberSpecialization());
1502 }
1503
1504
1505 template <typename ADNumberType, typename ScalarType, typename T>
1506 void
1507 TapelessDrivers<ADNumberType, ScalarType, T>::jacobian(
1508 const std::vector<ADNumberType> &,
1509 const std::vector<ADNumberType> &,
1510 FullMatrix<ScalarType> &) const
1511 {
1512 AssertThrow(false, ExcRequiresADNumberSpecialization());
1513 }
1514
1515
1516 namespace internal
1517 {
1522 template <typename ADNumberType>
1523 std::enable_if_t<!(ADNumberTraits<ADNumberType>::type_code ==
1524 NumberTypes::sacado_rad ||
1525 ADNumberTraits<ADNumberType>::type_code ==
1526 NumberTypes::sacado_rad_dfad)>
1527 reverse_mode_dependent_variable_activation(ADNumberType &)
1528 {}
1529
1530# ifdef DEAL_II_TRILINOS_WITH_SACADO
1531
1532
1541 template <typename ADNumberType>
1542 std::enable_if_t<
1543 ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_rad ||
1544 ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_rad_dfad>
1545 reverse_mode_dependent_variable_activation(
1546 ADNumberType &dependent_variable)
1547 {
1548 // Compute all gradients (adjoints) for this
1549 // reverse-mode Sacado dependent variable.
1550 // For reverse-mode Sacado numbers it is necessary to broadcast to
1551 // all independent variables that it is time to compute gradients.
1552 // For one dependent variable one would just need to call
1553 // ADNumberType::Gradcomp(), but since we have a more
1554 // generic implementation for vectors of dependent variables
1555 // (vector mode) we default to the complex case.
1556 ADNumberType::Outvar_Gradcomp(dependent_variable);
1557 }
1558
1559# endif
1560
1561
1566 template <typename ADNumberType>
1567 std::enable_if_t<!(ADNumberTraits<ADNumberType>::type_code ==
1568 NumberTypes::adolc_tapeless)>
1569 configure_tapeless_mode(const unsigned int)
1570 {}
1571
1572
1573# ifdef DEAL_II_WITH_ADOLC
1574
1575
1583 template <typename ADNumberType>
1584 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1585 NumberTypes::adolc_tapeless>
1586 configure_tapeless_mode(const unsigned int n_directional_derivatives)
1587 {
1588# ifdef DEAL_II_ADOLC_WITH_TAPELESS_REFCOUNTING
1589 // See ADOL-C manual section 7.1
1590 //
1591 // NOTE: It is critical that this is done for tapeless mode BEFORE
1592 // any adtl::adouble are created. If this is not done, then we see
1593 // this scary warning:
1594 //
1595 // "
1596 // ADOL-C Warning: Tapeless: Setting numDir could change memory
1597 // allocation of derivatives in existing adoubles and may lead to
1598 // erroneous results or memory corruption
1599 // "
1600 //
1601 // So we use this dummy function to configure this setting before
1602 // we create and initialize our class data
1603 const std::size_t n_live_variables = adtl::refcounter::getNumLiveVar();
1604 if (n_live_variables == 0)
1605 {
1606 adtl::setNumDir(n_directional_derivatives);
1607 }
1608 else
1609 {
1610 // So there are some live active variables floating around. Here we
1611 // check if we ask to increase the number of computable
1612 // directional derivatives. If this really is necessary then it's
1613 // absolutely vital that there exist no live variables before doing
1614 // so.
1615 const std::size_t n_set_directional_derivatives = adtl::getNumDir();
1616 if (n_directional_derivatives > n_set_directional_derivatives)
1617 AssertThrow(
1618 n_live_variables == 0,
1619 ExcMessage(
1620 "There are currently " + std::to_string(n_live_variables) +
1621 " live "
1622 "adtl::adouble variables in existence. They currently "
1623 "assume " +
1624 std::to_string(n_set_directional_derivatives) +
1625 " directional derivatives "
1626 "but you wish to increase this to " +
1627 std::to_string(n_directional_derivatives) +
1628 ". \n"
1629 "To safely change (or more specifically in this case, "
1630 "increase) the number of directional derivatives, there "
1631 "must be no tapeless doubles in local/global scope."));
1632 }
1633# else
1634 // If ADOL-C is not configured with tapeless number reference counting
1635 // then there is no way that we can guarantee that the following call
1636 // is safe. No comment... :-/
1637 adtl::setNumDir(n_directional_derivatives);
1638# endif
1639 }
1640
1641# else // DEAL_II_WITH_ADOLC
1642
1643 template <typename ADNumberType>
1644 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1645 NumberTypes::adolc_tapeless>
1646 configure_tapeless_mode(const unsigned int /*n_directional_derivatives*/)
1647 {
1648 AssertThrow(false, ExcRequiresADOLC());
1649 }
1650
1651# endif
1652
1653 } // namespace internal
1654
1655
1656
1657 // Specialization for auto-differentiable numbers that use
1658 // reverse mode to compute the first derivatives (and, if supported,
1659 // forward mode for the second).
1660
1661# ifndef DOXYGEN
1662 template <typename ADNumberType, typename ScalarType>
1663 TapelessDrivers<
1664 ADNumberType,
1665 ScalarType,
1666 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1667 NumberTypes::sacado_rad ||
1668 ADNumberTraits<ADNumberType>::type_code ==
1669 NumberTypes::sacado_rad_dfad>>::TapelessDrivers()
1670 : dependent_variable_marking_safe(false)
1671 {}
1672# endif
1673
1674
1675 template <typename ADNumberType, typename ScalarType>
1676 void
1677 TapelessDrivers<ADNumberType,
1678 ScalarType,
1679 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1680 NumberTypes::sacado_rad ||
1681 ADNumberTraits<ADNumberType>::type_code ==
1682 NumberTypes::sacado_rad_dfad>>::
1683 initialize_global_environment(const unsigned int n_independent_variables)
1684 {
1685 internal::configure_tapeless_mode<ADNumberType>(n_independent_variables);
1686 }
1687
1688
1689 template <typename ADNumberType, typename ScalarType>
1690 void
1691 TapelessDrivers<
1692 ADNumberType,
1693 ScalarType,
1694 std::enable_if_t<
1695 ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_rad ||
1696 ADNumberTraits<ADNumberType>::type_code ==
1697 NumberTypes::sacado_rad_dfad>>::allow_dependent_variable_marking()
1698 {
1699 dependent_variable_marking_safe = true;
1700 }
1701
1702
1703 template <typename ADNumberType, typename ScalarType>
1704 void
1705 TapelessDrivers<
1706 ADNumberType,
1707 ScalarType,
1708 std::enable_if_t<
1709 ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_rad ||
1710 ADNumberTraits<ADNumberType>::type_code ==
1711 NumberTypes::sacado_rad_dfad>>::prevent_dependent_variable_marking()
1712 {
1713 dependent_variable_marking_safe = false;
1714 }
1715
1716
1717 template <typename ADNumberType, typename ScalarType>
1718 bool
1719 TapelessDrivers<ADNumberType,
1720 ScalarType,
1721 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1722 NumberTypes::sacado_rad ||
1723 ADNumberTraits<ADNumberType>::type_code ==
1724 NumberTypes::sacado_rad_dfad>>::
1725 is_dependent_variable_marking_allowed() const
1726 {
1727 return dependent_variable_marking_safe;
1728 }
1729
1730
1731 template <typename ADNumberType, typename ScalarType>
1732 ScalarType
1733 TapelessDrivers<ADNumberType,
1734 ScalarType,
1735 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1736 NumberTypes::sacado_rad ||
1737 ADNumberTraits<ADNumberType>::type_code ==
1738 NumberTypes::sacado_rad_dfad>>::
1739 value(const std::vector<ADNumberType> &dependent_variables) const
1740 {
1741 Assert(dependent_variables.size() == 1,
1742 ExcDimensionMismatch(dependent_variables.size(), 1));
1743 return ADNumberTraits<ADNumberType>::get_scalar_value(
1744 dependent_variables[0]);
1745 }
1746
1747
1748 template <typename ADNumberType, typename ScalarType>
1749 void
1750 TapelessDrivers<ADNumberType,
1751 ScalarType,
1752 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1753 NumberTypes::sacado_rad ||
1754 ADNumberTraits<ADNumberType>::type_code ==
1755 NumberTypes::sacado_rad_dfad>>::
1756 gradient(const std::vector<ADNumberType> &independent_variables,
1757 const std::vector<ADNumberType> &dependent_variables,
1758 Vector<ScalarType> & gradient) const
1759 {
1760 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
1761 1,
1762 ExcSupportedDerivativeLevels(
1763 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
1764 1));
1765 Assert(dependent_variables.size() == 1,
1766 ExcDimensionMismatch(dependent_variables.size(), 1));
1767 Assert(gradient.size() == independent_variables.size(),
1768 ExcDimensionMismatch(gradient.size(),
1769 independent_variables.size()));
1770
1771 // In reverse mode, the gradients are computed from the
1772 // independent variables (i.e. the adjoint)
1773 internal::reverse_mode_dependent_variable_activation(
1774 const_cast<ADNumberType &>(dependent_variables[0]));
1775 const std::size_t n_independent_variables = independent_variables.size();
1776 for (unsigned int i = 0; i < n_independent_variables; ++i)
1777 gradient[i] = internal::NumberType<ScalarType>::value(
1778 ADNumberTraits<ADNumberType>::get_directional_derivative(
1779 independent_variables[i], 0 /*This number doesn't really matter*/));
1780 }
1781
1782
1783 template <typename ADNumberType, typename ScalarType>
1784 void
1785 TapelessDrivers<ADNumberType,
1786 ScalarType,
1787 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1788 NumberTypes::sacado_rad ||
1789 ADNumberTraits<ADNumberType>::type_code ==
1790 NumberTypes::sacado_rad_dfad>>::
1791 hessian(const std::vector<ADNumberType> &independent_variables,
1792 const std::vector<ADNumberType> &dependent_variables,
1793 FullMatrix<ScalarType> & hessian) const
1794 {
1795 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
1796 2,
1797 ExcSupportedDerivativeLevels(
1798 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
1799 2));
1800 Assert(dependent_variables.size() == 1,
1801 ExcDimensionMismatch(dependent_variables.size(), 1));
1802 Assert(hessian.m() == independent_variables.size(),
1803 ExcDimensionMismatch(hessian.m(), independent_variables.size()));
1804 Assert(hessian.n() == independent_variables.size(),
1805 ExcDimensionMismatch(hessian.n(), independent_variables.size()));
1806
1807 // In reverse mode, the gradients are computed from the
1808 // independent variables (i.e. the adjoint)
1809 internal::reverse_mode_dependent_variable_activation(
1810 const_cast<ADNumberType &>(dependent_variables[0]));
1811 const std::size_t n_independent_variables = independent_variables.size();
1812 for (unsigned int i = 0; i < n_independent_variables; ++i)
1813 {
1814 using derivative_type =
1815 typename ADNumberTraits<ADNumberType>::derivative_type;
1816 const derivative_type gradient_i =
1817 ADNumberTraits<ADNumberType>::get_directional_derivative(
1818 independent_variables[i], i);
1819
1820 for (unsigned int j = 0; j <= i; ++j) // Symmetry
1821 {
1822 // Extract higher-order directional derivatives. Depending on
1823 // the AD number type, the result may be another AD number or a
1824 // floating point value.
1825 const ScalarType hessian_ij =
1826 internal::NumberType<ScalarType>::value(
1827 ADNumberTraits<derivative_type>::get_directional_derivative(
1828 gradient_i, j));
1829 hessian[i][j] = hessian_ij;
1830 if (i != j)
1831 hessian[j][i] = hessian_ij; // Symmetry
1832 }
1833 }
1834 }
1835
1836
1837 template <typename ADNumberType, typename ScalarType>
1838 void
1839 TapelessDrivers<ADNumberType,
1840 ScalarType,
1841 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1842 NumberTypes::sacado_rad ||
1843 ADNumberTraits<ADNumberType>::type_code ==
1844 NumberTypes::sacado_rad_dfad>>::
1845 values(const std::vector<ADNumberType> &dependent_variables,
1846 Vector<ScalarType> & values) const
1847 {
1848 Assert(values.size() == dependent_variables.size(),
1849 ExcDimensionMismatch(values.size(), dependent_variables.size()));
1850
1851 const std::size_t n_dependent_variables = dependent_variables.size();
1852 for (unsigned int i = 0; i < n_dependent_variables; ++i)
1853 values[i] = ADNumberTraits<ADNumberType>::get_scalar_value(
1854 dependent_variables[i]);
1855 }
1856
1857
1858 template <typename ADNumberType, typename ScalarType>
1859 void
1860 TapelessDrivers<ADNumberType,
1861 ScalarType,
1862 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1863 NumberTypes::sacado_rad ||
1864 ADNumberTraits<ADNumberType>::type_code ==
1865 NumberTypes::sacado_rad_dfad>>::
1866 jacobian(const std::vector<ADNumberType> &independent_variables,
1867 const std::vector<ADNumberType> &dependent_variables,
1868 FullMatrix<ScalarType> & jacobian) const
1869 {
1870 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
1871 1,
1872 ExcSupportedDerivativeLevels(
1873 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
1874 1));
1875 Assert(jacobian.m() == dependent_variables.size(),
1876 ExcDimensionMismatch(jacobian.m(), dependent_variables.size()));
1877 Assert(jacobian.n() == independent_variables.size(),
1878 ExcDimensionMismatch(jacobian.n(), independent_variables.size()));
1879
1880 const std::size_t n_independent_variables = independent_variables.size();
1881 const std::size_t n_dependent_variables = dependent_variables.size();
1882
1883 // In reverse mode, the gradients are computed from the
1884 // independent variables (i.e. the adjoint).
1885 // For a demonstration of why this accumulation process is
1886 // required, see the unit tests
1887 // sacado/basic_01b.cc and sacado/basic_02b.cc
1888 // Here we also take into consideration the derivative type:
1889 // The Sacado number may be of the nested variety, in which
1890 // case the effect of the accumulation process on the
1891 // sensitivities of the nested number need to be accounted for.
1892 using accumulation_type =
1893 typename ADNumberTraits<ADNumberType>::derivative_type;
1894 std::vector<accumulation_type> rad_accumulation(
1895 n_independent_variables,
1896 ::internal::NumberType<accumulation_type>::value(0.0));
1897
1898 for (unsigned int i = 0; i < n_dependent_variables; ++i)
1899 {
1900 internal::reverse_mode_dependent_variable_activation(
1901 const_cast<ADNumberType &>(dependent_variables[i]));
1902 for (unsigned int j = 0; j < n_independent_variables; ++j)
1903 {
1904 const accumulation_type df_i_dx_j =
1905 ADNumberTraits<ADNumberType>::get_directional_derivative(
1906 independent_variables[j],
1907 i /*This number doesn't really matter*/) -
1908 rad_accumulation[j];
1909 jacobian[i][j] =
1910 internal::NumberType<ScalarType>::value(df_i_dx_j);
1911 rad_accumulation[j] += df_i_dx_j;
1912 }
1913 }
1914 }
1915
1916
1917
1918 // Specialization for auto-differentiable numbers that use
1919 // forward mode to compute the first (and, if supported, second)
1920 // derivatives.
1921
1922# ifndef DOXYGEN
1923 template <typename ADNumberType, typename ScalarType>
1924 TapelessDrivers<
1925 ADNumberType,
1926 ScalarType,
1927 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1928 NumberTypes::adolc_tapeless ||
1929 ADNumberTraits<ADNumberType>::type_code ==
1930 NumberTypes::sacado_dfad ||
1931 ADNumberTraits<ADNumberType>::type_code ==
1932 NumberTypes::sacado_dfad_dfad>>::TapelessDrivers()
1933 : dependent_variable_marking_safe(false)
1934 {}
1935# endif
1936
1937
1938 template <typename ADNumberType, typename ScalarType>
1939 void
1940 TapelessDrivers<ADNumberType,
1941 ScalarType,
1942 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1943 NumberTypes::adolc_tapeless ||
1944 ADNumberTraits<ADNumberType>::type_code ==
1945 NumberTypes::sacado_dfad ||
1946 ADNumberTraits<ADNumberType>::type_code ==
1947 NumberTypes::sacado_dfad_dfad>>::
1948 initialize_global_environment(const unsigned int n_independent_variables)
1949 {
1950 internal::configure_tapeless_mode<ADNumberType>(n_independent_variables);
1951 }
1952
1953
1954 template <typename ADNumberType, typename ScalarType>
1955 void
1956 TapelessDrivers<
1957 ADNumberType,
1958 ScalarType,
1959 std::enable_if_t<
1960 ADNumberTraits<ADNumberType>::type_code ==
1961 NumberTypes::adolc_tapeless ||
1962 ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_dfad ||
1963 ADNumberTraits<ADNumberType>::type_code ==
1964 NumberTypes::sacado_dfad_dfad>>::allow_dependent_variable_marking()
1965 {
1966 dependent_variable_marking_safe = true;
1967 }
1968
1969
1970 template <typename ADNumberType, typename ScalarType>
1971 void
1972 TapelessDrivers<
1973 ADNumberType,
1974 ScalarType,
1975 std::enable_if_t<
1976 ADNumberTraits<ADNumberType>::type_code ==
1977 NumberTypes::adolc_tapeless ||
1978 ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_dfad ||
1979 ADNumberTraits<ADNumberType>::type_code ==
1980 NumberTypes::sacado_dfad_dfad>>::prevent_dependent_variable_marking()
1981 {
1982 dependent_variable_marking_safe = false;
1983 }
1984
1985
1986 template <typename ADNumberType, typename ScalarType>
1987 bool
1988 TapelessDrivers<ADNumberType,
1989 ScalarType,
1990 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
1991 NumberTypes::adolc_tapeless ||
1992 ADNumberTraits<ADNumberType>::type_code ==
1993 NumberTypes::sacado_dfad ||
1994 ADNumberTraits<ADNumberType>::type_code ==
1995 NumberTypes::sacado_dfad_dfad>>::
1996 is_dependent_variable_marking_allowed() const
1997 {
1998 return dependent_variable_marking_safe;
1999 }
2000
2001
2002 template <typename ADNumberType, typename ScalarType>
2003 ScalarType
2004 TapelessDrivers<ADNumberType,
2005 ScalarType,
2006 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
2007 NumberTypes::adolc_tapeless ||
2008 ADNumberTraits<ADNumberType>::type_code ==
2009 NumberTypes::sacado_dfad ||
2010 ADNumberTraits<ADNumberType>::type_code ==
2011 NumberTypes::sacado_dfad_dfad>>::
2012 value(const std::vector<ADNumberType> &dependent_variables) const
2013 {
2014 Assert(dependent_variables.size() == 1,
2015 ExcDimensionMismatch(dependent_variables.size(), 1));
2016 return ADNumberTraits<ADNumberType>::get_scalar_value(
2017 dependent_variables[0]);
2018 }
2019
2020
2021 template <typename ADNumberType, typename ScalarType>
2022 void
2023 TapelessDrivers<ADNumberType,
2024 ScalarType,
2025 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
2026 NumberTypes::adolc_tapeless ||
2027 ADNumberTraits<ADNumberType>::type_code ==
2028 NumberTypes::sacado_dfad ||
2029 ADNumberTraits<ADNumberType>::type_code ==
2030 NumberTypes::sacado_dfad_dfad>>::
2031 gradient(const std::vector<ADNumberType> &independent_variables,
2032 const std::vector<ADNumberType> &dependent_variables,
2033 Vector<ScalarType> & gradient) const
2034 {
2035 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
2036 1,
2037 ExcSupportedDerivativeLevels(
2038 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
2039 1));
2040 Assert(dependent_variables.size() == 1,
2041 ExcDimensionMismatch(dependent_variables.size(), 1));
2042 Assert(gradient.size() == independent_variables.size(),
2043 ExcDimensionMismatch(gradient.size(),
2044 independent_variables.size()));
2045
2046 // In forward mode, the gradients are computed from the
2047 // dependent variables
2048 const std::size_t n_independent_variables = independent_variables.size();
2049 for (unsigned int i = 0; i < n_independent_variables; ++i)
2050 gradient[i] = internal::NumberType<ScalarType>::value(
2051 ADNumberTraits<ADNumberType>::get_directional_derivative(
2052 dependent_variables[0], i));
2053 }
2054
2055
2056 template <typename ADNumberType, typename ScalarType>
2057 void
2058 TapelessDrivers<ADNumberType,
2059 ScalarType,
2060 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
2061 NumberTypes::adolc_tapeless ||
2062 ADNumberTraits<ADNumberType>::type_code ==
2063 NumberTypes::sacado_dfad ||
2064 ADNumberTraits<ADNumberType>::type_code ==
2065 NumberTypes::sacado_dfad_dfad>>::
2066 hessian(const std::vector<ADNumberType> &independent_variables,
2067 const std::vector<ADNumberType> &dependent_variables,
2068 FullMatrix<ScalarType> & hessian) const
2069 {
2070 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
2071 2,
2072 ExcSupportedDerivativeLevels(
2073 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
2074 2));
2075 Assert(dependent_variables.size() == 1,
2076 ExcDimensionMismatch(dependent_variables.size(), 1));
2077 Assert(hessian.m() == independent_variables.size(),
2078 ExcDimensionMismatch(hessian.m(), independent_variables.size()));
2079 Assert(hessian.n() == independent_variables.size(),
2080 ExcDimensionMismatch(hessian.n(), independent_variables.size()));
2081
2082 // In forward mode, the gradients are computed from the
2083 // dependent variables
2084 const std::size_t n_independent_variables = independent_variables.size();
2085 for (unsigned int i = 0; i < n_independent_variables; ++i)
2086 {
2087 using derivative_type =
2088 typename ADNumberTraits<ADNumberType>::derivative_type;
2089 const derivative_type gradient_i =
2090 ADNumberTraits<ADNumberType>::get_directional_derivative(
2091 dependent_variables[0], i);
2092
2093 for (unsigned int j = 0; j <= i; ++j) // Symmetry
2094 {
2095 // Extract higher-order directional derivatives. Depending on
2096 // the AD number type, the result may be another AD number or a
2097 // floating point value.
2098 const ScalarType hessian_ij =
2099 internal::NumberType<ScalarType>::value(
2100 ADNumberTraits<derivative_type>::get_directional_derivative(
2101 gradient_i, j));
2102 hessian[i][j] = hessian_ij;
2103 if (i != j)
2104 hessian[j][i] = hessian_ij; // Symmetry
2105 }
2106 }
2107 }
2108
2109
2110 template <typename ADNumberType, typename ScalarType>
2111 void
2112 TapelessDrivers<ADNumberType,
2113 ScalarType,
2114 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
2115 NumberTypes::adolc_tapeless ||
2116 ADNumberTraits<ADNumberType>::type_code ==
2117 NumberTypes::sacado_dfad ||
2118 ADNumberTraits<ADNumberType>::type_code ==
2119 NumberTypes::sacado_dfad_dfad>>::
2120 values(const std::vector<ADNumberType> &dependent_variables,
2121 Vector<ScalarType> & values) const
2122 {
2123 Assert(values.size() == dependent_variables.size(),
2124 ExcDimensionMismatch(values.size(), dependent_variables.size()));
2125
2126 const std::size_t n_dependent_variables = dependent_variables.size();
2127 for (unsigned int i = 0; i < n_dependent_variables; ++i)
2128 values[i] = ADNumberTraits<ADNumberType>::get_scalar_value(
2129 dependent_variables[i]);
2130 }
2131
2132
2133 template <typename ADNumberType, typename ScalarType>
2134 void
2135 TapelessDrivers<ADNumberType,
2136 ScalarType,
2137 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
2138 NumberTypes::adolc_tapeless ||
2139 ADNumberTraits<ADNumberType>::type_code ==
2140 NumberTypes::sacado_dfad ||
2141 ADNumberTraits<ADNumberType>::type_code ==
2142 NumberTypes::sacado_dfad_dfad>>::
2143 jacobian(const std::vector<ADNumberType> &independent_variables,
2144 const std::vector<ADNumberType> &dependent_variables,
2145 FullMatrix<ScalarType> & jacobian) const
2146 {
2147 Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
2148 1,
2149 ExcSupportedDerivativeLevels(
2150 AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
2151 1));
2152 Assert(jacobian.m() == dependent_variables.size(),
2153 ExcDimensionMismatch(jacobian.m(), dependent_variables.size()));
2154 Assert(jacobian.n() == independent_variables.size(),
2155 ExcDimensionMismatch(jacobian.n(), independent_variables.size()));
2156
2157 const std::size_t n_independent_variables = independent_variables.size();
2158 const std::size_t n_dependent_variables = dependent_variables.size();
2159
2160 // In forward mode, the gradients are computed from the
2161 // dependent variables
2162 for (unsigned int i = 0; i < n_dependent_variables; ++i)
2163 for (unsigned int j = 0; j < n_independent_variables; ++j)
2164 jacobian[i][j] = internal::NumberType<ScalarType>::value(
2165 ADNumberTraits<ADNumberType>::get_directional_derivative(
2166 dependent_variables[i], j));
2167 }
2168
2169
2170 } // namespace AD
2171} // namespace Differentiation
2172
2173
2174/* --- Explicit instantiations --- */
2175# include "ad_drivers.inst"
2176# ifdef DEAL_II_WITH_ADOLC
2177# include "ad_drivers.inst1"
2178# endif
2179# ifdef DEAL_II_TRILINOS_WITH_SACADO
2180# include "ad_drivers.inst2"
2181# endif
2182
2183
2184DEAL_II_NAMESPACE_CLOSE
2185
2186
2187#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
ad_drivers.h
ad_number_traits.h
ad_number_types.h
adolc_number_types.h
FullMatrix::n
size_type n() const
FullMatrix::m
size_type m() const
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
full_matrix.h
Differentiation::AD::ExcSupportedDerivativeLevels
static ::ExceptionBase & ExcSupportedDerivativeLevels(std::size_t arg1, std::size_t arg2)
Differentiation::AD::ExcRequiresADNumberSpecialization
static ::ExceptionBase & ExcRequiresADNumberSpecialization()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
Differentiation::AD::internal::reverse_mode_dependent_variable_activation
std::enable_if_t<!(ADNumberTraits< ADNumberType >::type_code==NumberTypes::sacado_rad||ADNumberTraits< ADNumberType >::type_code==NumberTypes::sacado_rad_dfad)> reverse_mode_dependent_variable_activation(ADNumberType &)
Definition ad_drivers.cc:1527
Differentiation::AD::internal::configure_tapeless_mode
std::enable_if_t<!(ADNumberTraits< ADNumberType >::type_code==NumberTypes::adolc_tapeless)> configure_tapeless_mode(const unsigned int)
Definition ad_drivers.cc:1569
sacado_number_types.h
Differentiation::AD::TapedDrivers::gradient
void gradient(const typename Types< ADNumberType >::tape_index active_tape_index, const std::vector< ScalarType > &independent_variables, Vector< ScalarType > &gradient) const
Definition ad_drivers.cc:208
Differentiation::AD::TapedDrivers::is_registered_tape
bool is_registered_tape(const typename Types< ADNumberType >::tape_index tape_index) const
Definition ad_drivers.cc:72
Differentiation::AD::TapedDrivers::stop_taping
void stop_taping(const typename Types< ADNumberType >::tape_index active_tape_index, const bool write_tapes_to_file)
Definition ad_drivers.cc:113
Differentiation::AD::TapedDrivers::remove_tape
void remove_tape(const typename Types< ADNumberType >::tape_index tape_index)
Definition ad_drivers.cc:162
Differentiation::AD::TapedDrivers::start_taping
void start_taping(const typename Types< ADNumberType >::tape_index tape_index, const bool keep_independent_values)
Definition ad_drivers.cc:103
Differentiation::AD::TapedDrivers::hessian
void hessian(const typename Types< ADNumberType >::tape_index active_tape_index, const std::vector< ScalarType > &independent_variables, FullMatrix< ScalarType > &hessian) const
Definition ad_drivers.cc:219
Differentiation::AD::TapedDrivers::set_tape_buffer_sizes
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:91
Differentiation::AD::TapedDrivers::activate_tape
void activate_tape(const typename Types< ADNumberType >::tape_index tape_index)
Definition ad_drivers.cc:133
Differentiation::AD::TapedDrivers::requires_retaping
bool requires_retaping(const typename Types< ADNumberType >::tape_index tape_index) const
Definition ad_drivers.cc:142
Differentiation::AD::TapedDrivers::print
void print(std::ostream &stream) const
Definition ad_drivers.cc:179
Differentiation::AD::TapedDrivers::print_tape_stats
void print_tape_stats(const typename Types< ADNumberType >::tape_index tape_index, std::ostream &stream) const
Definition ad_drivers.cc:187
Differentiation::AD::TapedDrivers::value
ScalarType value(const typename Types< ADNumberType >::tape_index active_tape_index, const std::vector< ScalarType > &independent_variables) const
Definition ad_drivers.cc:197
Differentiation::AD::TapedDrivers::active_tape_index
Types< ADNumberType >::tape_index active_tape_index() const
Definition ad_drivers.cc:63
Differentiation::AD::TapedDrivers::reset
void reset(const bool clear_registered_tapes)
Definition ad_drivers.cc:171
Differentiation::AD::TapedDrivers::get_registered_tape_indices
std::vector< typename Types< ADNumberType >::tape_index > get_registered_tape_indices() const
Definition ad_drivers.cc:123
Differentiation::AD::TapedDrivers::values
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
Definition ad_drivers.cc:230
Differentiation::AD::TapedDrivers::jacobian
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
Definition ad_drivers.cc:242
internal::NumberType::value
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition numbers.h:702
Generated by doxygen 1.9.7