]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactor the AD drivers classes.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 6 Dec 2018 13:37:14 +0000 (14:37 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 24 Jan 2019 07:44:11 +0000 (08:44 +0100)
This commit significantly reworks the implementation of the AD driver
classes. Firstly, all class methods are moved to the associated source
file. The (logically) static class methods have been removed, so all
methods must be called through an instance of a driver class. The
previously implemented logic for existant ADOL-C tapes has been
superceeded by low-level calls to the library itself.

The ADHelper classes are updated to accommodate the changes to the
driver classes.

Fixes #7335
Fixes #7372

include/deal.II/differentiation/ad/ad_drivers.h
include/deal.II/differentiation/ad/ad_helpers.h
source/differentiation/ad/CMakeLists.txt
source/differentiation/ad/ad_drivers.cc [new file with mode: 0644]
source/differentiation/ad/ad_drivers.inst1.in [new file with mode: 0644]
source/differentiation/ad/ad_drivers.inst2.in [new file with mode: 0644]
source/differentiation/ad/ad_helpers.cc

index 05ec5584067338cfb67fa2f60027b1f87dd9ce7d..5af046665bf6799985b4b8f89a2e12046e9cb919 100644 (file)
 
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/types.h>
-#include <deal.II/base/utilities.h>
 
 #include <deal.II/differentiation/ad/ad_number_traits.h>
 #include <deal.II/differentiation/ad/ad_number_types.h>
-#include <deal.II/differentiation/ad/adolc_number_types.h>
-#include <deal.II/differentiation/ad/sacado_number_types.h>
 
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/vector.h>
 #ifdef DEAL_II_WITH_ADOLC
 
 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
-#  include <adolc/adolc_fatalerror.h>
-#  include <adolc/drivers/drivers.h>
 #  include <adolc/internal/usrparms.h>
-#  include <adolc/taping.h>
 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 
 #endif // DEAL_II_WITH_ADOLC
@@ -136,7 +130,7 @@ namespace Differentiation
 
 
     /**
-     * A driver class for taped auto-differentiable numbers.
+     * A prototype driver class for taped auto-differentiable numbers.
      *
      * It is intended that this class be specialized for the valid
      * combinations of auto-differentiable numbers and output scalar
@@ -258,6 +252,12 @@ namespace Differentiation
         const typename Types<ADNumberType>::tape_index active_tape_index,
         const bool                                     write_tapes_to_file);
 
+      /**
+       * Return a list of registered tape indices.
+       */
+      std::vector<typename Types<ADNumberType>::tape_index>
+      get_registered_tape_indices() const;
+
       /**
        * Select a tape to record to or read from.
        *
@@ -275,6 +275,39 @@ namespace Differentiation
       void
       activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
+      /**
+       * Return a flag that, when <code>true</code>, indicates that the retaping
+       * of the dependent function for the chosen @p tape_index is necessary for
+       * a reliable computation to be performed.
+       * This may be necessary a sign comparison within branched operations
+       * yields different results to those computed at the original tape
+       * evaluation point.
+       *
+       * @note The chosen tape index must be greater than
+       * Numbers<ADNumberType>::invalid_tape_index and less than
+       * Numbers<ADNumberType>::max_tape_index.
+       */
+      bool
+      requires_retaping(
+        const typename Types<ADNumberType>::tape_index tape_index) const;
+
+      /**
+       * Return a flag that, when <code>true</code>, indicates that the retaping
+       * of the dependent function is necessary for a reliable computation to be
+       * performed on the active tape.
+       * This may be necessary a sign comparison within branched operations
+       * yields different results to those computed at the original tape
+       * evaluation point.
+       */
+      bool
+      last_action_requires_retaping() const;
+
+      /**
+       * Completely erases the tape with the given @p tape_index.
+       */
+      void
+      remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
+
       /**
        * Reset the state of the class.
        *
@@ -302,10 +335,10 @@ namespace Differentiation
        * @param[out] stream The output stream to which the values are to be
        *            written.
        */
-      static void
+      void
       print_tape_stats(
         const typename Types<ADNumberType>::tape_index tape_index,
-        std::ostream &                                 stream);
+        std::ostream &                                 stream) const;
 
       //@}
 
@@ -324,9 +357,9 @@ namespace Differentiation
        *
        * @return The scalar value of the function.
        */
-      static ScalarType
+      ScalarType
       value(const typename Types<ADNumberType>::tape_index active_tape_index,
-            const std::vector<ScalarType> &independent_variables);
+            const std::vector<ScalarType> &independent_variables) const;
 
       /**
        * Compute the gradient of the scalar field with respect to all
@@ -341,10 +374,10 @@ namespace Differentiation
        *             correct size (with length
        *             <code>n_independent_variables</code>).
        */
-      static void
+      void
       gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
                const std::vector<ScalarType> &independent_variables,
-               Vector<ScalarType> &           gradient);
+               Vector<ScalarType> &           gradient) const;
 
       /**
        * Compute the Hessian of the scalar field with respect to all
@@ -359,10 +392,10 @@ namespace Differentiation
        *             size (with dimensions
        *             <code>n_independent_variables</code>$\times$<code>n_independent_variables</code>).
        */
-      static void
+      void
       hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
               const std::vector<ScalarType> &independent_variables,
-              FullMatrix<ScalarType> &       hessian);
+              FullMatrix<ScalarType> &       hessian) const;
 
       //@}
 
@@ -383,11 +416,11 @@ namespace Differentiation
        *             It is expected that this vector be of the correct size
        *             (with length <code>n_dependent_variables</code>).
        */
-      static void
+      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);
+             Vector<ScalarType> &           values) const;
 
       /**
        * Compute the Jacobian of the vector field.
@@ -407,11 +440,11 @@ namespace Differentiation
        *             size (with dimensions
        *             <code>n_dependent_variables</code>$\times$<code>n_independent_variables</code>).
        */
-      static void
+      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);
+               FullMatrix<ScalarType> &       jacobian) const;
 
       //@}
     };
@@ -510,8 +543,8 @@ namespace Differentiation
        *
        * @return The scalar value of the function.
        */
-      static ScalarType
-      value(const std::vector<ADNumberType> &dependent_variables);
+      ScalarType
+      value(const std::vector<ADNumberType> &dependent_variables) const;
 
       /**
        * Compute the gradient of the scalar field with respect to all
@@ -526,10 +559,10 @@ namespace Differentiation
        *             correct size (with length
        *             <code>n_independent_variables</code>).
        */
-      static void
+      void
       gradient(const std::vector<ADNumberType> &independent_variables,
                const std::vector<ADNumberType> &dependent_variables,
-               Vector<ScalarType> &             gradient);
+               Vector<ScalarType> &             gradient) const;
 
       /**
        * Compute the Hessian of the scalar field with respect to all
@@ -544,10 +577,10 @@ namespace Differentiation
        *             size (with dimensions
        *             <code>n_independent_variables</code>$\times$<code>n_independent_variables</code>).
        */
-      static void
+      void
       hessian(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
-              FullMatrix<ScalarType> &         hessian);
+              FullMatrix<ScalarType> &         hessian) const;
 
       //@}
 
@@ -565,9 +598,9 @@ namespace Differentiation
        *             It is expected that this vector be of the correct size
        *             (with length <code>n_dependent_variables</code>).
        */
-      static void
+      void
       values(const std::vector<ADNumberType> &dependent_variables,
-             Vector<ScalarType> &             values);
+             Vector<ScalarType> &             values) const;
 
       /**
        * Compute the Jacobian of the vector field.
@@ -586,10 +619,10 @@ namespace Differentiation
        *             size (with dimensions
        *             <code>n_dependent_variables</code>$\times$<code>n_independent_variables</code>).
        */
-      static void
+      void
       jacobian(const std::vector<ADNumberType> &independent_variables,
                const std::vector<ADNumberType> &dependent_variables,
-               FullMatrix<ScalarType> &         jacobian);
+               FullMatrix<ScalarType> &         jacobian) const;
 
       //@}
     };
@@ -608,170 +641,6 @@ namespace Differentiation
 {
   namespace AD
   {
-    // -------------   TapedDrivers   -------------
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    bool
-    TapedDrivers<ADNumberType, ScalarType, T>::is_recording() const
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-      return false;
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    typename Types<ADNumberType>::tape_index
-    TapedDrivers<ADNumberType, ScalarType, T>::active_tape_index() const
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-      return Numbers<ADNumberType>::invalid_tape_index;
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    bool
-    TapedDrivers<ADNumberType, ScalarType, T>::is_registered_tape(
-      const typename Types<ADNumberType>::tape_index) const
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-      return false;
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    bool
-    TapedDrivers<ADNumberType, ScalarType, T>::keep_independent_values() const
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-      return false;
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::set_tape_buffer_sizes(
-      const typename Types<ADNumberType>::tape_buffer_sizes,
-      const typename Types<ADNumberType>::tape_buffer_sizes,
-      const typename Types<ADNumberType>::tape_buffer_sizes,
-      const typename Types<ADNumberType>::tape_buffer_sizes)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::start_taping(
-      const typename Types<ADNumberType>::tape_index,
-      const bool)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::stop_taping(
-      const typename Types<ADNumberType>::tape_index,
-      const bool)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::activate_tape(
-      const typename Types<ADNumberType>::tape_index)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::reset(const bool)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::print(std::ostream &) const
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::print_tape_stats(
-      const typename Types<ADNumberType>::tape_index,
-      std::ostream &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    ScalarType
-    TapedDrivers<ADNumberType, ScalarType, T>::value(
-      const typename Types<ADNumberType>::tape_index,
-      const std::vector<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-      return ScalarType(0.0);
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::gradient(
-      const typename Types<ADNumberType>::tape_index,
-      const std::vector<ScalarType> &,
-      Vector<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::hessian(
-      const typename Types<ADNumberType>::tape_index,
-      const std::vector<ScalarType> &,
-      FullMatrix<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::values(
-      const typename Types<ADNumberType>::tape_index,
-      const unsigned int,
-      const std::vector<ScalarType> &,
-      Vector<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapedDrivers<ADNumberType, ScalarType, T>::jacobian(
-      const typename Types<ADNumberType>::tape_index,
-      const unsigned int,
-      const std::vector<ScalarType> &,
-      FullMatrix<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
 #  ifdef DEAL_II_WITH_ADOLC
 
     /**
@@ -846,309 +715,114 @@ namespace Differentiation
     {
       using scalar_type = double;
 
-      TapedDrivers()
-        : active_tape(Numbers<ADNumberType>::invalid_tape_index)
-        , keep_values(true)
-        , is_recording_flag(false)
-        , use_stored_taped_buffer_sizes(false)
-        , obufsize(0u)
-        , lbufsize(0u)
-        , vbufsize(0u)
-        , tbufsize(0u)
-      {}
+      /**
+       * Constructor
+       */
+      TapedDrivers();
+
 
-      // === Taping ===
+      /**
+       * @name Taping
+       */
+      //@{
 
       bool
-      is_recording() const
-      {
-        return is_recording_flag;
-      }
+      is_recording() const;
 
       typename Types<ADNumberType>::tape_index
-      active_tape_index() const
-      {
-        return active_tape;
-      }
+      active_tape_index() const;
 
       bool
-      keep_independent_values() const
-      {
-        return keep_values;
-      }
+      keep_independent_values() const;
 
       bool
       is_registered_tape(
-        const typename Types<ADNumberType>::tape_index tape_index) const
-      {
-        // See https://gitlab.com/adol-c/adol-c/issues/11#note_108341333
-        try
-          {
-            std::vector<std::size_t> counts(STAT_SIZE);
-            ::tapestats(tape_index, counts.data());
-            return true;
-          }
-        catch (const ::FatalError &exc)
-          {
-            return false;
-          }
-      }
+        const typename Types<ADNumberType>::tape_index tape_index) const;
 
       void
       set_tape_buffer_sizes(
         const typename Types<ADNumberType>::tape_buffer_sizes in_obufsize,
         const typename Types<ADNumberType>::tape_buffer_sizes in_lbufsize,
         const typename Types<ADNumberType>::tape_buffer_sizes in_vbufsize,
-        const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize)
-      {
-        // When valid for the chosen AD number type, these values will be used
-        // the next time start_recording_operations() is called.
-        obufsize                      = in_obufsize;
-        lbufsize                      = in_lbufsize;
-        vbufsize                      = in_vbufsize;
-        tbufsize                      = in_tbufsize;
-        use_stored_taped_buffer_sizes = true;
-      }
+        const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize);
 
       void
       start_taping(const typename Types<ADNumberType>::tape_index tape_index,
-                   const bool keep_independent_values)
-      {
-        if (use_stored_taped_buffer_sizes)
-          trace_on(tape_index,
-                   keep_independent_values,
-                   obufsize,
-                   lbufsize,
-                   vbufsize,
-                   tbufsize);
-        else
-          trace_on(tape_index, keep_independent_values);
-
-        // Set some other flags to their indicated / required values
-        keep_values       = keep_independent_values;
-        is_recording_flag = true;
-      }
+                   const bool keep_independent_values);
 
       void
       stop_taping(
         const typename Types<ADNumberType>::tape_index active_tape_index,
-        const bool                                     write_tapes_to_file)
-      {
-        if (write_tapes_to_file)
-          trace_off(active_tape_index); // Slow
-        else
-          trace_off(); // Fast(er)
-
-        // Now that we've turned tracing off, we've definitely
-        // stopped all tape recording.
-        is_recording_flag = false;
-
-        // If the keep_values flag is set, then we expect the user to use this
-        // tape immediately after recording it. There is therefore no need to
-        // invalidate it. However, there is now also no way to double-check
-        // that the newly recorded tape is indeed the active tape.
-        if (keep_independent_values() == false)
-          active_tape = Numbers<ADNumberType>::invalid_tape_index;
-      }
-
-      void
-      activate_tape(const typename Types<ADNumberType>::tape_index tape_index)
-      {
-        active_tape = tape_index;
-      }
-
-      void
-      reset(const bool clear_registered_tapes)
-      {
-        active_tape       = Numbers<ADNumberType>::invalid_tape_index;
-        is_recording_flag = false;
-        if (clear_registered_tapes)
-          {
-            const std::vector<typename Types<ADNumberType>::tape_index>
-              registered_tape_indices = get_registered_tape_indices();
-            for (const auto &tape_index : registered_tape_indices)
-              removeTape(tape_index, TapeRemovalType::ADOLC_REMOVE_COMPLETELY);
-          }
-      }
-
-      void
-      print(std::ostream &stream) const
-      {
-        const std::vector<typename Types<ADNumberType>::tape_index>
-          registered_tape_indices = get_registered_tape_indices();
-        stream << "Registered tapes: ";
-        auto it_registered_tape = registered_tape_indices.begin();
-        for (unsigned int i = 0; i < registered_tape_indices.size();
-             ++i, ++it_registered_tape)
-          stream << *it_registered_tape
-                 << (i < (registered_tape_indices.size() - 1) ? "," : "");
-        stream << "\n";
-
-        stream << "Keep values? " << keep_independent_values() << "\n";
-        stream << "Use stored tape buffer sizes? "
-               << use_stored_taped_buffer_sizes << "\n";
-      }
+        const bool                                     write_tapes_to_file);
 
-      static void
+      std::vector<typename Types<ADNumberType>::tape_index>
+      get_registered_tape_indices() const;
+
+      void
+      activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
+
+      bool
+      requires_retaping(
+        const typename Types<ADNumberType>::tape_index tape_index) const;
+
+      bool
+      last_action_requires_retaping() const;
+
+      void
+      remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
+
+      void
+      reset(const bool clear_registered_tapes);
+
+      void
+      print(std::ostream &stream) const;
+
+      void
       print_tape_stats(
         const typename Types<ADNumberType>::tape_index tape_index,
-        std::ostream &                                 stream)
-      {
-        // See ADOL-C manual section 2.1
-        // and adolc/taping.h
-        std::vector<std::size_t> counts(STAT_SIZE);
-        ::tapestats(tape_index, counts.data());
-        Assert(counts.size() >= 18, ExcInternalError());
-        stream
-          << "Tape index: " << tape_index << "\n"
-          << "Number of independent variables: " << counts[0] << "\n"
-          << "Number of dependent variables:   " << counts[1] << "\n"
-          << "Max number of live, active variables: " << counts[2] << "\n"
-          << "Size of taylor stack (number of overwrites): " << counts[3]
-          << "\n"
-          << "Operations buffer size: " << counts[4] << "\n"
-          << "Total number of recorded operations: " << counts[5] << "\n"
-          << "Operations file written or not: " << counts[6] << "\n"
-          << "Overall number of locations: " << counts[7] << "\n"
-          << "Locations file written or not: " << counts[8] << "\n"
-          << "Overall number of values: " << counts[9] << "\n"
-          << "Values file written or not: " << counts[10] << "\n"
-          << "Locations buffer size: " << counts[11] << "\n"
-          << "Values buffer size: " << counts[12] << "\n"
-          << "Taylor buffer size: " << counts[13] << "\n"
-          << "Number of eq_*_prod for sparsity pattern: " << counts[14] << "\n"
-          << "Use of 'min_op', deferred to 'abs_op' for piecewise calculations: "
-          << counts[15] << "\n"
-          << "Number of 'abs' calls that can switch branch: " << counts[16]
-          << "\n"
-          << "Number of parameters (doubles) interchangeable without retaping: "
-          << counts[17] << "\n"
-          << std::flush;
-      }
-
-
-      // === Scalar drivers ===
-
-      static scalar_type
-      value(const typename Types<ADNumberType>::tape_index active_tape_index,
-            const std::vector<scalar_type> &independent_variables)
-      {
-        scalar_type value = 0.0;
+        std::ostream &                                 stream) const;
+
+      //@}
 
-        ::function(active_tape_index,
-                   1, // Only one dependent variable
-                   independent_variables.size(),
-                   const_cast<double *>(independent_variables.data()),
-                   &value);
+      /**
+       * @name Drivers for scalar functions (one dependent variable)
+       */
+      //@{
 
-        return value;
-      }
+      scalar_type
+      value(const typename Types<ADNumberType>::tape_index active_tape_index,
+            const std::vector<scalar_type> &independent_variables) const;
 
-      static void
+      void
       gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
                const std::vector<scalar_type> &independent_variables,
-               Vector<scalar_type> &           gradient)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 1,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            1));
-        Assert(gradient.size() == independent_variables.size(),
-               ExcDimensionMismatch(gradient.size(),
-                                    independent_variables.size()));
-
-        // Note: ADOL-C's ::gradient function expects a *double as the last
-        // parameter. Here we take advantage of the fact that the data in the
-        // Vector class is aligned (e.g. stored as an Array)
-        ::gradient(active_tape_index,
-                   independent_variables.size(),
-                   const_cast<scalar_type *>(independent_variables.data()),
-                   gradient.data());
-      }
+               Vector<scalar_type> &           gradient) const;
 
