Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
ad_drivers.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #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 
35 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
36 # include <adolc/adolc_fatalerror.h>
37 # include <adolc/drivers/drivers.h>
38 # include <adolc/taping.h>
39 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
40 
41 # endif // DEAL_II_WITH_ADOLC
42 
43 # include <vector>
44 
45 
46 DEAL_II_NAMESPACE_OPEN
47 
48 
49 namespace 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 
2263 DEAL_II_NAMESPACE_CLOSE
2264 
2265 
2266 #endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
Differentiation::AD::TapedDrivers::last_action_requires_retaping
bool last_action_requires_retaping() const
Definition: ad_drivers.cc:156
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::internal::configure_tapeless_mode
std::enable_if< ADNumberTraits< ADNumberType >::type_code==NumberTypes::adolc_tapeless >::type configure_tapeless_mode(const unsigned int n_directional_derivatives)
Definition: ad_drivers.cc:1647
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
types.h
internal::NumberType::value
static constexpr const T & value(const T &t)
Definition: numbers.h:703
Differentiation::AD::TapedDrivers::is_registered_tape
bool is_registered_tape(const typename Types< ADNumberType >::tape_index tape_index) const
Definition: ad_drivers.cc:76
utilities.h
FullMatrix::m
size_type m() const
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::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
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
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::Types< ADNumberType >::tape_index
unsigned int tape_index
Definition: ad_drivers.h:99
adolc_number_types.h
FullMatrix::n
size_type n() const
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 &dependent_variable)
Definition: ad_drivers.cc:1606
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::NumberTypes::sacado_rad_dfad
@ sacado_rad_dfad
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::NumberTypes::sacado_rad
@ sacado_rad
Differentiation::AD::ExcRequiresADNumberSpecialization
static ::ExceptionBase & ExcRequiresADNumberSpecialization()
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:409
Differentiation::AD::ExcRequiresADOLC
static ::ExceptionBase & ExcRequiresADOLC()
ad_number_types.h
Differentiation::AD::TapedDrivers::keep_independent_values
bool keep_independent_values() const
Definition: ad_drivers.cc:86
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
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::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
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
ad_drivers.h
Differentiation
Definition: numbers.h:645
Differentiation::AD::NumberTypes::sacado_dfad_dfad
@ sacado_dfad_dfad
Differentiation::AD::TapedDrivers::is_recording
bool is_recording() const
Definition: ad_drivers.cc:58
Differentiation::AD::ExcSupportedDerivativeLevels
static ::ExceptionBase & ExcSupportedDerivativeLevels(std::size_t arg1, std::size_t arg2)
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::activate_tape
void activate_tape(const typename Types< ADNumberType >::tape_index tape_index)
Definition: ad_drivers.cc:137
scalar_type
Differentiation::AD::NumberTypes::adolc_taped
@ adolc_taped
exceptions.h
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
value
static const bool value
Definition: dof_tools_constraints.cc:433
Differentiation::AD::Types
Definition: ad_drivers.h:94
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Differentiation::AD::TapedDrivers::print
void print(std::ostream &stream) const
Definition: ad_drivers.cc:183
sacado_number_types.h
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
ad_number_traits.h
internal::NumberType
Definition: numbers.h:700
vector.h
Differentiation::AD::ADNumberTraits
Definition: ad_number_traits.h:72
Differentiation::AD::NumberTypes::adolc_tapeless
@ adolc_tapeless
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
config.h
Differentiation::AD::TapelessDrivers
Definition: ad_drivers.h:522
FullMatrix
Definition: full_matrix.h:71
internal
Definition: aligned_vector.h:369
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Differentiation::AD::NumberTypes::sacado_dfad
@ sacado_dfad
Vector< double >
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
Definition: ad_drivers.h:150
Differentiation::AD::Numbers
Definition: ad_drivers.h:113
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
Differentiation::AD::Types< ADNumberType >::tape_buffer_sizes
unsigned int tape_buffer_sizes
Definition: ad_drivers.h:104
full_matrix.h
Differentiation::AD::TapedDrivers::remove_tape
void remove_tape(const typename Types< ADNumberType >::tape_index tape_index)
Definition: ad_drivers.cc:166
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:372
Generated by   doxygen 1.8.17