]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement driver classes for internally supported AD libraries.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Sat, 18 Nov 2017 11:13:50 +0000 (12:13 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 26 Jul 2018 08:02:41 +0000 (10:02 +0200)
doc/doxygen/headers/automatic_and_symbolic_differentiation.h
doc/news/changes/minor/20171118Jean-PaulPelteret [new file with mode: 0644]
include/deal.II/differentiation/ad/ad_drivers.h [new file with mode: 0644]

index 2e8b6f701c76448bce79ab3ed9c3d7f69771dada..15ea93c26ead72f98485e341966526daa0d139c3 100644 (file)
  * A summary of the files that implement the interface to the supported auto-differentiable
  * numbers is as follows:
  *
+ * - ad_drivers.h: Provides classes that act as drivers to the interface of internally supported 
+ *   automatic differentiation libraries. These are used internally as an intermediary to the
+ *   helper classes that we provide.
  * - ad_number_types.h: Introduces an enumeration (called a type code) for the
  *   auto-differentiable number combinations that will be supported by the driver classes.
  *   The rationale behind the use of this somewhat restrictive mechanism is discussed below.
diff --git a/doc/news/changes/minor/20171118Jean-PaulPelteret b/doc/news/changes/minor/20171118Jean-PaulPelteret
new file mode 100644 (file)
index 0000000..49bc0b2
--- /dev/null
@@ -0,0 +1,5 @@
+New: Classes that act as an interface to the drivers of internally supported automatic
+differentiation libraries (ADOL-C and Sacado) have been implemented. A brief summary
+of what these classes are used for is given in the @ref auto_symb_diff module.
+<br>
+(Jean-Paul Pelteret, 2017/11/18)
diff --git a/include/deal.II/differentiation/ad/ad_drivers.h b/include/deal.II/differentiation/ad/ad_drivers.h
new file mode 100644 (file)
index 0000000..32f24d5
--- /dev/null
@@ -0,0 +1,1337 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_differentiation_ad_ad_drivers_h
+#define dealii_differentiation_ad_ad_drivers_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/differentiation/ad/ad_number_types.h>
+#include <deal.II/differentiation/ad/ad_number_traits.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/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
+  {
+
+    /**
+     * @addtogroup Exceptions
+     */
+    //@{
+
+    /**
+     * Exception denoting that a class requires some specialization
+     * in order to be used.
+     */
+    DeclExceptionMsg (ExcRequiresADNumberSpecialization,
+                      "This function is called in a class that is expected to be specialized "
+                      "for auto-differentiable numbers.");
+
+    /**
+     * Exception denoting that Adol-C is a required feature.
+     */
+    DeclExceptionMsg (ExcRequiresAdolC,
+                      "This function is only available if deal.II is compiled with ADOL-C.");
+
+    /**
+     * This exception is raised whenever the an auto-differentiable number does not
+     * support the required number of derivative operations
+     *
+     * The first parameter to the constructor is the number of derivative operations
+     * that it provides, and the second is the minimum number that are required.
+     * Both parameters are of type <tt>int</tt>.
+     */
+    DeclException2 (ExcSupportedDerivativeLevels,
+                    std::size_t, std::size_t,
+                    << "The number of derivative levels that this auto-differentiable number supports is "
+                    << arg1 << ", but it is required that it supports at least " << arg2 << " levels.");
+
+    //@}
+
+    /**
+     * A 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
+     * number types.
+     *
+     * @tparam ADNumberType A type corresponding to a supported
+     *         auto-differentiable number.
+     * @tparam ScalarType A real or complex floating point number.
+     * @tparam T An arbitrary type resulting from the application of
+     *         the SFINAE idiom to selectively specialize this class.
+     *
+     * @author Jean-Paul Pelteret, 2017
+     */
+    template<typename ADNumberType, typename ScalarType, typename T = void>
+    struct TapedDrivers
+    {
+      // This dummy class definition safely supports compilation
+      // against tapeless numbers or taped number types that have
+      // not yet been implemented.
+
+      /**
+       * @name Taping
+       */
+      //@{
+
+      /**
+       * Enable the recording mode for a given tape.
+       *
+       * @param[in] tape_index The index of the tape to be written
+       * @param[in] keep_independent_values Determines whether the numerical
+       *            values of all independent variables are recorded in the
+       *            tape buffer.
+       */
+      static void
+      enable_taping(const unsigned int &tape_index,
+                    const bool         &keep_independent_values);
+
+      /**
+       * Disable the recording mode for a given tape.
+       *
+       * @param[in] active_tape_index The index of the (currently active) tape
+       *            to be finalized and potentially written to file.
+       * @param[in] write_tapes_to_file A flag that specified whether the tape
+       *            should be written to file or kept in memory.
+       */
+      static void
+      disable_taping(const unsigned int &active_tape_index,
+                     const bool         &write_tapes_to_file);
+
+      /**
+       * Prints the statistics regarding the usage of the tapes.
+       *
+       * @param[in] stream The output stream to which the values are to be written.
+       * @param[in] tape_index The index of the tape to get the statistics of.
+       */
+      template<typename Stream>
+      static void
+      print_tape_stats(Stream             &stream,
+                       const unsigned int &tape_index);
+
+      //@}
+
+      /**
+       * @name Drivers for scalar functions (one dependent variable)
+       */
+      //@{
+
+      /**
+       * Computes the value of the scalar field.
+       *
+       * @param[in] active_tape_index The index of the tape on which the dependent
+       *            function is recorded.
+       * @param[in] n_independent_variables The number of independent variables
+       *            whose sensitivities were tracked.
+       * @param[in] independent_variables The scalar values of the independent
+       *            variables whose sensitivities were tracked.
+       *
+       * @return The scalar values of the function.
+       */
+      static ScalarType
+      value (const unsigned int            &active_tape_index,
+             const unsigned int            &n_independent_variables,
+             const std::vector<ScalarType> &independent_variables);
+
+      /**
+       * Computes the gradient of the scalar field with respect to all independent
+       * variables.
+       *
+       * @param[in] active_tape_index The index of the tape on which the dependent
+       *            function is recorded.
+       * @param[in] n_independent_variables The number of independent variables
+       *            whose sensitivities were tracked.
+       * @param[in] independent_variables The scalar values of the independent
+       *            variables whose sensitivities were tracked.
+       * @param[out] gradient The scalar values of the dependent function's gradients.
+       */
+      static void
+      gradient (const unsigned int            &active_tape_index,
+                const unsigned int            &n_independent_variables,
+                const std::vector<ScalarType> &independent_variables,
+                Vector<ScalarType>            &gradient);
+
+      /**
+       * Computes the hessian of the scalar field with respect to all independent
+       * variables.
+       *
+       * @param[in] active_tape_index The index of the tape on which the dependent
+       *            function is recorded.
+       * @param[in] n_independent_variables The number of independent variables
+       *            whose sensitivities were tracked.
+       * @param[in] independent_variables The scalar values of the independent
+       *            variables whose sensitivities were tracked.
+       * @param[out] hessian The scalar values of the dependent function's hessian.
+       */
+      static void
+      hessian (const unsigned int            &active_tape_index,
+               const unsigned int            &n_independent_variables,
+               const std::vector<ScalarType> &independent_variables,
+               FullMatrix<ScalarType>        &hessian);
+
+      //@}
+
+      /**
+       * @name Drivers for vector functions (multiple dependent variables)
+       */
+      //@{
+
+      /**
+       * Computes the values of the vector field.
+       *
+       * @param[in] active_tape_index The index of the tape on which the dependent
+       *            function is recorded.
+       * @param[in] n_dependent_variables The number of dependent variables.
+       * @param[in] n_independent_variables The number of independent variables
+       *            whose sensitivities were tracked.
+       * @param[in] independent_variables The scalar values of the independent
+       *            variables whose sensitivities were tracked.
+       * @param[out] values The scalar values of the dependent functions.
+       */
+      static void
+      values (const unsigned int            &active_tape_index,
+              const unsigned int            &n_dependent_variables,
+              const unsigned int            &n_independent_variables,
+              const std::vector<ScalarType> &independent_variables,
+              Vector<ScalarType>            &values);
+
+      /**
+       * Computes the gradient of the vector field.
+       *
+       * @param[in] active_tape_index The index of the tape on which the dependent
+       *            function is recorded.
+       * @param[in] n_dependent_variables The number of dependent variables.
+       * @param[in] n_independent_variables The number of independent variables
+       *            whose sensitivities were tracked.
+       * @param[in] independent_variables The scalar values of the independent
+       *            variables whose sensitivities were tracked.
+       * @param[out] jacobian The scalar values of the dependent function's jacobian.
+       */
+      static void
+      jacobian (const unsigned int            &active_tape_index,
+                const unsigned int            &n_dependent_variables,
+                const unsigned int            &n_independent_variables,
+                const std::vector<ScalarType> &independent_variables,
+                FullMatrix<ScalarType>        &jacobian);
+
+      //@}
+
+    };
+
+
+
+    /**
+     * A prototype driver class for tapeless auto-differentiable numbers.
+     *
+     * It is intended that this class be specialized for the valid
+     * combinations of auto-differentiable numbers and output scalar
+     * number types.
+     *
+     * @tparam ADNumberType A type corresponding to a supported
+     *         auto-differentiable number.
+     * @tparam ScalarType A real or complex floating point number.
+     * @tparam T An arbitrary type resulting from the application of
+     *         the SFINAE idiom to selectively specialize this class.
+     *
+     * @author Jean-Paul Pelteret, 2017
+     */
+    template<typename ADNumberType, typename ScalarType, typename T = void>
+    struct TapelessDrivers
+    {
+      // This dummy class definition safely supports compilation
+      // against taped numbers or tapeless number types that have
+      // not yet been implemented.
+
+      /**
+       * @name Configuration
+       */
+      //@{
+
+      /**
+       * In the event that the tapeless mode requires <i>a priori</i> knowledge
+       * of how many directional derivatives might need to be computed, this function
+       * informs the auto-differentiable library of what this number is.
+       *
+       * @param[in] n_independent_variables The number of independent variables
+       *            that will be used in the entire duration of the
+       *            simulation.
+       *
+       * @warning For Adol-C tapeless numbers, the value given to @p n_independent_variables
+       *          should be the <b>maximum</b> number of independent variables that will be
+       *          used in the entire duration of the simulation. This is important in the
+       *          context of, for example, hp-FEM and for multiple constitutive models with
+       *          a different number of fields from which a linearization must be computed.
+       */
+      static void
+      initialize (const unsigned int &n_independent_variables);
+
+      //@}
+
+      /**
+       * @name Drivers for scalar functions
+       */
+      //@{
+
+      /**
+       * Computes the value of the scalar field.
+       *
+       * @param[in] dependent_variables The dependent variables whose values are to
+       *            be extracted.
+       *
+       * @return The scalar values of the function.
+       */
+      static ScalarType
+      value (const std::vector<ADNumberType> &dependent_variables);
+
+      /**
+       * Computes the gradient of the scalar field with respect to all independent
+       * variables.
+       *
+       * @param[in] independent_variables The independent variables whose sensitivities
+       *            were tracked.
+       * @param[in] dependent_variables The (single) dependent variable whose gradients
+       *            are to be extracted.
+       * @param[out] gradient The scalar values of the dependent function's gradients.
+       */
+      static void
+      gradient (const std::vector<ADNumberType> &independent_variables,
+                const std::vector<ADNumberType> &dependent_variables,
+                Vector<ScalarType>              &gradient);
+
+      /**
+       * Computes the hessian of the scalar field with respect to all independent
+       * variables.
+       *
+       * @param[in] independent_variables The independent variables whose sensitivities
+       *            were tracked.
+       * @param[in] dependent_variables The (single) dependent variable whose hessians
+       *            are to be extracted.
+       * @param[out] hessian The scalar values of the function's hessian.
+       */
+      static void
+      hessian (const std::vector<ADNumberType> &independent_variables,
+               const std::vector<ADNumberType> &dependent_variables,
+               FullMatrix<ScalarType>          &hessian);
+
+      //@}
+
+      /**
+       * @name Drivers for vector functions
+       */
+      //@{
+
+      /**
+       * Computes the values of the vector field.
+       *
+       * @param[in] dependent_variables The dependent variables  whose hessians
+       *            are to be extracted.
+       * @param[out] values The scalar values of the dependent functions.
+       */
+      static void
+      values (const std::vector<ADNumberType> &dependent_variables,
+              Vector<ScalarType>              &values);
+
+      /**
+       * Computes the gradient of the vector field.
+       *
+       * @param[in] independent_variables The independent variables whose sensitivities
+       *            were tracked.
+       * @param[in] dependent_variables The dependent variables whose jacobian
+       *            are to be extracted.
+       * @param[out] jacobian The scalar values of the function's jacobian.
+       */
+      static void
+      jacobian (const std::vector<ADNumberType> &independent_variables,
+                const std::vector<ADNumberType> &dependent_variables,
+                FullMatrix<ScalarType>          &jacobian);
+
+      //@}
+
+    };
+
+  }
+}
+
+
+
+/* --------------------------- inline and template functions ------------------------- */
+
+
+#ifndef DOXYGEN
+
+namespace Differentiation
+{
+  namespace AD
+  {
+
+    // -------------   TapedDrivers   -------------
+
+
+    template<typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType,ScalarType,T>::enable_taping(
+      const unsigned int &,
+      const bool &)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template<typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType,ScalarType,T>::disable_taping(
+      const unsigned int &,
+      const bool &)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template<typename ADNumberType, typename ScalarType, typename T>
+    template<typename Stream>
+    void
+    TapedDrivers<ADNumberType,ScalarType,T>::print_tape_stats(
+      Stream             &stream,
+      const unsigned int &tape_index)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template<typename ADNumberType, typename ScalarType, typename T>
+    ScalarType
+    TapedDrivers<ADNumberType,ScalarType,T>::value (
+      const unsigned int &,
+      const unsigned int &,
+      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 unsigned int &,
+      const unsigned int &,
+      const std::vector<ScalarType> &,
+      Vector<ScalarType> &)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template<typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType,ScalarType,T>::hessian (
+      const unsigned int &,
+      const unsigned int &,
+      const std::vector<ScalarType> &,
+      FullMatrix<ScalarType> &)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    template<typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapedDrivers<ADNumberType,ScalarType,T>::values (
+      const unsigned int &,
+      const unsigned int &,
+      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 unsigned int &,
+      const unsigned int &,
+      const unsigned int &,
+      const std::vector<ScalarType> &,
+      FullMatrix<ScalarType> &)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    // Specialization for taped Adol-C auto-differentiable numbers.
+    //
+    // Note: In the case of Adol-C taped numbers, the associated scalar
+    // type is always expected to be a double. So we need to make a further
+    // specialization when ScalarType is a float.
+    template<typename ADNumberType>
+    struct TapedDrivers<ADNumberType,double,typename std::enable_if<
+      ADNumberTraits<ADNumberType>::type_code == NumberTypes::adolc_taped
+      >::type>
+    {
+      typedef double scalar_type;
+
+      // === Taping ===
+
+      static void
+      enable_taping(const unsigned int &tape_index,
+                    const bool         &keep)
+      {
+        trace_on(tape_index,keep);
+      }
+
+      static void
+      disable_taping(const unsigned int &active_tape_index,
+                     const bool         &write_tapes_to_file)
+      {
+#ifdef DEAL_II_WITH_ADOLC
+        if (write_tapes_to_file)
+          {
+            trace_off(active_tape_index); // Slow
+            std::vector<std::size_t> counts (STAT_SIZE);
+            ::tapestats(active_tape_index, counts.data());
+          }
+        else
+          trace_off(); // Fast(er)
+#else
+        AssertThrow(false, ExcRequiresAdolC());
+#endif
+      }
+
+      template<typename Stream>
+      static void
+      print_tape_stats(Stream             &stream,
+                       const unsigned int &tape_index)
+      {
+        // 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) interchangable without retaping: " << counts[17] << "\n"
+            << std::flush;
+      }
+
+      // === Scalar drivers ===
+
+      static scalar_type
+      value (const unsigned int             &active_tape_index,
+             const unsigned int             &n_independent_variables,
+             const std::vector<scalar_type> &independent_variables)
+      {
+
+        scalar_type *f = new scalar_type();
+
+#ifdef DEAL_II_WITH_ADOLC
+        ::function(active_tape_index,
+                   1, // Only one dependent variable
+                   n_independent_variables,
+                   const_cast<scalar_type *>(independent_variables.data()),
+                   f);
+#else
+        AssertThrow(false, ExcRequiresAdolC());
+#endif
+
+        const scalar_type value = f[0];
+
+        // Cleanup :-/
+        delete f;
+        f = nullptr;
+
+        return value;
+      }
+
+      static void
+      gradient (const unsigned int             &active_tape_index,
+                const unsigned int             &n_independent_variables,
+                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()));
+
+        scalar_type *g = new scalar_type[n_independent_variables];
+
+#ifdef DEAL_II_WITH_ADOLC
+        ::gradient(active_tape_index,
+                   n_independent_variables,
+                   const_cast<scalar_type *>(independent_variables.data()),
+                   g);
+#else
+        AssertThrow(false, ExcRequiresAdolC());
+#endif
+
+        for (unsigned int i=0; i<n_independent_variables; ++i)
+          gradient[i] = g[i];
+
+        // Cleanup :-/
+        delete[] g;
+        g = nullptr;
+      }
+
+      static void
+      hessian (const unsigned int             &active_tape_index,
+               const unsigned int             &n_independent_variables,
+               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()));
+
+        scalar_type **H = new scalar_type*[n_independent_variables];
+        for (unsigned int i=0; i<n_independent_variables; ++i)
+          H[i] = new scalar_type[i+1]; // Symmetry
+
+#ifdef DEAL_II_WITH_ADOLC
+        ::hessian(active_tape_index,
+                  n_independent_variables,
+                  const_cast<scalar_type *>(independent_variables.data()),
+                  H);
+#else
+        AssertThrow(false, ExcRequiresAdolC());
+#endif
+
+        for (unsigned int i=0; i<n_independent_variables; i++)
+          for (unsigned int j=0; j<i+1; j++)
+            {
+              hessian[i][j] = H[i][j];
+              if (i != j)
+                hessian[j][i] = H[i][j]; // Symmetry
+            }
+
+        // Cleanup :-/
+        for (unsigned int i=0; i<n_independent_variables; i++)
+          delete[] H[i];
+        delete[] H;
+        H = nullptr;
+      }
+
+      // === Vector drivers ===
+
+      static void
+      values (const unsigned int             &active_tape_index,
+              const unsigned int             &n_dependent_variables,
+              const unsigned int             &n_independent_variables,
+              const std::vector<scalar_type> &independent_variables,
+              Vector<scalar_type>            &values)
+      {
+        Assert(values.size() == n_dependent_variables,
+               ExcDimensionMismatch(values.size(),n_dependent_variables));
+
+        scalar_type *f = new scalar_type[n_dependent_variables];
+
+#ifdef DEAL_II_WITH_ADOLC
+        ::function(active_tape_index,
+                   n_dependent_variables,
+                   n_independent_variables,
+                   const_cast<scalar_type *>(independent_variables.data()),
+                   f);
+#else
+        AssertThrow(false, ExcRequiresAdolC());
+#endif
+
+        for (unsigned int i=0; i<n_dependent_variables; i++)
+          values[i] = f[i];
+
+        // Cleanup :-/
+        delete[] f;
+        f = nullptr;
+      }
+
+      static void
+      jacobian (const unsigned int             &active_tape_index,
+                const unsigned int             &n_dependent_variables,
+                const unsigned int             &n_independent_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()));
+
+        scalar_type **J = new scalar_type*[n_dependent_variables];
+        for (unsigned int i=0; i<n_dependent_variables; ++i)
+          J[i] = new scalar_type[n_independent_variables];
+
+#ifdef DEAL_II_WITH_ADOLC
+        ::jacobian(active_tape_index,
+                   n_dependent_variables,
+                   n_independent_variables,
+                   independent_variables.data(),
+                   J);
+#else
+        AssertThrow(false, ExcRequiresAdolC());
+#endif
+
+        for (unsigned int i=0; i<n_dependent_variables; i++)
+          for (unsigned int j=0; j<n_independent_variables; j++)
+            jacobian[i][j] = J[i][j];
+
+        // Cleanup :-/
+        for (unsigned int i=0; i<n_dependent_variables; i++)
+          delete[] J[i];
+        delete[] J;
+        J = nullptr;
+      }
+    };
+
+
+    // Specialization for Adol-C taped numbers. It is expected that the
+    // scalar return type for this class is a float.
+    //
+    // Note: Adol-C only has drivers for doubles, and so floats are
+    // not intrinsically supported. This wrapper struct works around
+    // the issue when necessary.
+    template<typename ADNumberType>
+    struct TapedDrivers<ADNumberType,float,typename std::enable_if<
+      ADNumberTraits<ADNumberType>::type_code == NumberTypes::adolc_taped
+      >::type>
+    {
+      typedef float scalar_type;
+
+      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;
+      }
+
+      // === Taping ===
+
+      static void
+      enable_taping(const unsigned int &tape_index,
+                    const bool         &keep)
+      {
+        TapedDrivers<ADNumberType,double>::enable_taping(tape_index,keep);
+      }
+
+      static void
+      disable_taping(const unsigned int &active_tape_index,
+                     const bool         &write_tapes_to_file)
+      {
+        TapedDrivers<ADNumberType,double>::disable_taping(active_tape_index,write_tapes_to_file);
+      }
+
+      template<typename Stream>
+      static void
+      print_tape_stats(Stream             &stream,
+                       const unsigned int &tape_index)
+      {
+        TapedDrivers<ADNumberType,double>::print_tape_stats(stream,tape_index);
+      }
+
+      // === Scalar drivers ===
+
+      static scalar_type
+      value (const unsigned int             &active_tape_index,
+             const unsigned int             &n_independent_variables,
+             const std::vector<scalar_type> &independent_variables)
+      {
+
+        return TapedDrivers<ADNumberType,double>::value(
+                 active_tape_index,
+                 n_independent_variables,
+                 vector_float_to_double(independent_variables));
+      }
+
+      static void
+      gradient (const unsigned int             &active_tape_index,
+                const unsigned int             &n_independent_variables,
+                const std::vector<scalar_type> &independent_variables,
+                Vector<scalar_type>            &gradient)
+      {
+        Vector<double> gradient_double (gradient.size());
+        TapedDrivers<ADNumberType,double>::gradient(
+          active_tape_index,
+          n_independent_variables,
+          vector_float_to_double(independent_variables),
+          gradient_double);
+        gradient = gradient_double;
+      }
+
+      static void
+      hessian (const unsigned int             &active_tape_index,
+               const unsigned int             &n_independent_variables,
+               const std::vector<scalar_type> &independent_variables,
+               FullMatrix<scalar_type>        &hessian)
+      {
+        FullMatrix<double> hessian_double (hessian.m(), hessian.n());
+        TapedDrivers<ADNumberType,double>::hessian(
+          active_tape_index,
+          n_independent_variables,
+          vector_float_to_double(independent_variables),
+          hessian_double);
+        hessian = hessian_double;
+      }
+
+      // === Vector drivers ===
+
+      static void
+      values (const unsigned int             &active_tape_index,
+              const unsigned int             &n_dependent_variables,
+              const unsigned int             &n_independent_variables,
+              const std::vector<scalar_type> &independent_variables,
+              Vector<scalar_type>            &values)
+      {
+        Vector<double> values_double (values.size());
+        TapedDrivers<ADNumberType,double>::values(
+          active_tape_index,
+          n_dependent_variables,
+          n_independent_variables,
+          vector_float_to_double(independent_variables),
+          values_double);
+        values = values_double;
+      }
+
+      static void
+      jacobian (const unsigned int             &active_tape_index,
+                const unsigned int             &n_dependent_variables,
+                const unsigned int             &n_independent_variables,
+                const std::vector<scalar_type> &independent_variables,
+                FullMatrix<scalar_type>        &jacobian)
+      {
+        FullMatrix<double> jacobian_double (jacobian.m(), jacobian.n());
+        TapedDrivers<ADNumberType,double>::jacobian(
+          active_tape_index,
+          n_dependent_variables,
+          n_independent_variables,
+          vector_float_to_double(independent_variables),
+          jacobian_double);
+        jacobian = jacobian_double;
+      }
+    };
+
+
+    // -------------   TapelessDrivers   -------------
+
+
+    template<typename ADNumberType, typename ScalarType, typename T>
+    void
+    TapelessDrivers<ADNumberType,ScalarType,T>::initialize (
+      const unsigned int &)
+    {
+      AssertThrow(false, ExcRequiresADNumberSpecialization());
+    }
+
+
+    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>
+      static typename std::enable_if<
+      !is_sacado_rad_number<ADNumberType>::value
+      >::type
+      reverse_mode_dependent_variable_activation (ADNumberType &)
+      {
+
+      }
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+
+      // 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>
+      static typename std::enable_if<
+      is_sacado_rad_number<ADNumberType>::value
+      >::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 all
+        // 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>
+      static typename std::enable_if<
+      !(is_adolc_number<ADNumberType>::value &&
+        is_tapeless_ad_number<ADNumberType>::value)
+      >::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>
+      static typename std::enable_if<
+      is_adolc_number<ADNumberType>::value &&
+      is_tapeless_ad_number<ADNumberType>::value
+      >::type
+      configure_tapeless_mode (const unsigned int n_directional_derivatives)
+      {
+        // 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 number of computable
+            // directional derivatives. If this really is necessary then its
+            // 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."));
+          }
+      }
+
+#endif
+
+    }
+
+
+    // 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>
+    struct TapelessDrivers<ADNumberType,ScalarType,typename std::enable_if<
+      ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_rad ||
+      ADNumberTraits<ADNumberType>::type_code == NumberTypes::sacado_rad_dfad
+      >::type>
+    {
+
+      // === Configuration ===
+
+      static void
+      initialize (const unsigned int &n_independent_variables)
+      {
+        internal::configure_tapeless_mode<ADNumberType>(n_independent_variables);
+      }
+
+      // === 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]);
+      }
+
+      static 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*/));
+      }
+
+      static 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++)
+          {
+            typedef typename ADNumberTraits<ADNumberType>::derivative_type 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 ===
+
+      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]);
+      }
+
+      static 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.
+        typedef typename ADNumberTraits<ADNumberType>::derivative_type AccumulationType;
+        std::vector<AccumulationType> rad_accumulation (
+          n_independent_variables,
+          dealii::internal::NumberType<AccumulationType>::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 AccumulationType 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>
+    struct 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>
+    {
+
+      // === Configuration ===
+
+      static void
+      initialize (const unsigned int &n_independent_variables)
+      {
+        internal::configure_tapeless_mode<ADNumberType>(n_independent_variables);
+      }
+
+      // === 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]);
+      }
+
+      static 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));
+      }
+
+      static 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++)
+          {
+            typedef typename ADNumberTraits<ADNumberType>::derivative_type 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 ===
+
+      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]);
+      }
+
+      static 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));
+      }
+
+    };
+
+
+  } // namespace AD
+} // namespace Differentiation
+
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif

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.