-      static void
+      void
       hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
               const std::vector<scalar_type> &independent_variables,
-              FullMatrix<scalar_type> &       hessian)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            2));
-        Assert(hessian.m() == independent_variables.size(),
-               ExcDimensionMismatch(hessian.m(), independent_variables.size()));
-        Assert(hessian.n() == independent_variables.size(),
-               ExcDimensionMismatch(hessian.n(), independent_variables.size()));
-
-        const unsigned int n_independent_variables =
-          independent_variables.size();
-        std::vector<scalar_type *> H(n_independent_variables);
-        for (unsigned int i = 0; i < n_independent_variables; ++i)
-          H[i] = &hessian[i][0];
-
-        ::hessian(active_tape_index,
-                  n_independent_variables,
-                  const_cast<scalar_type *>(independent_variables.data()),
-                  H.data());
-
-        // ADOL-C builds only the lower-triangular part of the
-        // symmetric Hessian, so we should copy the relevant
-        // entries into the upper triangular part.
-        for (unsigned int i = 0; i < n_independent_variables; i++)
-          for (unsigned int j = 0; j < i; j++)
-            hessian[j][i] = hessian[i][j]; // Symmetry
-      }
-
-      // === Vector drivers ===
+              FullMatrix<scalar_type> &       hessian) const;
 
-      static void
+      //@}
+
+      /**
+       * @name Drivers for vector functions (multiple dependent variables)
+       */
+      //@{
+
+      void
       values(const typename Types<ADNumberType>::tape_index active_tape_index,
              const unsigned int              n_dependent_variables,
              const std::vector<scalar_type> &independent_variables,
-             Vector<scalar_type> &           values)
-      {
-        Assert(values.size() == n_dependent_variables,
-               ExcDimensionMismatch(values.size(), n_dependent_variables));
-
-        // Note: ADOL-C's ::function function expects a *double as the last
-        // parameter. Here we take advantage of the fact that the data in the
-        // Vector class is aligned (e.g. stored as an Array)
-        ::function(active_tape_index,
-                   n_dependent_variables,
-                   independent_variables.size(),
-                   const_cast<scalar_type *>(independent_variables.data()),
-                   values.data());
-      }
+             Vector<scalar_type> &           values) const;
 
-      static void
+      void
       jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
                const unsigned int              n_dependent_variables,
                const std::vector<scalar_type> &independent_variables,
-               FullMatrix<scalar_type> &       jacobian)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 1,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            1));
-        Assert(jacobian.m() == n_dependent_variables,
-               ExcDimensionMismatch(jacobian.m(), n_dependent_variables));
-        Assert(jacobian.n() == independent_variables.size(),
-               ExcDimensionMismatch(jacobian.n(),
-                                    independent_variables.size()));
-
-        std::vector<scalar_type *> J(n_dependent_variables);
-        for (unsigned int i = 0; i < n_dependent_variables; ++i)
-          J[i] = &jacobian[i][0];
-
-        ::jacobian(active_tape_index,
-                   n_dependent_variables,
-                   independent_variables.size(),
-                   independent_variables.data(),
-                   J.data());
-      }
+               FullMatrix<scalar_type> &       jacobian) const;
+
+      //@}
 
     protected:
       /**
@@ -1171,6 +845,33 @@ namespace Differentiation
        */
       bool is_recording_flag;
 
+      /**
+       * The status of the last function or derivative evaluation performed on
+       * the selected tape. As quoted from the ADOL-C manual, this can take
+       * on one of six different values with the following interpretation:
+       *
+       * - <b>+3:</b> The function is locally analytic.
+       * - <b>+2:</b> The function is locally analytic but the sparsity
+       *     structure (compared to the situation at the taping point) may have
+       *     changed, e.g. while at taping arguments <code>fmax(a,b)</code>
+       *     returned <code>a</code>, we get <code>b</code> at the argument
+       *     currently used.
+       * - <b>+1:</b> At least one of the functions <code>fmin</code>,
+       *     <code>fmax</code> or <code>fabs</code> is evaluated at a tie or
+       *     zero, respectively. Hence, the function to be differentiated is
+       *     Lipschitz-continuous but possibly non-differentiable.
+       * - <b>0:</b>  Some arithmetic comparison involving adoubles yields a
+       *     tie. Hence, the function to be differentiated may be discontinuous.
+       * - <b>-1:</b> An <code>adouble</code> comparison yields different
+       *     results from the evaluation point at which the tape was generated.
+       * - <b>-2:</b> The argument of a user-defined quadrature has changed from
+       *     the evaluation point at which the tape was generated.
+       *
+       * When the @p status variable takes a negative value, retaping of the dependent
+       * function is necessary for a reliable computation to be performed.
+       */
+      mutable std::map<typename Types<ADNumberType>::tape_index, int> status;
+
       /**
        * A flag indicating that we should preferentially use the user-defined
        * taped buffer sizes as opposed to either the default values selected
@@ -1198,30 +899,9 @@ namespace Differentiation
        * ADOL-C Taylor buffer size.
        */
       typename Types<ADNumberType>::tape_buffer_sizes tbufsize;
-
-      /**
-       * Return a list of registered tape indices
-       */
-      std::vector<typename Types<ADNumberType>::tape_index>
-      get_registered_tape_indices() const
-      {
-        // We've chosen to use unsigned shorts for the tape
-        // index type (a safety precaution) so we need to
-        // perform a conversion betwwen ADOL-C's native tape
-        // index type and that chosen by us.
-        std::vector<short> registered_tape_indices_s;
-        cachedTraceTags(registered_tape_indices_s);
-
-        std::vector<typename Types<ADNumberType>::tape_index>
-          registered_tape_indices(registered_tape_indices_s.size());
-        std::copy(registered_tape_indices_s.begin(),
-                  registered_tape_indices_s.end(),
-                  registered_tape_indices.begin());
-        return registered_tape_indices;
-      }
     };
 
-#  else // DEAL_II_WITH_ADOLC
+#  else
 
     /**
      * Specialization for taped ADOL-C auto-differentiable numbers.
@@ -1229,7 +909,7 @@ namespace Differentiation
      * Although we could revert to the default definition for the
      * unspecialized TapedDrivers class, we add this specialization
      * to provide a more descriptive error message if any of its
-     * static member functions are called.
+     * member functions are called.
      */
     template <typename ADNumberType>
     struct TapedDrivers<
@@ -1240,134 +920,107 @@ namespace Differentiation
     {
       using scalar_type = double;
 
-      // === Taping ===
+      /**
+       * @name Taping
+       */
+      //@{
 
       bool
-      is_recording() const
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-        return false;
-      }
+      is_recording() const;
 
       typename Types<ADNumberType>::tape_index
-      active_tape_index() const
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-        return Numbers<ADNumberType>::invalid_tape_index;
-      }
+      active_tape_index() const;
 
       bool
-      keep_independent_values() const
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-        return false;
-      }
+      keep_independent_values() const;
 
       bool
-      is_registered_tape(const typename Types<ADNumberType>::tape_index) const
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-        return false;
-      }
+      is_registered_tape(
+        const typename Types<ADNumberType>::tape_index tape_index) const;
 
       void
       set_tape_buffer_sizes(
         const typename Types<ADNumberType>::tape_buffer_sizes,
         const typename Types<ADNumberType>::tape_buffer_sizes,
         const typename Types<ADNumberType>::tape_buffer_sizes,
-        const typename Types<ADNumberType>::tape_buffer_sizes)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+        const typename Types<ADNumberType>::tape_buffer_sizes);
 
       void
-      start_taping(const typename Types<ADNumberType>::tape_index, const bool)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+      start_taping(const typename Types<ADNumberType>::tape_index, const bool);
 
       void
-      stop_taping(const typename Types<ADNumberType>::tape_index, const bool)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+      stop_taping(const typename Types<ADNumberType>::tape_index, const bool);
+
+      std::vector<typename Types<ADNumberType>::tape_index>
+      get_registered_tape_indices() const;
 
       void
-      activate_tape(const typename Types<ADNumberType>::tape_index)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+      activate_tape(const typename Types<ADNumberType>::tape_index);
+
+      bool
+      requires_retaping(const typename Types<ADNumberType>::tape_index) const;
+
+      bool
+      last_action_requires_retaping() const;
 
       void
-      reset(const bool)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+      remove_tape(const typename Types<ADNumberType>::tape_index);
 
       void
-      print(std::ostream &) const
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+      reset(const bool);
 
-      static void
+      void
+      print(std::ostream &stream) const;
+
+      void
       print_tape_stats(const typename Types<ADNumberType>::tape_index,
-                       std::ostream &)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+                       std::ostream &) const;
 
+      //@}
 
-      // === Scalar drivers ===
+      /**
+       * @name Drivers for scalar functions (one dependent variable)
+       */
+      //@{
 
-      static scalar_type
+      scalar_type
       value(const typename Types<ADNumberType>::tape_index,
-            const std::vector<scalar_type> &)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-        return 0.0;
-      }
+            const std::vector<scalar_type> &) const;
 
