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