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