-      static void
+      void
       gradient(const typename Types<ADNumberType>::tape_index,
                const std::vector<scalar_type> &,
-               Vector<scalar_type> &)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+               Vector<scalar_type> &) const;
 
-      static void
+      void
       hessian(const typename Types<ADNumberType>::tape_index,
               const std::vector<scalar_type> &,
-              FullMatrix<scalar_type> &)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+              FullMatrix<scalar_type> &) const;
 
-      // === Vector drivers ===
+      //@}
 
-      static void
+      /**
+       * @name Drivers for vector functions (multiple dependent variables)
+       */
+      //@{
+
+      void
       values(const typename Types<ADNumberType>::tape_index,
              const unsigned int,
              const std::vector<scalar_type> &,
-             Vector<scalar_type> &)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+             Vector<scalar_type> &) const;
 
-      static void
+      void
       jacobian(const typename Types<ADNumberType>::tape_index,
                const unsigned int,
                const std::vector<scalar_type> &,
-               FullMatrix<scalar_type> &)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
+               FullMatrix<scalar_type> &) const;
+
+      //@}
     };
 
 #  endif // DEAL_II_WITH_ADOLC
 
-
     /**
      * Specialization for ADOL-C taped numbers. It is expected that the
      * scalar return type for this class is a float.
@@ -1385,191 +1038,115 @@ namespace Differentiation
     {
       using scalar_type = float;
 
-      // === Taping ===
+      /**
+       * @name Taping
+       */
+      //@{
 
       bool
-      is_recording() const
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        return taped_driver.is_recording();
-      }
+      is_recording() const;
 
       typename Types<ADNumberType>::tape_index
-      active_tape_index() const
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        return taped_driver.active_tape_index();
-      }
+      active_tape_index() const;
 
       bool
-      keep_independent_values() const
-      {
-        return taped_driver.keep_independent_values();
-      }
+      keep_independent_values() const;
 
       bool
       is_registered_tape(
-        const typename Types<ADNumberType>::tape_index tape_index) const
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        return taped_driver.is_registered_tape(tape_index);
-      }
+        const typename Types<ADNumberType>::tape_index tape_index) const;
 
       void
       set_tape_buffer_sizes(
         const typename Types<ADNumberType>::tape_buffer_sizes obufsize,
         const typename Types<ADNumberType>::tape_buffer_sizes lbufsize,
         const typename Types<ADNumberType>::tape_buffer_sizes vbufsize,
-        const typename Types<ADNumberType>::tape_buffer_sizes tbufsize)
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        taped_driver.set_tape_buffer_sizes(obufsize,
-                                           lbufsize,
-                                           vbufsize,
-                                           tbufsize);
-      }
+        const typename Types<ADNumberType>::tape_buffer_sizes tbufsize);
 
       void
       start_taping(const typename Types<ADNumberType>::tape_index tape_index,
-                   const bool keep_independent_values)
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        taped_driver.start_taping(tape_index, keep_independent_values);
-      }
+                   const bool keep_independent_values);
 
       void
       stop_taping(
         const typename Types<ADNumberType>::tape_index active_tape_index,
-        const bool                                     write_tapes_to_file)
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        taped_driver.stop_taping(active_tape_index, write_tapes_to_file);
-      }
+        const bool                                     write_tapes_to_file);
+
+      std::vector<typename Types<ADNumberType>::tape_index>
+      get_registered_tape_indices() const;
+
+      void
+      activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
+
+      bool
+      requires_retaping(
+        const typename Types<ADNumberType>::tape_index tape_index) const;
+
+      bool
+      last_action_requires_retaping() const;
 
       void
-      activate_tape(const typename Types<ADNumberType>::tape_index tape_index)
-      {
-        taped_driver.activate_tape(tape_index);
-      }
+      remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
       void
-      reset(const bool clear_registered_tapes)
-      {
-        taped_driver.reset(clear_registered_tapes);
-      }
+      reset(const bool clear_registered_tapes);
 
       void
-      print(std::ostream &stream) const
-      {
-        taped_driver.print(stream);
-      }
+      print(std::ostream &stream) const;
 
-      static void
+      void
       print_tape_stats(
         const typename Types<ADNumberType>::tape_index tape_index,
-        std::ostream &                                 stream)
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        TapedDrivers<ADNumberType, double>::print_tape_stats(tape_index,
-                                                             stream);
-      }
+        std::ostream &                                 stream) const;
+
+      //@}
 
-      // === Scalar drivers ===
+      /**
+       * @name Drivers for scalar functions (one dependent variable)
+       */
+      //@{
 
-      static scalar_type
+      scalar_type
       value(const typename Types<ADNumberType>::tape_index active_tape_index,
-            const std::vector<scalar_type> &independent_variables)
-      {
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        return TapedDrivers<ADNumberType, double>::value(
-          active_tape_index, vector_float_to_double(independent_variables));
-      }
+            const std::vector<scalar_type> &independent_variables) const;
 
-      static void
+      void
       gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
                const std::vector<scalar_type> &independent_variables,
-               Vector<scalar_type> &           gradient)
-      {
-        Vector<double> gradient_double(gradient.size());
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        TapedDrivers<ADNumberType, double>::gradient(active_tape_index,
-                                                     vector_float_to_double(
-                                                       independent_variables),
-                                                     gradient_double);
-        gradient = gradient_double;
-      }
+               Vector<scalar_type> &           gradient) const;
 
-      static void
+      void
       hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
               const std::vector<scalar_type> &independent_variables,
-              FullMatrix<scalar_type> &       hessian)
-      {
-        FullMatrix<double> hessian_double(hessian.m(), hessian.n());
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        TapedDrivers<ADNumberType, double>::hessian(active_tape_index,
-                                                    vector_float_to_double(
-                                                      independent_variables),
-                                                    hessian_double);
-        hessian = hessian_double;
-      }
-
-      // === Vector drivers ===
+              FullMatrix<scalar_type> &       hessian) const;
 
-      static void
+      //@}
+
+      /**
+       * @name Drivers for vector functions (multiple dependent variables)
+       */
+      //@{
+
+      void
       values(const typename Types<ADNumberType>::tape_index active_tape_index,
              const unsigned int              n_dependent_variables,
              const std::vector<scalar_type> &independent_variables,
-             Vector<scalar_type> &           values)
-      {
-        Vector<double> values_double(values.size());
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        TapedDrivers<ADNumberType, double>::values(active_tape_index,
-                                                   n_dependent_variables,
-                                                   vector_float_to_double(
-                                                     independent_variables),
-                                                   values_double);
-        values = values_double;
-      }
+             Vector<scalar_type> &           values) const;
 
-      static void
+      void
       jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
                const unsigned int              n_dependent_variables,
                const std::vector<scalar_type> &independent_variables,
-               FullMatrix<scalar_type> &       jacobian)
-      {
-        FullMatrix<double> jacobian_double(jacobian.m(), jacobian.n());
-        // ADOL-C only supports 'double', not 'float', so we can forward to
-        // the 'double' implementation of this function
-        TapedDrivers<ADNumberType, double>::jacobian(active_tape_index,
-                                                     n_dependent_variables,
-                                                     vector_float_to_double(
-                                                       independent_variables),
-                                                     jacobian_double);
-        jacobian = jacobian_double;
-      }
+               FullMatrix<scalar_type> &       jacobian) const;
+
+      //@}
 
     private:
       /**
        * Copy a vector of floats into a vector of doubles
        */
-      static std::vector<double>
-      vector_float_to_double(const std::vector<float> &in)
-      {
-        std::vector<double> out(in.size());
-        std::copy(in.begin(), in.end(), out.begin());
-        return out;
-      }
+      std::vector<double>
+      vector_float_to_double(const std::vector<float> &in) const;
 
       /**
        * The object that actually takes care of the taping
@@ -1581,235 +1158,6 @@ namespace Differentiation
     // -------------   TapelessDrivers   -------------
 
 
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapelessDrivers<ADNumberType, ScalarType, T>::initialize_global_environment(
-      const unsigned int)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapelessDrivers<ADNumberType, ScalarType, T>::
-      allow_dependent_variable_marking()
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapelessDrivers<ADNumberType, ScalarType, T>::
-      prevent_dependent_variable_marking()
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    bool
-    TapelessDrivers<ADNumberType, ScalarType, T>::
-      is_dependent_variable_marking_allowed() const
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-      return false;
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    ScalarType
-    TapelessDrivers<ADNumberType, ScalarType, T>::value(
-      const std::vector<ADNumberType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-      return ScalarType(0.0);
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapelessDrivers<ADNumberType, ScalarType, T>::gradient(
-      const std::vector<ADNumberType> &,
-      const std::vector<ADNumberType> &,
-      Vector<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapelessDrivers<ADNumberType, ScalarType, T>::hessian(
-      const std::vector<ADNumberType> &,
-      const std::vector<ADNumberType> &,
-      FullMatrix<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapelessDrivers<ADNumberType, ScalarType, T>::values(
-      const std::vector<ADNumberType> &,
-      Vector<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    template <typename ADNumberType, typename ScalarType, typename T>
-    void
-    TapelessDrivers<ADNumberType, ScalarType, T>::jacobian(
-      const std::vector<ADNumberType> &,
-      const std::vector<ADNumberType> &,
-      FullMatrix<ScalarType> &)
-    {
-      AssertThrow(false, ExcRequiresADNumberSpecialization());
-    }
-
-
-    namespace internal
-    {
-      /**
-       * A dummy function to define the active dependent variable when using
-       * reverse-mode AD.
-       */
-      template <typename ADNumberType>
-      typename 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 &)
-      {}
-
-#  ifdef DEAL_II_TRILINOS_WITH_SACADO
-
-
-      /**
-       * Define the active dependent variable when using reverse-mode AD.
-       *
-       * If there are multiple dependent variables then it is necessary to
-       * inform the independent variables, from which the adjoints are computed,
-       * which dependent variable they are computing the gradients with respect
-       * to. This function broadcasts this information.
-       */
-      template <typename ADNumberType>
-      typename 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)
-      {
-        // Compute all gradients (adjoints) for this
-        // reverse-mode Sacado dependent variable.
-        // For reverse-mode Sacado numbers it is necessary to broadcast to
-        // all independent variables that it is time to compute gradients.
-        // For one dependent variable one would just need to call
-        // ADNumberType::Gradcomp(), but since we have a more
-        // generic implementation for vectors of dependent variables
-        // (vector mode) we default to the complex case.
-        ADNumberType::Outvar_Gradcomp(dependent_variable);
-      }
-
-#  endif
-
-
-      /**
-       * A dummy function to enable vector mode for tapeless
-       * auto-differentiable numbers.
-       */
-      template <typename ADNumberType>
-      typename std::enable_if<!(ADNumberTraits<ADNumberType>::type_code ==
-                                NumberTypes::adolc_tapeless)>::type
-      configure_tapeless_mode(const unsigned int)
-      {}
-
-
-#  ifdef DEAL_II_WITH_ADOLC
-
-
-      /**
-       * Enable vector mode for ADOL-C tapeless numbers.
-       *
-       * This function checks to see if its legal to increase the maximum
-       * number of directional derivatives to be considered during calculations.
-       * If not then it throws an error.
-       */
-      template <typename ADNumberType>
-      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
-                              NumberTypes::adolc_tapeless>::type
-      configure_tapeless_mode(const unsigned int n_directional_derivatives)
-      {
-#    ifdef DEAL_II_ADOLC_WITH_TAPELESS_REFCOUNTING
-        // See ADOL-C manual section 7.1
-        //
-        // NOTE: It is critical that this is done for tapeless mode BEFORE
-        // any adtl::adouble are created. If this is not done, then we see
-        // this scary warning:
-        //
-        // "
-        // ADOL-C Warning: Tapeless: Setting numDir could change memory
-        // allocation of derivatives in existing adoubles and may lead to
-        // erroneous results or memory corruption
-        // "
-        //
-        // So we use this dummy function to configure this setting before
-        // we create and initialize our class data
-        const std::size_t n_live_variables = adtl::refcounter::getNumLiveVar();
-        if (n_live_variables == 0)
-          {
-            adtl::setNumDir(n_directional_derivatives);
-          }
-        else
-          {
-            // So there are some live active variables floating around. Here we
-            // check if we ask to increase the number of computable
-            // directional derivatives. If this really is necessary then it's
-            // absolutely vital that there exist no live variables before doing
-            // so.
-            const std::size_t n_set_directional_derivatives = adtl::getNumDir();
-            if (n_directional_derivatives > n_set_directional_derivatives)
-              AssertThrow(
-                n_live_variables == 0,
-                ExcMessage(
-                  "There are currently " +
-                  Utilities::to_string(n_live_variables) +
-                  " live "
-                  "adtl::adouble variables in existence. They currently "
-                  "assume " +
-                  Utilities::to_string(n_set_directional_derivatives) +
-                  " directional derivatives "
-                  "but you wish to increase this to " +
-                  Utilities::to_string(n_directional_derivatives) +
-                  ". \n"
-                  "To safely change (or more specifically in this case, "
-                  "increase) the number of directional derivatives, there "
-                  "must be no tapeless doubles in local/global scope."));
-          }
-#    else
-        // If ADOL-C is not configured with tapeless number reference counting
-        // then there is no way that we can guarantee that the following call
-        // is safe. No comment... :-/
-        adtl::setNumDir(n_directional_derivatives);
-#    endif
-      }
-
-#  else // DEAL_II_WITH_ADOLC
-
-      template <typename ADNumberType>
-      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
-                              NumberTypes::adolc_tapeless>::type
-      configure_tapeless_mode(const unsigned int /*n_directional_derivatives*/)
-      {
-        AssertThrow(false, ExcRequiresADOLC());
-      }
-
-#  endif
-
-    } // namespace internal
-
-
     /**
      * Specialization for auto-differentiable numbers that use
      * reverse mode to compute the first derivatives (and, if supported,
@@ -1824,192 +1172,72 @@ namespace Differentiation
                               ADNumberTraits<ADNumberType>::type_code ==
                                 NumberTypes::sacado_rad_dfad>::type>
     {
-      TapelessDrivers()
-        : dependent_variable_marking_safe(false)
-      {}
+      /**
+       * Constructor
+       */
+      TapelessDrivers();
 
