\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)
mul_level_adolc.cpp¶
View page sourceUsing Adolc with Multiple Levels of Taping: Example and Test¶
Purpose¶
In this example, we use AD< adouble> >
(level two taping),
the compute values of the function \(f : \B{R}^n \rightarrow \B{R}\) where
We then use Adolc’s adouble
(level one taping) to compute
the directional derivative
where \(v \in \B{R}^n\).
We then use double
(no taping) to compute
This is only meant as an example of multiple levels of taping. The example hes_times_dir.cpp computes the same value more efficiently by using the identity:
The example mul_level.cpp computes the same values using
AD< AD<double> >
and AD<double>
.
Memory Management¶
Adolc uses raw memory arrays that depend on the number of dependent and independent variables. The memory management utility thread_alloc is used to manage this memory allocation.
Configuration Requirement¶
This example will be compiled and tested provided include_adolc is true on the cmake command line.
Source¶
// suppress conversion warnings before other includes
# include <cppad/wno_conversion.hpp>
//
# include <adolc/adouble.h>
# include <adolc/taping.h>
# include <adolc/interfaces.h>
// adouble definitions not in Adolc distribution and
// required in order to use CppAD::AD<adouble>
# include <cppad/example/base_adolc.hpp>
# include <cppad/cppad.hpp>
namespace {
// f(x) = |x|^2 / 2 = .5 * ( x[0]^2 + ... + x[n-1]^2 )
template <class Type>
Type f(const CPPAD_TESTVECTOR(Type)& x)
{ Type sum;
sum = 0.;
size_t i = size_t(x.size());
for(i = 0; i < size_t(x.size()); i++)
sum += x[i] * x[i];
return .5 * sum;
}
}
bool mul_level_adolc(void)
{ bool ok = true; // initialize test result
using CppAD::thread_alloc; // The CppAD memory allocator
typedef adouble a1type; // for first level of taping
typedef CppAD::AD<a1type> a2type; // for second level of taping
size_t n = 5; // number independent variables
size_t j;
// 10 times machine epsilon
double eps = 10. * std::numeric_limits<double>::epsilon();
CPPAD_TESTVECTOR(double) x(n);
CPPAD_TESTVECTOR(a1type) a1x(n);
CPPAD_TESTVECTOR(a2type) a2x(n);
// Values for the independent variables while taping the function f(x)
for(j = 0; j < n; j++)
a2x[j] = double(j);
// Declare the independent variable for taping f(x)
CppAD::Independent(a2x);
// Use AD<adouble> to tape the evaluation of f(x)
CPPAD_TESTVECTOR(a2type) a2y(1);
a2y[0] = f(a2x);
// Declare a1f as the corresponding ADFun<adouble> function f(x)
// (make sure we do not run zero order forward during constructor)
CppAD::ADFun<a1type> a1f;
a1f.Dependent(a2x, a2y);
// Value of the independent variables whitle taping f'(x) * v
short tag = 0;
int keep = 1;
trace_on(tag, keep);
for(j = 0; j < n; j++)
a1x[j] <<= double(j);
// set the argument value x for computing f'(x) * v
a1f.Forward(0, a1x);
// compute f'(x) * v
CPPAD_TESTVECTOR(a1type) a1v(n);
CPPAD_TESTVECTOR(a1type) a1df(1);
for(j = 0; j < n; j++)
a1v[j] = double(n - j);
a1df = a1f.Forward(1, a1v);
// declare Adolc function corresponding to f'(x) * v
double df;
a1df[0] >>= df;
trace_off();
// compute the d/dx of f'(x) * v = f''(x) * v
size_t m = 1; // # dependent in f'(x) * v
// w = new double[capacity] where capacity >= m
size_t capacity;
double* w = thread_alloc::create_array<double>(m, capacity);
// dw = new double[capacity] where capacity >= n
double* dw = thread_alloc::create_array<double>(n, capacity);
w[0] = 1.;
fos_reverse(tag, int(m), int(n), w, dw);
for(j = 0; j < n; j++)
{ double vj = a1v[j].value();
ok &= CppAD::NearEqual(dw[j], vj, eps, eps);
}
// make memory avaialble for other use by this thread
thread_alloc::delete_array(w);
thread_alloc::delete_array(dw);
return ok;
}