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