-      // === Configuration ===
+      /**
+       * @name Configuration
+       */
+      //@{
 
       static void
-      initialize_global_environment(const unsigned int n_independent_variables)
-      {
-        internal::configure_tapeless_mode<ADNumberType>(
-          n_independent_variables);
-      }
+      initialize_global_environment(const unsigned int n_independent_variables);
+
+      //@}
+
+      /**
+       * Operation status
+       */
+      //@{
 
-      // === Operation status ===
       void
-      allow_dependent_variable_marking()
-      {
-        dependent_variable_marking_safe = true;
-      }
+      allow_dependent_variable_marking();
 
       void
-      prevent_dependent_variable_marking()
-      {
-        dependent_variable_marking_safe = false;
-      }
+      prevent_dependent_variable_marking();
 
       bool
-      is_dependent_variable_marking_allowed() const
-      {
-        return dependent_variable_marking_safe;
-      }
-
-      // === Scalar drivers ===
-
-      static ScalarType
-      value(const std::vector<ADNumberType> &dependent_variables)
-      {
-        Assert(dependent_variables.size() == 1,
-               ExcDimensionMismatch(dependent_variables.size(), 1));
-        return ADNumberTraits<ADNumberType>::get_scalar_value(
-          dependent_variables[0]);
-      }
+      is_dependent_variable_marking_allowed() const;
 
-      static void
+      //@}
+
+      /**
+       * @name Drivers for scalar functions
+       */
+      //@{
+
+      ScalarType
+      value(const std::vector<ADNumberType> &dependent_variables) const;
+
+      void
       gradient(const std::vector<ADNumberType> &independent_variables,
                const std::vector<ADNumberType> &dependent_variables,
-               Vector<ScalarType> &             gradient)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 1,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            1));
-        Assert(dependent_variables.size() == 1,
-               ExcDimensionMismatch(dependent_variables.size(), 1));
-        Assert(gradient.size() == independent_variables.size(),
-               ExcDimensionMismatch(gradient.size(),
-                                    independent_variables.size()));
-
-        // In reverse mode, the gradients are computed from the
-        // independent variables (i.e. the adjoint)
-        internal::reverse_mode_dependent_variable_activation(
-          const_cast<ADNumberType &>(dependent_variables[0]));
-        const std::size_t n_independent_variables =
-          independent_variables.size();
-        for (unsigned int i = 0; i < n_independent_variables; i++)
-          gradient[i] = internal::NumberType<ScalarType>::value(
-            ADNumberTraits<ADNumberType>::get_directional_derivative(
-              independent_variables[i],
-              0 /*This number doesn't really matter*/));
-      }
+               Vector<ScalarType> &             gradient) const;
 
-      static void
+      void
       hessian(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
-              FullMatrix<ScalarType> &         hessian)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            2));
-        Assert(dependent_variables.size() == 1,
-               ExcDimensionMismatch(dependent_variables.size(), 1));
-        Assert(hessian.m() == independent_variables.size(),
-               ExcDimensionMismatch(hessian.m(), independent_variables.size()));
-        Assert(hessian.n() == independent_variables.size(),
-               ExcDimensionMismatch(hessian.n(), independent_variables.size()));
-
-        // In reverse mode, the gradients are computed from the
-        // independent variables (i.e. the adjoint)
-        internal::reverse_mode_dependent_variable_activation(
-          const_cast<ADNumberType &>(dependent_variables[0]));
-        const std::size_t n_independent_variables =
-          independent_variables.size();
-        for (unsigned int i = 0; i < n_independent_variables; i++)
-          {
-            using derivative_type =
-              typename ADNumberTraits<ADNumberType>::derivative_type;
-            const derivative_type gradient_i =
-              ADNumberTraits<ADNumberType>::get_directional_derivative(
-                independent_variables[i], i);
-
-            for (unsigned int j = 0; j <= i; ++j) // Symmetry
-              {
-                // Extract higher-order directional derivatives. Depending on
-                // the AD number type, the result may be another AD number or a
-                // floating point value.
-                const ScalarType hessian_ij =
-                  internal::NumberType<ScalarType>::value(
-                    ADNumberTraits<derivative_type>::get_directional_derivative(
-                      gradient_i, j));
-                hessian[i][j] = hessian_ij;
-                if (i != j)
-                  hessian[j][i] = hessian_ij; // Symmetry
-              }
-          }
-      }
-
-      // === Vector drivers ===
+              FullMatrix<ScalarType> &         hessian) const;
 
-      static void
-      values(const std::vector<ADNumberType> &dependent_variables,
-             Vector<ScalarType> &             values)
-      {
-        Assert(values.size() == dependent_variables.size(),
-               ExcDimensionMismatch(values.size(), dependent_variables.size()));
+      //@}
 
-        const std::size_t n_dependent_variables = dependent_variables.size();
-        for (unsigned int i = 0; i < n_dependent_variables; i++)
-          values[i] = ADNumberTraits<ADNumberType>::get_scalar_value(
-            dependent_variables[i]);
-      }
+      /**
+       * @name Drivers for vector functions
+       */
+      //@{
 
-      static void
+      void
+      values(const std::vector<ADNumberType> &dependent_variables,
+             Vector<ScalarType> &             values) const;
+
+      void
       jacobian(const std::vector<ADNumberType> &independent_variables,
                const std::vector<ADNumberType> &dependent_variables,
-               FullMatrix<ScalarType> &         jacobian)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 1,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            1));
-        Assert(jacobian.m() == dependent_variables.size(),
-               ExcDimensionMismatch(jacobian.m(), dependent_variables.size()));
-        Assert(jacobian.n() == independent_variables.size(),
-               ExcDimensionMismatch(jacobian.n(),
-                                    independent_variables.size()));
-
-        const std::size_t n_independent_variables =
-          independent_variables.size();
-        const std::size_t n_dependent_variables = dependent_variables.size();
-
-        // In reverse mode, the gradients are computed from the
-        // independent variables (i.e. the adjoint).
-        // For a demonstration of why this accumulation process is
-        // required, see the unit tests
-        // sacado/basic_01b.cc and sacado/basic_02b.cc
-        // Here we also take into consideration the derivative type:
-        // The Sacado number may be of the nested variety, in which
-        // case the effect of the accumulation process on the
-        // sensitivities of the nested number need to be accounted for.
-        using accumulation_type =
-          typename ADNumberTraits<ADNumberType>::derivative_type;
-        std::vector<accumulation_type> rad_accumulation(
-          n_independent_variables,
-          dealii::internal::NumberType<accumulation_type>::value(0.0));
-
-        for (unsigned int i = 0; i < n_dependent_variables; i++)
-          {
-            internal::reverse_mode_dependent_variable_activation(
-              const_cast<ADNumberType &>(dependent_variables[i]));
-            for (unsigned int j = 0; j < n_independent_variables; j++)
-              {
-                const accumulation_type df_i_dx_j =
-                  ADNumberTraits<ADNumberType>::get_directional_derivative(
-                    independent_variables[j],
-                    i /*This number doesn't really matter*/) -
-                  rad_accumulation[j];
-                jacobian[i][j] =
-                  internal::NumberType<ScalarType>::value(df_i_dx_j);
-                rad_accumulation[j] += df_i_dx_j;
-              }
-          }
-      }
+               FullMatrix<ScalarType> &         jacobian) const;
+
+      //@}
 
     private:
       /**
@@ -2036,163 +1264,72 @@ namespace Differentiation
                               ADNumberTraits<ADNumberType>::type_code ==
                                 NumberTypes::sacado_dfad_dfad>::type>
     {
-      TapelessDrivers()
-        : dependent_variable_marking_safe(false)
-      {}
+      /**
+       * Constructor
+       */
+      TapelessDrivers();
 
-      // === Configuration ===
+      /**
+       * @name Configuration
+       */
+      //@{
 
       static void
-      initialize_global_environment(const unsigned int n_independent_variables)
-      {
-        internal::configure_tapeless_mode<ADNumberType>(
-          n_independent_variables);
-      }
+      initialize_global_environment(const unsigned int n_independent_variables);
+
+      //@}
+
+      /**
+       * Operation status
+       */
+      //@{
 
-      // === Operation status ===
       void
-      allow_dependent_variable_marking()
-      {
-        dependent_variable_marking_safe = true;
-      }
+      allow_dependent_variable_marking();
 
       void
-      prevent_dependent_variable_marking()
-      {
-        dependent_variable_marking_safe = false;
-      }
+      prevent_dependent_variable_marking();
 
       bool
-      is_dependent_variable_marking_allowed() const
-      {
-        return dependent_variable_marking_safe;
-      }
-
-      // === Scalar drivers ===
-
-      static ScalarType
-      value(const std::vector<ADNumberType> &dependent_variables)
-      {
-        Assert(dependent_variables.size() == 1,
-               ExcDimensionMismatch(dependent_variables.size(), 1));
-        return ADNumberTraits<ADNumberType>::get_scalar_value(
-          dependent_variables[0]);
-      }
+      is_dependent_variable_marking_allowed() const;
 
-      static void
+      //@}
+
+      /**
+       * @name Drivers for scalar functions
+       */
+      //@{
+
+      ScalarType
+      value(const std::vector<ADNumberType> &dependent_variables) const;
+
+      void
       gradient(const std::vector<ADNumberType> &independent_variables,
                const std::vector<ADNumberType> &dependent_variables,
-               Vector<ScalarType> &             gradient)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 1,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            1));
-        Assert(dependent_variables.size() == 1,
-               ExcDimensionMismatch(dependent_variables.size(), 1));
-        Assert(gradient.size() == independent_variables.size(),
-               ExcDimensionMismatch(gradient.size(),
-                                    independent_variables.size()));
-
-        // In forward mode, the gradients are computed from the
-        // dependent variables
-        const std::size_t n_independent_variables =
-          independent_variables.size();
-        for (unsigned int i = 0; i < n_independent_variables; i++)
-          gradient[i] = internal::NumberType<ScalarType>::value(
-            ADNumberTraits<ADNumberType>::get_directional_derivative(
-              dependent_variables[0], i));
-      }
+               Vector<ScalarType> &             gradient) const;
 
-      static void
+      void
       hessian(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
-              FullMatrix<ScalarType> &         hessian)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 2,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            2));
-        Assert(dependent_variables.size() == 1,
-               ExcDimensionMismatch(dependent_variables.size(), 1));
-        Assert(hessian.m() == independent_variables.size(),
-               ExcDimensionMismatch(hessian.m(), independent_variables.size()));
-        Assert(hessian.n() == independent_variables.size(),
-               ExcDimensionMismatch(hessian.n(), independent_variables.size()));
-
-        // In forward mode, the gradients are computed from the
-        // dependent variables
-        const std::size_t n_independent_variables =
-          independent_variables.size();
-        for (unsigned int i = 0; i < n_independent_variables; i++)
-          {
-            using derivative_type =
-              typename ADNumberTraits<ADNumberType>::derivative_type;
-            const derivative_type gradient_i =
-              ADNumberTraits<ADNumberType>::get_directional_derivative(
-                dependent_variables[0], i);
-
-            for (unsigned int j = 0; j <= i; ++j) // Symmetry
-              {
-                // Extract higher-order directional derivatives. Depending on
-                // the AD number type, the result may be another AD number or a
-                // floating point value.
-                const ScalarType hessian_ij =
-                  internal::NumberType<ScalarType>::value(
-                    ADNumberTraits<derivative_type>::get_directional_derivative(
-                      gradient_i, j));
-                hessian[i][j] = hessian_ij;
-                if (i != j)
-                  hessian[j][i] = hessian_ij; // Symmetry
-              }
-          }
-      }
-
-      // === Vector drivers ===
+              FullMatrix<ScalarType> &         hessian) const;
 
-      static void
-      values(const std::vector<ADNumberType> &dependent_variables,
-             Vector<ScalarType> &             values)
-      {
-        Assert(values.size() == dependent_variables.size(),
-               ExcDimensionMismatch(values.size(), dependent_variables.size()));
+      //@}
 
-        const std::size_t n_dependent_variables = dependent_variables.size();
-        for (unsigned int i = 0; i < n_dependent_variables; i++)
-          values[i] = ADNumberTraits<ADNumberType>::get_scalar_value(
-            dependent_variables[i]);
-      }
+      /**
+       * @name Drivers for vector functions
+       */
+      //@{
 
-      static void
+      void
+      values(const std::vector<ADNumberType> &dependent_variables,
+             Vector<ScalarType> &             values) const;
+
+      void
       jacobian(const std::vector<ADNumberType> &independent_variables,
                const std::vector<ADNumberType> &dependent_variables,
-               FullMatrix<ScalarType> &         jacobian)
-      {
-        Assert(
-          AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >= 1,
-          ExcSupportedDerivativeLevels(
-            AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
-            1));
-        Assert(jacobian.m() == dependent_variables.size(),
-               ExcDimensionMismatch(jacobian.m(), dependent_variables.size()));
-        Assert(jacobian.n() == independent_variables.size(),
-               ExcDimensionMismatch(jacobian.n(),
-                                    independent_variables.size()));
-
-        const std::size_t n_independent_variables =
-          independent_variables.size();
-        const std::size_t n_dependent_variables = dependent_variables.size();
-
-        // In forward mode, the gradients are computed from the
-        // dependent variables
-        for (unsigned int i = 0; i < n_dependent_variables; i++)
-          for (unsigned int j = 0; j < n_independent_variables; j++)
-            jacobian[i][j] = internal::NumberType<ScalarType>::value(
-              ADNumberTraits<ADNumberType>::get_directional_derivative(
-                dependent_variables[i], j));
-      }
+               FullMatrix<ScalarType> &         jacobian) const;
+
+      //@}
 
     private:
       /**
index 83e8a899c57120ec719af87c17868a60af80c5d8..7aad5de7edab5966a916a7e6979a2b39260ecdf8 100644 (file)
@@ -328,6 +328,9 @@ namespace Differentiation
        * image space.
        * @param[in] clear_registered_tapes A flag that indicates the that
        * list of @p registered_tapes must be cleared.
+       * If set to <tt>true</tt> then the data structure that tracks which
+       * tapes have been recorded is cleared as well. It is then expected that
+       * any preexisting tapes be re-recorded.
        *
        * @note This also resets the active tape number to an invalid number, and
        * deactivates the recording mode for taped variables.
index 81e83bb6339c6426ee2c905c2d67c64ed5c84461..9a2f31e4dc6770f5432a9a6e0cc797a08a19ed93 100644 (file)
@@ -1,6 +1,6 @@
 ## ---------------------------------------------------------------------
 ##
-## Copyright (C) 2017 by the deal.II authors
+## Copyright (C) 2017 - 2018 by the deal.II authors
 ##
 ## This file is part of the deal.II library.
 ##
 INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 
 SET(_src
+  ad_drivers.cc
   ad_helpers.cc
   adolc_number_types.cc
   sacado_number_types.cc
   )
 
 SET(_inst
+  ad_drivers.inst1.in
+  ad_drivers.inst2.in
   ad_helpers.inst1.in
   ad_helpers.inst2.in
   adolc_number_types.inst.in
diff --git a/source/differentiation/ad/ad_drivers.cc b/source/differentiation/ad/ad_drivers.cc
new file mode 100644 (file)
index 0000000..a0447c8
--- /dev/null
@@ -0,0 +1,2258 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/config.h>
+
+#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
+
+#  include <deal.II/base/exceptions.h>
+#  include <deal.II/base/types.h>
+#  include <deal.II/base/utilities.h>
+
+#  include <deal.II/differentiation/ad/ad_drivers.h>
+#  include <deal.II/differentiation/ad/ad_number_traits.h>
+#  include <deal.II/differentiation/ad/ad_number_types.h>
+#  include <deal.II/differentiation/ad/adolc_number_types.h>
+#  include <deal.II/differentiation/ad/sacado_number_types.h>
+
+#  include <deal.II/lac/full_matrix.h>
+#  include <deal.II/lac/vector.h>
+
+#  ifdef DEAL_II_WITH_ADOLC
+
+DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+#    include <adolc/adolc_fatalerror.h>
+#    include <adolc/drivers/drivers.h>
+#    include <adolc/taping.h>
+DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
+
+#  endif // DEAL_II_WITH_ADOLC
+
+#  include <vector>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace Differentiation
+{
+  namespace AD
+  {
+    // -------------   TapedDrivers   -------------
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    bool
+    TapedDrivers<ADNumberType, ScalarType, T>::is_recording() const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    typename Types<ADNumberType>::tape_index
+    TapedDrivers<ADNumberType, ScalarType, T>::active_tape_index() const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return Numbers<ADNumberType>::invalid_tape_index;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    bool
+    TapedDrivers<ADNumberType, ScalarType, T>::is_registered_tape(
+      const typename Types<ADNumberType>::tape_index) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    bool
+    TapedDrivers<ADNumberType, ScalarType, T>::keep_independent_values() const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::set_tape_buffer_sizes(
+      const typename Types<ADNumberType>::tape_buffer_sizes,
+      const typename Types<ADNumberType>::tape_buffer_sizes,
+      const typename Types<ADNumberType>::tape_buffer_sizes,
+      const typename Types<ADNumberType>::tape_buffer_sizes)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::start_taping(
+      const typename Types<ADNumberType>::tape_index,
+      const bool)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::stop_taping(
+      const typename Types<ADNumberType>::tape_index,
+      const bool)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    std::vector<typename Types<ADNumberType>::tape_index>
+    TapedDrivers<ADNumberType, ScalarType, T>::get_registered_tape_indices()
+      const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return std::vector<typename Types<ADNumberType>::tape_index>();
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::activate_tape(
+      const typename Types<ADNumberType>::tape_index)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    bool
+    TapedDrivers<ADNumberType, ScalarType, T>::requires_retaping(
+      const typename Types<ADNumberType>::tape_index) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    bool
+    TapedDrivers<ADNumberType, ScalarType, T>::last_action_requires_retaping()
+      const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::remove_tape(
+      const typename Types<ADNumberType>::tape_index)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::reset(const bool)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::print(std::ostream &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::print_tape_stats(
+      const typename Types<ADNumberType>::tape_index,
+      std::ostream &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    ScalarType
+    TapedDrivers<ADNumberType, ScalarType, T>::value(
+      const typename Types<ADNumberType>::tape_index,
+      const std::vector<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return ScalarType(0.0);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::gradient(
+      const typename Types<ADNumberType>::tape_index,
+      const std::vector<ScalarType> &,
+      Vector<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::hessian(
+      const typename Types<ADNumberType>::tape_index,
+      const std::vector<ScalarType> &,
+      FullMatrix<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::values(
+      const typename Types<ADNumberType>::tape_index,
+      const unsigned int,
+      const std::vector<ScalarType> &,
+      Vector<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType, ScalarType, T>::jacobian(
+      const typename Types<ADNumberType>::tape_index,
+      const unsigned int,
+      const std::vector<ScalarType> &,
+      FullMatrix<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+
+#  ifdef DEAL_II_WITH_ADOLC
+
+    // Specialization for taped ADOL-C auto-differentiable numbers.
+
+    template <typename ADNumberType>
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::TapedDrivers()
+      : active_tape(Numbers<ADNumberType>::invalid_tape_index)
+      , keep_values(true)
+      , is_recording_flag(false)
+      , use_stored_taped_buffer_sizes(false)
+      , obufsize(0u)
+      , lbufsize(0u)
+      , vbufsize(0u)
+      , tbufsize(0u)
+    {}
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::is_recording()
+      const
+    {
+      return is_recording_flag;
+    }
+
+
+    template <typename ADNumberType>
+    typename Types<ADNumberType>::tape_index
+    TapedDrivers<ADNumberType,
+                 double,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::active_tape_index() const
+    {
+      return active_tape;
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<ADNumberType,
+                 double,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::keep_independent_values()
+      const
+    {
+      return keep_values;
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      is_registered_tape(
+        const typename Types<ADNumberType>::tape_index tape_index) const
+    {
+      // Sigh... This is a mess :-/
+      // The most succinct way to get this piece of information, would be to
+      // use the getTapeInfos() function, but would come at the expense of
+      // creating an inactive tape data within ADOL-C's global store. For
+      // hints as to why this is the way it is, see the getTapeInfos()
+      // function in
+      // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/src/tape_handling.cpp
+      // An alternative solution would be to manually access the tape data;
+      // see the removeTape() function in
+      // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/src/tape_handling.cpp
+      // with the consideration of the #defines in
+      // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/src/taping_p.h
+      // as to how this would be performed.
+      // Doing things "manually" (the second way) without creating the
+      // additional data object would be a lot more work...
+      //
+      // Both of the above solutions would be possible IF ADOL-C exposed
+      // this data object or a method to access it to the outside world.
+      // But they don't :-(
+      // Instead, what we'll have to do is get the statistics for this tape
+      // and make our own determination as to whether or not this tape exists.
+      // This effectively executes the first solution, with even more
+      // overhead! If the tape is in existence, we can assume that it should a
+      // non-zero number of dependent and independent variables. Those are
+      // stored in the zeroth and first entries of the statistics vector.
+      //
+      // But, oh wait... Surprise! This will trigger an error if the tape
+      // doesn't exist at all! So lets first check their tape info cache to
+      // see if the tape REALLY exists (i.e. has been touched, even if nothing
+      // has been written to it) before trying to access it. It'll only take
+      // an O(n) search, but at this point who really cares about efficiency?
+      //
+      // It has been suggested in
+      // https://gitlab.com/adol-c/adol-c/issues/11
+      // that a simply try-catch block around ::tapestats is the solution that
+      // we want here. Unfortunately this results in unwanted pollution of
+      // the terminal, of the form
+      // ADOL-C error: reading integer tape number 4!
+      //               >>> File or directory not found! <<<
+      // , every time a query is made about a non-existant tape.
+      // So either way we have to guard that check with something more
+      // conervative so that we don't output useless messages for our users.
+      const std::vector<typename Types<ADNumberType>::tape_index>
+                 registered_tape_indices = get_registered_tape_indices();
+      const auto it = std::find(registered_tape_indices.begin(),
+                                registered_tape_indices.end(),
+                                tape_index);
+      if (it == registered_tape_indices.end())
+        return false;
+
+      // See https://gitlab.com/adol-c/adol-c/issues/11#note_108341333
+      try
+        {
+          std::vector<std::size_t> counts(STAT_SIZE);
+          ::tapestats(tape_index, counts.data());
+          return true;
+        }
+      catch (const ::FatalError &exc)
+        {
+          return false;
+        }
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      set_tape_buffer_sizes(
+        const typename Types<ADNumberType>::tape_buffer_sizes in_obufsize,
+        const typename Types<ADNumberType>::tape_buffer_sizes in_lbufsize,
+        const typename Types<ADNumberType>::tape_buffer_sizes in_vbufsize,
+        const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize)
+    {
+      // When valid for the chosen AD number type, these values will be used
+      // the next time start_recording_operations() is called.
+      obufsize                      = in_obufsize;
+      lbufsize                      = in_lbufsize;
+      vbufsize                      = in_vbufsize;
+      tbufsize                      = in_tbufsize;
+      use_stored_taped_buffer_sizes = true;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      start_taping(const typename Types<ADNumberType>::tape_index tape_index,
+                   const bool keep_independent_values)
+    {
+      if (use_stored_taped_buffer_sizes)
+        trace_on(tape_index,
+                 keep_independent_values,
+                 obufsize,
+                 lbufsize,
+                 vbufsize,
+                 tbufsize);
+      else
+        trace_on(tape_index, keep_independent_values);
+
+      // Set some other flags to their indicated / required values
+      keep_values       = keep_independent_values;
+      is_recording_flag = true;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      stop_taping(
+        const typename Types<ADNumberType>::tape_index active_tape_index,
+        const bool                                     write_tapes_to_file)
+    {
+      if (write_tapes_to_file)
+        trace_off(active_tape_index); // Slow
+      else
+        trace_off(); // Fast(er)
+
+      // Now that we've turned tracing off, we've definitely
+      // stopped all tape recording.
+      is_recording_flag = false;
+
+      // If the keep_values flag is set, then we expect the user to use this
+      // tape immediately after recording it. There is therefore no need to
+      // invalidate it. However, there is now also no way to double-check
+      // that the newly recorded tape is indeed the active tape.
+      if (keep_independent_values() == false)
+        active_tape = Numbers<ADNumberType>::invalid_tape_index;
+    }
+
+
+    template <typename ADNumberType>
+    std::vector<typename Types<ADNumberType>::tape_index>
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      get_registered_tape_indices() const
+    {
+      // We've chosen to use unsigned shorts for the tape
+      // index type (a safety precaution) so we need to
+      // perform a conversion betwwen ADOL-C's native tape
+      // index type and that chosen by us.
+      std::vector<short> registered_tape_indices_s;
+      cachedTraceTags(registered_tape_indices_s);
+
+      std::vector<typename Types<ADNumberType>::tape_index>
+        registered_tape_indices(registered_tape_indices_s.size());
+      std::copy(registered_tape_indices_s.begin(),
+                registered_tape_indices_s.end(),
+                registered_tape_indices.begin());
+      return registered_tape_indices;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      activate_tape(const typename Types<ADNumberType>::tape_index tape_index)
+    {
+      active_tape = tape_index;
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      requires_retaping(
+        const typename Types<ADNumberType>::tape_index tape_index) const
+    {
+      Assert(
+        is_registered_tape(tape_index) == true,
+        ExcMessage(
+          "Cannot ask for the status of a tape that has not yet been recorded and used."));
+
+      const auto it_status_tape = status.find(tape_index);
+
+      // This tape's status has not been found in the map. This could be
+      // because a non-existant tape_index has been used as an argument, or
+      // because the tape exists but has not been used (in a way that
+      // initiated a status update). For example, on can create a tape in one
+      // section of code and then query the status of all existing tapes in a
+      // completely different section of code that knows nothing about the
+      // first tape. There is no prerequisite that the first tape is ever
+      // used, and it can therefore not have a status. So, in this case
+      // there's not much we can do other than to report that the tape does
+      // not require retaping.
+      if (it_status_tape == status.end())
+        return false;
+
+      const auto status_tape = it_status_tape->second;
+
+      // See ADOL-C manual section 1.7 and comments in last paragraph of
+      // section 3.1
+      Assert(status_tape < 4 && status_tape >= -2,
+             ExcIndexRange(status_tape, -2, 4));
+      return (status_tape < 0);
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      last_action_requires_retaping() const
+    {
+      return requires_retaping(active_tape);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      remove_tape(const typename Types<ADNumberType>::tape_index tape_index)
+    {
+      Assert(is_registered_tape(tape_index),
+             ExcMessage(
+               "This tape does not exist, and therefore cannot be cleared."));
+      removeTape(tape_index, TapeRemovalType::ADOLC_REMOVE_COMPLETELY);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      reset(const bool clear_registered_tapes)
+    {
+      active_tape       = Numbers<ADNumberType>::invalid_tape_index;
+      is_recording_flag = false;
+      status.clear();
+      if (clear_registered_tapes)
+        {
+          const std::vector<typename Types<ADNumberType>::tape_index>
+            registered_tape_indices = get_registered_tape_indices();
+          for (const auto &tape_index : registered_tape_indices)
+            remove_tape(tape_index);
+        }
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<ADNumberType,
+                 double,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::print(std::ostream &stream)
+      const
+    {
+      const std::vector<typename Types<ADNumberType>::tape_index>
+        registered_tape_indices = get_registered_tape_indices();
+      stream << "Registered tapes and their status: ";
+      auto it_registered_tape = registered_tape_indices.begin();
+      for (unsigned int i = 0; i < registered_tape_indices.size();
+           ++i, ++it_registered_tape)
+        {
+          const auto tape_index     = *it_registered_tape;
+          const auto it_status_tape = status.find(tape_index);
+          Assert(it_status_tape != status.end(),
+                 ExcMessage(
+                   "This tape's status has not been found in the map."));
+          const auto status_tape = it_status_tape->second;
+
+          stream << tape_index << "->" << status_tape
+                 << (i < (registered_tape_indices.size() - 1) ? "," : "");
+        }
+      stream << "\n";
+
+      stream << "Keep values? " << keep_independent_values() << "\n";
+      stream << "Use stored tape buffer sizes? "
+             << use_stored_taped_buffer_sizes << "\n";
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      print_tape_stats(
+        const typename Types<ADNumberType>::tape_index tape_index,
+        std::ostream &                                 stream) const
+    {
+      // See ADOL-C manual section 2.1
+      // and adolc/taping.h
+      std::vector<std::size_t> counts(STAT_SIZE);
+      ::tapestats(tape_index, counts.data());
+      Assert(counts.size() >= 18, ExcInternalError());
+      stream
+        << "Tape index: " << tape_index << "\n"
+        << "Number of independent variables: " << counts[0] << "\n"
+        << "Number of dependent variables:   " << counts[1] << "\n"
+        << "Max number of live, active variables: " << counts[2] << "\n"
+        << "Size of taylor stack (number of overwrites): " << counts[3] << "\n"
+        << "Operations buffer size: " << counts[4] << "\n"
+        << "Total number of recorded operations: " << counts[5] << "\n"
+        << "Operations file written or not: " << counts[6] << "\n"
+        << "Overall number of locations: " << counts[7] << "\n"
+        << "Locations file written or not: " << counts[8] << "\n"
+        << "Overall number of values: " << counts[9] << "\n"
+        << "Values file written or not: " << counts[10] << "\n"
+        << "Locations buffer size: " << counts[11] << "\n"
+        << "Values buffer size: " << counts[12] << "\n"
+        << "Taylor buffer size: " << counts[13] << "\n"
+        << "Number of eq_*_prod for sparsity pattern: " << counts[14] << "\n"
+        << "Use of 'min_op', deferred to 'abs_op' for piecewise calculations: "
+        << counts[15] << "\n"
+        << "Number of 'abs' calls that can switch branch: " << counts[16]
+        << "\n"
+        << "Number of parameters (doubles) interchangeable without retaping: "
+        << counts[17] << "\n"
+        << std::flush;
+    }
+
+
+    template <typename ADNumberType>
+    typename TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::scalar_type
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      value(const typename Types<ADNumberType>::tape_index active_tape_index,
+            const std::vector<scalar_type> &independent_variables) const
+    {
+      Assert(is_registered_tape(active_tape_index),
+             ExcMessage("This tape has not yet been recorded."));
+
+      scalar_type value = 0.0;
+
+      status[active_tape_index] =
+        ::function(active_tape_index,
+                   1, // Only one dependent variable
+                   independent_variables.size(),
+                   const_cast<double *>(independent_variables.data()),
+                   &value);
+
+      return value;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
+               const std::vector<scalar_type> &independent_variables,
+               Vector<scalar_type> &           gradient) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               1,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               1));
+      Assert(gradient.size() == independent_variables.size(),
+             ExcDimensionMismatch(gradient.size(),
+                                  independent_variables.size()));
+      Assert(is_registered_tape(active_tape_index),
+             ExcMessage("This tape has not yet been recorded."));
+
+      // Note: ADOL-C's ::gradient function expects a *double as the last
+      // parameter. Here we take advantage of the fact that the data in the
+      // Vector class is aligned (e.g. stored as an Array)
+      status[active_tape_index] =
+        ::gradient(active_tape_index,
+                   independent_variables.size(),
+                   const_cast<scalar_type *>(independent_variables.data()),
+                   gradient.data());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
+              const std::vector<scalar_type> &independent_variables,
+              FullMatrix<scalar_type> &       hessian) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               2,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               2));
+      Assert(hessian.m() == independent_variables.size(),
+             ExcDimensionMismatch(hessian.m(), independent_variables.size()));
+      Assert(hessian.n() == independent_variables.size(),
+             ExcDimensionMismatch(hessian.n(), independent_variables.size()));
+      Assert(is_registered_tape(active_tape_index),
+             ExcMessage("This tape has not yet been recorded."));
+
+      const unsigned int n_independent_variables = independent_variables.size();
+      std::vector<scalar_type *> H(n_independent_variables);
+      for (unsigned int i = 0; i < n_independent_variables; ++i)
+        H[i] = &hessian[i][0];
+
+      status[active_tape_index] =
+        ::hessian(active_tape_index,
+                  n_independent_variables,
+                  const_cast<scalar_type *>(independent_variables.data()),
+                  H.data());
+
+      // ADOL-C builds only the lower-triangular part of the
+      // symmetric Hessian, so we should copy the relevant
+      // entries into the upper triangular part.
+      for (unsigned int i = 0; i < n_independent_variables; i++)
+        for (unsigned int j = 0; j < i; j++)
+          hessian[j][i] = hessian[i][j]; // Symmetry
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      values(const typename Types<ADNumberType>::tape_index active_tape_index,
+             const unsigned int              n_dependent_variables,
+             const std::vector<scalar_type> &independent_variables,
+             Vector<scalar_type> &           values) const
+    {
+      Assert(values.size() == n_dependent_variables,
+             ExcDimensionMismatch(values.size(), n_dependent_variables));
+      Assert(is_registered_tape(active_tape_index),
+             ExcMessage("This tape has not yet been recorded."));
+
+      // Note: ADOL-C's ::function function expects a *double as the last
+      // parameter. Here we take advantage of the fact that the data in the
+      // Vector class is aligned (e.g. stored as an Array)
+      status[active_tape_index] =
+        ::function(active_tape_index,
+                   n_dependent_variables,
+                   independent_variables.size(),
+                   const_cast<scalar_type *>(independent_variables.data()),
+                   values.data());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
+               const unsigned int              n_dependent_variables,
+               const std::vector<scalar_type> &independent_variables,
+               FullMatrix<scalar_type> &       jacobian) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               1,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               1));
+      Assert(jacobian.m() == n_dependent_variables,
+             ExcDimensionMismatch(jacobian.m(), n_dependent_variables));
+      Assert(jacobian.n() == independent_variables.size(),
+             ExcDimensionMismatch(jacobian.n(), independent_variables.size()));
+      Assert(is_registered_tape(active_tape_index),
+             ExcMessage("This tape has not yet been recorded."));
+
+      std::vector<scalar_type *> J(n_dependent_variables);
+      for (unsigned int i = 0; i < n_dependent_variables; ++i)
+        J[i] = &jacobian[i][0];
+
+      status[active_tape_index] = ::jacobian(active_tape_index,
+                                             n_dependent_variables,
+                                             independent_variables.size(),
+                                             independent_variables.data(),
+                                             J.data());
+    }
+
+#  else
+
+    // Specialization for taped ADOL-C auto-differentiable numbers.
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::is_recording()
+      const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return false;
+    }
+
+
+    template <typename ADNumberType>
+    typename Types<ADNumberType>::tape_index
+    TapedDrivers<ADNumberType,
+                 double,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::active_tape_index() const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return Numbers<ADNumberType>::invalid_tape_index;
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<ADNumberType,
+                 double,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::keep_independent_values()
+      const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return false;
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      is_registered_tape(
+        const typename Types<ADNumberType>::tape_index tape_index) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return false;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      set_tape_buffer_sizes(
+        const typename Types<ADNumberType>::tape_buffer_sizes,
+        const typename Types<ADNumberType>::tape_buffer_sizes,
+        const typename Types<ADNumberType>::tape_buffer_sizes,
+        const typename Types<ADNumberType>::tape_buffer_sizes)
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      start_taping(const typename Types<ADNumberType>::tape_index, const bool)
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      stop_taping(const typename Types<ADNumberType>::tape_index, const bool)
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+    std::vector<typename Types<ADNumberType>::tape_index>
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      get_registered_tape_indices() const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return std::vector<typename Types<ADNumberType>::tape_index>();
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      activate_tape(const typename Types<ADNumberType>::tape_index)
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      requires_retaping(const typename Types<ADNumberType>::tape_index) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return false;
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      last_action_requires_retaping() const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return false;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      remove_tape(const typename Types<ADNumberType>::tape_index)
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<ADNumberType,
+                 double,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::reset(const bool)
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<ADNumberType,
+                 double,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::print(std::ostream &stream)
+      const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      print_tape_stats(const typename Types<ADNumberType>::tape_index,
+                       std::ostream &) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::scalar_type
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      value(const typename Types<ADNumberType>::tape_index,
+            const std::vector<scalar_type> &) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+      return 0.0;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      gradient(const typename Types<ADNumberType>::tape_index,
+               const std::vector<scalar_type> &,
+               Vector<scalar_type> &) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      hessian(const typename Types<ADNumberType>::tape_index,
+              const std::vector<scalar_type> &,
+              FullMatrix<scalar_type> &) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      values(const typename Types<ADNumberType>::tape_index,
+             const unsigned int,
+             const std::vector<scalar_type> &,
+             Vector<scalar_type> &) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      double,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      jacobian(const typename Types<ADNumberType>::tape_index,
+               const unsigned int,
+               const std::vector<scalar_type> &,
+               FullMatrix<scalar_type> &) const
+    {
+      AssertThrow(false, ExcRequiresADOLC());
+    }
+
+#  endif // DEAL_II_WITH_ADOLC
+
+
+    // Specialization for ADOL-C taped numbers. It is expected that the
+    // scalar return type for this class is a float.
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::is_recording()
+      const
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      return taped_driver.is_recording();
+    }
+
+
+    template <typename ADNumberType>
+    typename Types<ADNumberType>::tape_index
+    TapedDrivers<ADNumberType,
+                 float,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::active_tape_index() const
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      return taped_driver.active_tape_index();
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<ADNumberType,
+                 float,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::keep_independent_values()
+      const
+    {
+      return taped_driver.keep_independent_values();
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      is_registered_tape(
+        const typename Types<ADNumberType>::tape_index tape_index) const
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      return taped_driver.is_registered_tape(tape_index);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      set_tape_buffer_sizes(
+        const typename Types<ADNumberType>::tape_buffer_sizes obufsize,
+        const typename Types<ADNumberType>::tape_buffer_sizes lbufsize,
+        const typename Types<ADNumberType>::tape_buffer_sizes vbufsize,
+        const typename Types<ADNumberType>::tape_buffer_sizes tbufsize)
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.set_tape_buffer_sizes(obufsize,
+                                         lbufsize,
+                                         vbufsize,
+                                         tbufsize);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      start_taping(const typename Types<ADNumberType>::tape_index tape_index,
+                   const bool keep_independent_values)
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.start_taping(tape_index, keep_independent_values);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      stop_taping(
+        const typename Types<ADNumberType>::tape_index active_tape_index,
+        const bool                                     write_tapes_to_file)
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.stop_taping(active_tape_index, write_tapes_to_file);
+    }
+
+
+    template <typename ADNumberType>
+    std::vector<typename Types<ADNumberType>::tape_index>
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      get_registered_tape_indices() const
+    {
+      return taped_driver.get_registered_tape_indices();
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      activate_tape(const typename Types<ADNumberType>::tape_index tape_index)
+    {
+      taped_driver.activate_tape(tape_index);
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      requires_retaping(
+        const typename Types<ADNumberType>::tape_index tape_index) const
+    {
+      return taped_driver.requires_retaping(tape_index);
+    }
+
+
+    template <typename ADNumberType>
+    bool
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      last_action_requires_retaping() const
+    {
+      return taped_driver.last_action_requires_retaping();
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      remove_tape(const typename Types<ADNumberType>::tape_index tape_index)
+    {
+      taped_driver.remove_tape(tape_index);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      reset(const bool clear_registered_tapes)
+    {
+      taped_driver.reset(clear_registered_tapes);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<ADNumberType,
+                 float,
+                 typename std::enable_if<
+                   ADNumberTraits<ADNumberType>::type_code ==
+                   NumberTypes::adolc_taped>::type>::print(std::ostream &stream)
+      const
+    {
+      taped_driver.print(stream);
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      print_tape_stats(
+        const typename Types<ADNumberType>::tape_index tape_index,
+        std::ostream &                                 stream) const
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.print_tape_stats(tape_index, stream);
+    }
+
+
+    template <typename ADNumberType>
+    typename TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::scalar_type
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      value(const typename Types<ADNumberType>::tape_index active_tape_index,
+            const std::vector<scalar_type> &independent_variables) const
+    {
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      return taped_driver.value(active_tape_index,
+                                vector_float_to_double(independent_variables));
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
+               const std::vector<scalar_type> &independent_variables,
+               Vector<scalar_type> &           gradient) const
+    {
+      Vector<double> gradient_double(gradient.size());
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.gradient(active_tape_index,
+                            vector_float_to_double(independent_variables),
+                            gradient_double);
+      gradient = gradient_double;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
+              const std::vector<scalar_type> &independent_variables,
+              FullMatrix<scalar_type> &       hessian) const
+    {
+      FullMatrix<double> hessian_double(hessian.m(), hessian.n());
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.hessian(active_tape_index,
+                           vector_float_to_double(independent_variables),
+                           hessian_double);
+      hessian = hessian_double;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      values(const typename Types<ADNumberType>::tape_index active_tape_index,
+             const unsigned int              n_dependent_variables,
+             const std::vector<scalar_type> &independent_variables,
+             Vector<scalar_type> &           values) const
+    {
+      Vector<double> values_double(values.size());
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.values(active_tape_index,
+                          n_dependent_variables,
+                          vector_float_to_double(independent_variables),
+                          values_double);
+      values = values_double;
+    }
+
+
+    template <typename ADNumberType>
+    void
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
+               const unsigned int              n_dependent_variables,
+               const std::vector<scalar_type> &independent_variables,
+               FullMatrix<scalar_type> &       jacobian) const
+    {
+      FullMatrix<double> jacobian_double(jacobian.m(), jacobian.n());
+      // ADOL-C only supports 'double', not 'float', so we can forward to
+      // the 'double' implementation of this function
+      taped_driver.jacobian(active_tape_index,
+                            n_dependent_variables,
+                            vector_float_to_double(independent_variables),
+                            jacobian_double);
+      jacobian = jacobian_double;
+    }
+
+
+    template <typename ADNumberType>
+    std::vector<double>
+    TapedDrivers<
+      ADNumberType,
+      float,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_taped>::type>::
+      vector_float_to_double(const std::vector<float> &in) const
+    {
+      std::vector<double> out(in.size());
+      std::copy(in.begin(), in.end(), out.begin());
+      return out;
+    }
+
+
+    // -------------   TapelessDrivers   -------------
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType, ScalarType, T>::initialize_global_environment(
+      const unsigned int)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType, ScalarType, T>::
+      allow_dependent_variable_marking()
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType, ScalarType, T>::
+      prevent_dependent_variable_marking()
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    bool
+    TapelessDrivers<ADNumberType, ScalarType, T>::
+      is_dependent_variable_marking_allowed() const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    ScalarType
+    TapelessDrivers<ADNumberType, ScalarType, T>::value(
+      const std::vector<ADNumberType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+      return ScalarType(0.0);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType, ScalarType, T>::gradient(
+      const std::vector<ADNumberType> &,
+      const std::vector<ADNumberType> &,
+      Vector<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType, ScalarType, T>::hessian(
+      const std::vector<ADNumberType> &,
+      const std::vector<ADNumberType> &,
+      FullMatrix<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType, ScalarType, T>::values(
+      const std::vector<ADNumberType> &,
+      Vector<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template <typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType, ScalarType, T>::jacobian(
+      const std::vector<ADNumberType> &,
+      const std::vector<ADNumberType> &,
+      FullMatrix<ScalarType> &) const
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    namespace internal
+    {
+      /**
+       * A dummy function to define the active dependent variable when using
+       * reverse-mode AD.
+       */
+      template <typename ADNumberType>
+      typename 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 &)
+      {}
+
+#  ifdef DEAL_II_TRILINOS_WITH_SACADO
+
+
+      /**
+       * Define the active dependent variable when using reverse-mode AD.
+       *
+       * If there are multiple dependent variables then it is necessary to
+       * inform the independent variables, from which the adjoints are computed,
+       * which dependent variable they are computing the gradients with respect
+       * to. This function broadcasts this information.
+       */
+      template <typename ADNumberType>
+      typename 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)
+      {
+        // Compute all gradients (adjoints) for this
+        // reverse-mode Sacado dependent variable.
+        // For reverse-mode Sacado numbers it is necessary to broadcast to
+        // all independent variables that it is time to compute gradients.
+        // For one dependent variable one would just need to call
+        // ADNumberType::Gradcomp(), but since we have a more
+        // generic implementation for vectors of dependent variables
+        // (vector mode) we default to the complex case.
+        ADNumberType::Outvar_Gradcomp(dependent_variable);
+      }
+
+#  endif
+
+
+      /**
+       * A dummy function to enable vector mode for tapeless
+       * auto-differentiable numbers.
+       */
+      template <typename ADNumberType>
+      typename std::enable_if<!(ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless)>::type
+      configure_tapeless_mode(const unsigned int)
+      {}
+
+
+#  ifdef DEAL_II_WITH_ADOLC
+
+
+      /**
+       * Enable vector mode for ADOL-C tapeless numbers.
+       *
+       * This function checks to see if its legal to increase the maximum
+       * number of directional derivatives to be considered during calculations.
+       * If not then it throws an error.
+       */
+      template <typename ADNumberType>
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_tapeless>::type
+      configure_tapeless_mode(const unsigned int n_directional_derivatives)
+      {
+#    ifdef DEAL_II_ADOLC_WITH_TAPELESS_REFCOUNTING
+        // See ADOL-C manual section 7.1
+        //
+        // NOTE: It is critical that this is done for tapeless mode BEFORE
+        // any adtl::adouble are created. If this is not done, then we see
+        // this scary warning:
+        //
+        // "
+        // ADOL-C Warning: Tapeless: Setting numDir could change memory
+        // allocation of derivatives in existing adoubles and may lead to
+        // erroneous results or memory corruption
+        // "
+        //
+        // So we use this dummy function to configure this setting before
+        // we create and initialize our class data
+        const std::size_t n_live_variables = adtl::refcounter::getNumLiveVar();
+        if (n_live_variables == 0)
+          {
+            adtl::setNumDir(n_directional_derivatives);
+          }
+        else
+          {
+            // So there are some live active variables floating around. Here we
+            // check if we ask to increase the number of computable
+            // directional derivatives. If this really is necessary then it's
+            // absolutely vital that there exist no live variables before doing
+            // so.
+            const std::size_t n_set_directional_derivatives = adtl::getNumDir();
+            if (n_directional_derivatives > n_set_directional_derivatives)
+              AssertThrow(
+                n_live_variables == 0,
+                ExcMessage(
+                  "There are currently " +
+                  Utilities::to_string(n_live_variables) +
+                  " live "
+                  "adtl::adouble variables in existence. They currently "
+                  "assume " +
+                  Utilities::to_string(n_set_directional_derivatives) +
+                  " directional derivatives "
+                  "but you wish to increase this to " +
+                  Utilities::to_string(n_directional_derivatives) +
+                  ". \n"
+                  "To safely change (or more specifically in this case, "
+                  "increase) the number of directional derivatives, there "
+                  "must be no tapeless doubles in local/global scope."));
+          }
+#    else
+        // If ADOL-C is not configured with tapeless number reference counting
+        // then there is no way that we can guarantee that the following call
+        // is safe. No comment... :-/
+        adtl::setNumDir(n_directional_derivatives);
+#    endif
+      }
+
+#  else // DEAL_II_WITH_ADOLC
+
+      template <typename ADNumberType>
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                              NumberTypes::adolc_tapeless>::type
+      configure_tapeless_mode(const unsigned int /*n_directional_derivatives*/)
+      {
+        AssertThrow(false, ExcRequiresADOLC());
+      }
+
+#  endif
+
+    } // namespace internal
+
+
+
+    // Specialization for auto-differentiable numbers that use
+    // reverse mode to compute the first derivatives (and, if supported,
+    // forward mode for the second).
+
+    template <typename ADNumberType, typename ScalarType>
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<
+        ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_rad ||
+        ADNumberTraits<ADNumberType>::type_code ==
+          NumberTypes::sacado_rad_dfad>::type>::TapelessDrivers()
+      : dependent_variable_marking_safe(false)
+    {}
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      initialize_global_environment(const unsigned int n_independent_variables)
+    {
+      internal::configure_tapeless_mode<ADNumberType>(n_independent_variables);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      allow_dependent_variable_marking()
+    {
+      dependent_variable_marking_safe = true;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      prevent_dependent_variable_marking()
+    {
+      dependent_variable_marking_safe = false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    bool
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      is_dependent_variable_marking_allowed() const
+    {
+      return dependent_variable_marking_safe;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    ScalarType
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      value(const std::vector<ADNumberType> &dependent_variables) const
+    {
+      Assert(dependent_variables.size() == 1,
+             ExcDimensionMismatch(dependent_variables.size(), 1));
+      return ADNumberTraits<ADNumberType>::get_scalar_value(
+        dependent_variables[0]);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      gradient(const std::vector<ADNumberType> &independent_variables,
+               const std::vector<ADNumberType> &dependent_variables,
+               Vector<ScalarType> &             gradient) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               1,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               1));
+      Assert(dependent_variables.size() == 1,
+             ExcDimensionMismatch(dependent_variables.size(), 1));
+      Assert(gradient.size() == independent_variables.size(),
+             ExcDimensionMismatch(gradient.size(),
+                                  independent_variables.size()));
+
+      // In reverse mode, the gradients are computed from the
+      // independent variables (i.e. the adjoint)
+      internal::reverse_mode_dependent_variable_activation(
+        const_cast<ADNumberType &>(dependent_variables[0]));
+      const std::size_t n_independent_variables = independent_variables.size();
+      for (unsigned int i = 0; i < n_independent_variables; i++)
+        gradient[i] = internal::NumberType<ScalarType>::value(
+          ADNumberTraits<ADNumberType>::get_directional_derivative(
+            independent_variables[i], 0 /*This number doesn't really matter*/));
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      hessian(const std::vector<ADNumberType> &independent_variables,
+              const std::vector<ADNumberType> &dependent_variables,
+              FullMatrix<ScalarType> &         hessian) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               2,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               2));
+      Assert(dependent_variables.size() == 1,
+             ExcDimensionMismatch(dependent_variables.size(), 1));
+      Assert(hessian.m() == independent_variables.size(),
+             ExcDimensionMismatch(hessian.m(), independent_variables.size()));
+      Assert(hessian.n() == independent_variables.size(),
+             ExcDimensionMismatch(hessian.n(), independent_variables.size()));
+
+      // In reverse mode, the gradients are computed from the
+      // independent variables (i.e. the adjoint)
+      internal::reverse_mode_dependent_variable_activation(
+        const_cast<ADNumberType &>(dependent_variables[0]));
+      const std::size_t n_independent_variables = independent_variables.size();
+      for (unsigned int i = 0; i < n_independent_variables; i++)
+        {
+          using derivative_type =
+            typename ADNumberTraits<ADNumberType>::derivative_type;
+          const derivative_type gradient_i =
+            ADNumberTraits<ADNumberType>::get_directional_derivative(
+              independent_variables[i], i);
+
+          for (unsigned int j = 0; j <= i; ++j) // Symmetry
+            {
+              // Extract higher-order directional derivatives. Depending on
+              // the AD number type, the result may be another AD number or a
+              // floating point value.
+              const ScalarType hessian_ij =
+                internal::NumberType<ScalarType>::value(
+                  ADNumberTraits<derivative_type>::get_directional_derivative(
+                    gradient_i, j));
+              hessian[i][j] = hessian_ij;
+              if (i != j)
+                hessian[j][i] = hessian_ij; // Symmetry
+            }
+        }
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      values(const std::vector<ADNumberType> &dependent_variables,
+             Vector<ScalarType> &             values) const
+    {
+      Assert(values.size() == dependent_variables.size(),
+             ExcDimensionMismatch(values.size(), dependent_variables.size()));
+
+      const std::size_t n_dependent_variables = dependent_variables.size();
+      for (unsigned int i = 0; i < n_dependent_variables; i++)
+        values[i] = ADNumberTraits<ADNumberType>::get_scalar_value(
+          dependent_variables[i]);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_rad_dfad>::type>::
+      jacobian(const std::vector<ADNumberType> &independent_variables,
+               const std::vector<ADNumberType> &dependent_variables,
+               FullMatrix<ScalarType> &         jacobian) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               1,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               1));
+      Assert(jacobian.m() == dependent_variables.size(),
+             ExcDimensionMismatch(jacobian.m(), dependent_variables.size()));
+      Assert(jacobian.n() == independent_variables.size(),
+             ExcDimensionMismatch(jacobian.n(), independent_variables.size()));
+
+      const std::size_t n_independent_variables = independent_variables.size();
+      const std::size_t n_dependent_variables   = dependent_variables.size();
+
+      // In reverse mode, the gradients are computed from the
+      // independent variables (i.e. the adjoint).
+      // For a demonstration of why this accumulation process is
+      // required, see the unit tests
+      // sacado/basic_01b.cc and sacado/basic_02b.cc
+      // Here we also take into consideration the derivative type:
+      // The Sacado number may be of the nested variety, in which
+      // case the effect of the accumulation process on the
+      // sensitivities of the nested number need to be accounted for.
+      using accumulation_type =
+        typename ADNumberTraits<ADNumberType>::derivative_type;
+      std::vector<accumulation_type> rad_accumulation(
+        n_independent_variables,
+        dealii::internal::NumberType<accumulation_type>::value(0.0));
+
+      for (unsigned int i = 0; i < n_dependent_variables; i++)
+        {
+          internal::reverse_mode_dependent_variable_activation(
+            const_cast<ADNumberType &>(dependent_variables[i]));
+          for (unsigned int j = 0; j < n_independent_variables; j++)
+            {
+              const accumulation_type df_i_dx_j =
+                ADNumberTraits<ADNumberType>::get_directional_derivative(
+                  independent_variables[j],
+                  i /*This number doesn't really matter*/) -
+                rad_accumulation[j];
+              jacobian[i][j] =
+                internal::NumberType<ScalarType>::value(df_i_dx_j);
+              rad_accumulation[j] += df_i_dx_j;
+            }
+        }
+    }
+
+
+
+    // Specialization for auto-differentiable numbers that use
+    // forward mode to compute the first (and, if supported, second)
+    // derivatives.
+
+    template <typename ADNumberType, typename ScalarType>
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<
+        ADNumberTraits<ADNumberType>::type_code ==
+          NumberTypes::adolc_tapeless ||
+        ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_dfad ||
+        ADNumberTraits<ADNumberType>::type_code ==
+          NumberTypes::sacado_dfad_dfad>::type>::TapelessDrivers()
+      : dependent_variable_marking_safe(false)
+    {}
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      initialize_global_environment(const unsigned int n_independent_variables)
+    {
+      internal::configure_tapeless_mode<ADNumberType>(n_independent_variables);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      allow_dependent_variable_marking()
+    {
+      dependent_variable_marking_safe = true;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      prevent_dependent_variable_marking()
+    {
+      dependent_variable_marking_safe = false;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    bool
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      is_dependent_variable_marking_allowed() const
+    {
+      return dependent_variable_marking_safe;
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    ScalarType
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      value(const std::vector<ADNumberType> &dependent_variables) const
+    {
+      Assert(dependent_variables.size() == 1,
+             ExcDimensionMismatch(dependent_variables.size(), 1));
+      return ADNumberTraits<ADNumberType>::get_scalar_value(
+        dependent_variables[0]);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      gradient(const std::vector<ADNumberType> &independent_variables,
+               const std::vector<ADNumberType> &dependent_variables,
+               Vector<ScalarType> &             gradient) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               1,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               1));
+      Assert(dependent_variables.size() == 1,
+             ExcDimensionMismatch(dependent_variables.size(), 1));
+      Assert(gradient.size() == independent_variables.size(),
+             ExcDimensionMismatch(gradient.size(),
+                                  independent_variables.size()));
+
+      // In forward mode, the gradients are computed from the
+      // dependent variables
+      const std::size_t n_independent_variables = independent_variables.size();
+      for (unsigned int i = 0; i < n_independent_variables; i++)
+        gradient[i] = internal::NumberType<ScalarType>::value(
+          ADNumberTraits<ADNumberType>::get_directional_derivative(
+            dependent_variables[0], i));
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      hessian(const std::vector<ADNumberType> &independent_variables,
+              const std::vector<ADNumberType> &dependent_variables,
+              FullMatrix<ScalarType> &         hessian) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               2,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               2));
+      Assert(dependent_variables.size() == 1,
+             ExcDimensionMismatch(dependent_variables.size(), 1));
+      Assert(hessian.m() == independent_variables.size(),
+             ExcDimensionMismatch(hessian.m(), independent_variables.size()));
+      Assert(hessian.n() == independent_variables.size(),
+             ExcDimensionMismatch(hessian.n(), independent_variables.size()));
+
+      // In forward mode, the gradients are computed from the
+      // dependent variables
+      const std::size_t n_independent_variables = independent_variables.size();
+      for (unsigned int i = 0; i < n_independent_variables; i++)
+        {
+          using derivative_type =
+            typename ADNumberTraits<ADNumberType>::derivative_type;
+          const derivative_type gradient_i =
+            ADNumberTraits<ADNumberType>::get_directional_derivative(
+              dependent_variables[0], i);
+
+          for (unsigned int j = 0; j <= i; ++j) // Symmetry
+            {
+              // Extract higher-order directional derivatives. Depending on
+              // the AD number type, the result may be another AD number or a
+              // floating point value.
+              const ScalarType hessian_ij =
+                internal::NumberType<ScalarType>::value(
+                  ADNumberTraits<derivative_type>::get_directional_derivative(
+                    gradient_i, j));
+              hessian[i][j] = hessian_ij;
+              if (i != j)
+                hessian[j][i] = hessian_ij; // Symmetry
+            }
+        }
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      values(const std::vector<ADNumberType> &dependent_variables,
+             Vector<ScalarType> &             values) const
+    {
+      Assert(values.size() == dependent_variables.size(),
+             ExcDimensionMismatch(values.size(), dependent_variables.size()));
+
+      const std::size_t n_dependent_variables = dependent_variables.size();
+      for (unsigned int i = 0; i < n_dependent_variables; i++)
+        values[i] = ADNumberTraits<ADNumberType>::get_scalar_value(
+          dependent_variables[i]);
+    }
+
+
+    template <typename ADNumberType, typename ScalarType>
+    void
+    TapelessDrivers<
+      ADNumberType,
+      ScalarType,
+      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::adolc_tapeless ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad ||
+                              ADNumberTraits<ADNumberType>::type_code ==
+                                NumberTypes::sacado_dfad_dfad>::type>::
+      jacobian(const std::vector<ADNumberType> &independent_variables,
+               const std::vector<ADNumberType> &dependent_variables,
+               FullMatrix<ScalarType> &         jacobian) const
+    {
+      Assert(AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels >=
+               1,
+             ExcSupportedDerivativeLevels(
+               AD::ADNumberTraits<ADNumberType>::n_supported_derivative_levels,
+               1));
+      Assert(jacobian.m() == dependent_variables.size(),
+             ExcDimensionMismatch(jacobian.m(), dependent_variables.size()));
+      Assert(jacobian.n() == independent_variables.size(),
+             ExcDimensionMismatch(jacobian.n(), independent_variables.size()));
+
+      const std::size_t n_independent_variables = independent_variables.size();
+      const std::size_t n_dependent_variables   = dependent_variables.size();
+
+      // In forward mode, the gradients are computed from the
+      // dependent variables
+      for (unsigned int i = 0; i < n_dependent_variables; i++)
+        for (unsigned int j = 0; j < n_independent_variables; j++)
+          jacobian[i][j] = internal::NumberType<ScalarType>::value(
+            ADNumberTraits<ADNumberType>::get_directional_derivative(
+              dependent_variables[i], j));
+    }
+
+
+  } // namespace AD
+} // namespace Differentiation
+
+
+/* --- Explicit instantiations --- */
+#  ifdef DEAL_II_WITH_ADOLC
+#    include "ad_drivers.inst1"
+#  endif
+#  ifdef DEAL_II_TRILINOS_WITH_SACADO
+#    include "ad_drivers.inst2"
+#  endif
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
diff --git a/source/differentiation/ad/ad_drivers.inst1.in b/source/differentiation/ad/ad_drivers.inst1.in
new file mode 100644 (file)
index 0000000..f082939
--- /dev/null
@@ -0,0 +1,73 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// TODO: Include complex types
+
+for (number : REAL_SCALARS)
+{
+  namespace Differentiation
+  \{
+  namespace AD
+  \{
+
+    // -------------------------- TapedDrivers ----------------------
+
+    template
+    struct TapedDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::adolc_taped>::ad_type,number>;
+
+    template
+    struct TapedDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::adolc_tapeless>::ad_type,number>;
+
+    // -------------------------- TapelessDrivers ----------------------
+
+    template
+    struct TapelessDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::adolc_taped>::ad_type,number>;
+
+    template
+    struct TapelessDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::adolc_tapeless>::ad_type,number>;
+
+    \}
+    \}
+}
+
+// Instantiations for ADHelpers for which the underlying number type is fixed
+for ()
+{
+    namespace Differentiation
+    \{
+    namespace AD
+    \{
+
+
+
+    // -------------------------- Types ----------------------
+
+    template
+    struct Types<typename Differentiation::AD::NumberTraits<double,NumberTypes::adolc_taped>::ad_type>;
+
+    template
+    struct Types<typename Differentiation::AD::NumberTraits<double,NumberTypes::adolc_tapeless>::ad_type>;
+
+    // -------------------------- Numbers ----------------------
+
+    template
+    struct Numbers<typename Differentiation::AD::NumberTraits<double,NumberTypes::adolc_taped>::ad_type>;
+
+    template
+    struct Numbers<typename Differentiation::AD::NumberTraits<double,NumberTypes::adolc_tapeless>::ad_type>;
+
+    \}
+    \}
+}
\ No newline at end of file
diff --git a/source/differentiation/ad/ad_drivers.inst2.in b/source/differentiation/ad/ad_drivers.inst2.in
new file mode 100644 (file)
index 0000000..f8aa39e
--- /dev/null
@@ -0,0 +1,83 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// TODO: Include complex types
+
+for (number : REAL_SCALARS)
+{
+  namespace Differentiation
+  \{
+  namespace AD
+  \{
+
+    // -------------------------- Types ----------------------
+
+    template
+    struct Types<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad_dfad>::ad_type>;
+
+    template
+    struct Types<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad>::ad_type>;
+
+    template
+    struct Types<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad>::ad_type>;
+
+    template
+    struct Types<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad_dfad>::ad_type>;
+    
+    // -------------------------- Numbers ----------------------
+
+    template
+    struct Numbers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad_dfad>::ad_type>;
+
+    template
+    struct Numbers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad>::ad_type>;
+
+    template
+    struct Numbers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad>::ad_type>;
+
+    template
+    struct Numbers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad_dfad>::ad_type>;
+
+    // -------------------------- TapedDrivers ----------------------
+
+    template
+    struct TapedDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad_dfad>::ad_type,number>;
+
+    template
+    struct TapedDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad>::ad_type,number>;
+
+    template
+    struct TapedDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad>::ad_type,number>;
+
+    template
+    struct TapedDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad_dfad>::ad_type,number>;
+
+    // -------------------------- TapelessDrivers ----------------------
+
+    template
+    struct TapelessDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad_dfad>::ad_type,number>;
+
+    template
+    struct TapelessDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_dfad>::ad_type,number>;
+
+    template
+    struct TapelessDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad>::ad_type,number>;
+
+    template
+    struct TapelessDrivers<typename Differentiation::AD::NumberTraits<number,NumberTypes::sacado_rad_dfad>::ad_type,number>;
+
+  \}
+  \}
+}
index 941f6c830fc2a77bdf630e609def9e9ac32f8775..c1a29cb1f38916180a7aad211c6b584324308228 100644 (file)
@@ -377,7 +377,7 @@ namespace Differentiation
       Assert(is_registered_tape(tape_index),
              ExcMessage("Tape number not registered"));
 
-      TapedDrivers<ad_type, scalar_type>::print_tape_stats(tape_index, stream);
+      this->taped_driver.print_tape_stats(tape_index, stream);
     }
 
 
@@ -829,8 +829,8 @@ namespace Differentiation
                  ExcDimensionMismatch(this->independent_variable_values.size(),
                                       this->n_independent_variables()));
 
-          return TapedDrivers<ad_type, scalar_type>::value(
-            this->active_tape_index(), this->independent_variable_values);
+          return this->taped_driver.value(this->active_tape_index(),
+                                          this->independent_variable_values);
         }
       else
         {
@@ -840,8 +840,7 @@ namespace Differentiation
                    this->n_independent_variables(),
                  ExcInternalError());
 
-          return TapelessDrivers<ad_type, scalar_type>::value(
-            this->dependent_variables);
+          return this->tapeless_driver.value(this->dependent_variables);
         }
     }
 
@@ -891,10 +890,9 @@ namespace Differentiation
                  ExcDimensionMismatch(this->independent_variable_values.size(),
                                       this->n_independent_variables()));
 
-          TapedDrivers<ad_type, scalar_type>::gradient(
-            this->active_tape_index(),
-            this->independent_variable_values,
-            gradient);
+          this->taped_driver.gradient(this->active_tape_index(),
+                                      this->independent_variable_values,
+                                      gradient);
         }
       else
         {
@@ -904,8 +902,9 @@ namespace Differentiation
                    this->n_independent_variables(),
                  ExcInternalError());
 
-          TapelessDrivers<ad_type, scalar_type>::gradient(
-            this->independent_variables, this->dependent_variables, gradient);
+          this->tapeless_driver.gradient(this->independent_variables,
+                                         this->dependent_variables,
+                                         gradient);
         }
     }
 
@@ -961,10 +960,9 @@ namespace Differentiation
                  ExcDimensionMismatch(this->independent_variable_values.size(),
                                       this->n_independent_variables()));
 
-          TapedDrivers<ad_type, scalar_type>::hessian(
-            this->active_tape_index(),
-            this->independent_variable_values,
-            hessian);
+          this->taped_driver.hessian(this->active_tape_index(),
+                                     this->independent_variable_values,
+                                     hessian);
         }
       else
         {
@@ -973,8 +971,10 @@ namespace Differentiation
           Assert(this->independent_variables.size() ==
                    this->n_independent_variables(),
                  ExcInternalError());
-          TapelessDrivers<ad_type, scalar_type>::hessian(
-            this->independent_variables, this->dependent_variables, hessian);
+
+          this->tapeless_driver.hessian(this->independent_variables,
+                                        this->dependent_variables,
+                                        hessian);
         }
     }
 
@@ -1048,18 +1048,16 @@ namespace Differentiation
                  ExcDimensionMismatch(this->independent_variable_values.size(),
                                       this->n_independent_variables()));
 
-          TapedDrivers<ad_type, scalar_type>::values(
-            this->active_tape_index(),
-            this->n_dependent_variables(),
-            this->independent_variable_values,
-            values);
+          this->taped_driver.values(this->active_tape_index(),
+                                    this->n_dependent_variables(),
+                                    this->independent_variable_values,
+                                    values);
         }
       else
         {
           Assert(ADNumberTraits<ad_type>::is_tapeless == true,
                  ExcInternalError());
-          TapelessDrivers<ad_type, scalar_type>::values(
-            this->dependent_variables, values);
+          this->tapeless_driver.values(this->dependent_variables, values);
         }
     }
 
@@ -1106,11 +1104,10 @@ namespace Differentiation
                  ExcDimensionMismatch(this->independent_variable_values.size(),
                                       this->n_independent_variables()));
 
-          TapedDrivers<ad_type, scalar_type>::jacobian(
-            this->active_tape_index(),
-            this->n_dependent_variables(),
-            this->independent_variable_values,
-            jacobian);
+          this->taped_driver.jacobian(this->active_tape_index(),
+                                      this->n_dependent_variables(),
+                                      this->independent_variable_values,
+                                      jacobian);
         }
       else
         {
@@ -1119,8 +1116,10 @@ namespace Differentiation
           Assert(this->independent_variables.size() ==
                    this->n_independent_variables(),
                  ExcInternalError());
-          TapelessDrivers<ad_type, scalar_type>::jacobian(
-            this->independent_variables, this->dependent_variables, jacobian);
+
+          this->tapeless_driver.jacobian(this->independent_variables,
+                                         this->dependent_variables,
+                                         jacobian);
         }
     }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.