Tutorials

Tutorials

Each example is compiled and checked by the documentation build.

If you want to...	Start with...
Differentiate a scalar objective	1. Record and evaluate a gradient
Differentiate a vector-valued function	2. Compute a Jacobian
Compute derivatives beyond first order	3. Use higher-order scalar sweeps
Avoid recording a tape	4. Use tapeless forward mode
Evaluate one tape from several threads	5. Evaluate a tape concurrently
Work with piecewise-smooth functions	6. Build an abs-normal form
Reduce memory for repeated time steps	7. Add checkpointing
Differentiate code outside the tape interpreter	8. Add an external differentiated function

1. Record and evaluate a gradient

Goal: Compute the gradient of \(f(x_0,x_1)=x_0x_1+\sin(x_0)\) at \((2,3)\).

Tape-based ADOL-C separates a computation into two phases:

1. Record: trace_on() starts a tape. The <<= operators mark independent variables, ordinary C++ expressions record operations, and >>= marks a dependent result. trace_off() completes the tape.
2. Evaluate: gradient() interprets that tape at the requested point. The same tape can be evaluated repeatedly without recording it again, provided the recorded control flow still represents the computation.

createNewTape() supplies a process-wide unique tape identifier, avoiding hard-coded tag collisions when several components use ADOL-C.

  const short tapeId = createNewTape();
  const std::array<double, 2> point{2.0, 3.0};
 
  trace_on(tapeId);
  {
    std::array<adouble, 2> x;
    double value = 0.0;
    x[0] <<= point[0];
    x[1] <<= point[1];
    adouble y = x[0] * x[1] + sin(x[0]);
    y >>= value;
  }
  trace_off();
 
  std::array<double, 2> gradientValue{};
  const int status = gradient(tapeId, static_cast<int>(point.size()),
                              point.data(), gradientValue.data());

For this function, \(\nabla f=(x_1+\cos(x_0),x_0)\). The executable checks both entries against that expression before printing:

gradient = [2.583853, 2.000000]

2. Compute a Jacobian

Goal: Differentiate the vector function \(F(x_0,x_1)=(x_0x_1,\sin(x_0)+x_1^2)\).

Recording multiple independent and dependent variables follows the same pattern as the scalar case. Use the array overloads of <<= and >>= to mark all elements in order.

jacobian() writes an \(m\times n\) matrix: one row for each dependent and one column for each independent. Its C-compatible interface accepts double**, so this example stores the values in nested std::array objects and passes row pointers.

  constexpr int inputs = 2;
  constexpr int outputs = 2;
  const short tapeId = createNewTape();
  std::array<double, inputs> point{2.0, 3.0};
 
  trace_on(tapeId);
  {
    std::array<adouble, inputs> x;
    std::array<double, outputs> values{};
    x <<= point;
    std::array<adouble, outputs> y{x[0] * x[1], sin(x[0]) + x[1] * x[1]};
    y >>= values;
  }
  trace_off();
 
  std::array<std::array<double, inputs>, outputs> storage{};
  std::array<double *, outputs> rows{storage[0].data(), storage[1].data()};
  const int status =
      jacobian(tapeId, outputs, inputs, point.data(), rows.data());

The test compares all four entries with

\[ J_F(x)= \left( \begin{array}{cc} x_1 & x_0 \\ \cos(x_0) & 2x_1 \end{array} \right). \]

Jacobian row 0 = [3.000000, 2.000000]
Jacobian row 1 = [-0.416147, 6.000000]

3. Use higher-order scalar sweeps

Goal: Compute the first three derivatives of \(f(x)=x^4\) at \(x=2\) with both hos_forward() and hos_reverse().

Higher-order forward mode propagates a Taylor polynomial through the tape. For degree \(d\), each independent row contains \(d\) coefficients. The input below represents

\[ x(t)=2+1t+0t^2+0t^3. \]

hos_forward() returns Taylor coefficients of \(f(x(t))\). For the input \(x(t)=x+t\) used here, \(Y_{k-1}=f^{(k)}(x)/k!\). Multiply by \(k!\) to recover the derivatives. For other input polynomials, the coefficients describe derivatives of the composition \(f(x(t))\).

A reverse sweep of degree \(d\) requires a preceding forward sweep with keep = d + 1. Here, keep = 3 prepares a reverse sweep of degree two. With \(x(t)=x+t\) and output weight one, the reverse coefficients are \(Z_k=f^{(k+1)}(x)/k!\).

  const short tapeId = createNewTape();
 
  trace_on(tapeId);
  {
    adouble x;
    double value = 0.0;
    x <<= 2.0;
    adouble y = x * x * x * x;
    y >>= value;
  }
  trace_off();
 
  constexpr int degree = 3;
  const std::array<double, 1> point{2.0};
 
  // X stores the coefficients of x(t) = 2 + t.
  const std::array<double, degree> inputTaylor{1.0, 0.0, 0.0};
  const double *inputRows[]{inputTaylor.data()};
  std::array<double, 1> value{};
  std::array<double, degree> outputTaylor{};
  double *outputRows[]{outputTaylor.data()};
 
  // keep=degree saves the Taylor data needed by the following reverse sweep.
  const int forwardStatus =
      hos_forward(tapeId, 1, 1, degree, degree, point.data(), inputRows,
                  value.data(), outputRows);
 
  // outputTaylor[k-1] = f^(k)(2) / k! for x(t) = 2 + t.
  const std::array<double, degree> forwardDerivatives{
      outputTaylor[0], 2.0 * outputTaylor[1], 6.0 * outputTaylor[2]};
 
  const std::array<double, 1> weight{1.0};
  constexpr int reverseDegree = degree - 1;
  std::array<double, reverseDegree + 1> gradientTaylor{};
  double *gradientRows[]{gradientTaylor.data()};
  const int reverseStatus =
      hos_reverse(tapeId, 1, 1, reverseDegree, weight.data(), gradientRows);
 
  // Z[k] is f^(k+1)(x) / k! for this scalar weighted output.
  const std::array<double, degree> reverseDerivatives{
      gradientTaylor[0], gradientTaylor[1], 2.0 * gradientTaylor[2]};

The executable checks both paths against \(f'(2)=32\), \(f''(2)=48\), and \(f'''(2)=48\):

hos_forward derivatives = [32.000000, 48.000000, 48.000000]
hos_reverse derivatives = [32.000000, 48.000000, 48.000000]

4. Use tapeless forward mode

Goal: Compute a directional derivative while evaluating the primal function, without recording a reusable tape.

Tapeless mode uses adtl::adouble. adtl::setNumDir() first selects the number of derivative directions carried by every active value. setADValue() seeds the direction, and each overloaded operation immediately propagates it.

Set the direction count before constructing active values and keep it fixed until they are destroyed. Use taped mode when reverse sweeps are required.

  adtl::setNumDir(1);
 
  adtl::adouble x = 2.0;
  const double direction = 1.0;
  x.setADValue(&direction);
 
  const adtl::adouble y = x * x + sin(x);
  const double value = y.getValue();
  const double derivative = *y.getADValue();

With direction \(\dot{x}=1\), the derivative of \(x^2+\sin(x)\) is \(2x+\cos(x)\). The executable validates both the primal value and derivative:

value = 4.909297
derivative = 3.583853

5. Evaluate a tape concurrently

Goal: Reuse one recorded scalar tape at two points on different threads.

Tape evaluation is exclusive by default and reuses owned buffers. After recording, ValueTape::setSharedMode() permits concurrent no-keep sweeps. Each task below uses fos_forward() with direction one, so it obtains the function value and first derivative independently.

Sweeps that save Taylor values require exclusive access. External differentiated functions are unsupported in shared no-keep sweeps. This example uses low-level sweeps; shared mode does not make all high-level drivers concurrent. Drivers that cache tape-owned scratch data need separate tapes or application-level synchronization.

  const short tapeId = createNewTape();
 
  trace_on(tapeId);
  {
    adouble x;
    double value = 0.0;
    x <<= 0.0;
    adouble y = x * x + sin(x);
    y >>= value;
  }
  trace_off();
 
  findTape(tapeId).setSharedMode();
 
  auto derivativeAt = [tapeId](double point) {
    double direction = 1.0;
    double value = 0.0;
    double derivative = 0.0;
    const int status =
        fos_forward(tapeId, 1, 1, 0, &point, &direction, &value, &derivative);
    return std::tuple{status, value, derivative};
  };
 
  auto first = std::async(std::launch::async, derivativeAt, 2.0);
  auto second = std::async(std::launch::async, derivativeAt, 3.0);
  const auto [firstStatus, valueAt2, derivativeAt2] = first.get();
  const auto [secondStatus, valueAt3, derivativeAt3] = second.get();

Both tasks are joined before their results are checked against \(2x+\cos(x)\):

f'(2) = 3.583853
f'(3) = 5.010008

6. Build an abs-normal form

Goal: Extract an abs-normal representation of the piecewise-smooth function \(y=x_0+x_1-|x_0|-|x_1|\).

Record the absolute-value operations, then query the switch count. enableMinMaxUsingAbs() also represents min/max operations as absolute-value switches; this example uses fabs() directly.

void taping(ADProblem &problem) {
  findTape(problem.tapeId).enableMinMaxUsingAbs();
  trace_on(problem.tapeId);
 
  {
    std::vector<adouble> ax(ADProblem::dimIn);
    std::vector<adouble> ay(ADProblem::dimOut);
 
    for (size_t i = 0; i < ADProblem::dimIn; ++i)
      ax[i] <<= problem.x[i];
 
    // Record the switches in a fixed order on every compiler.
    adouble abs0 = fabs(ax[0]);
    adouble abs1 = fabs(ax[1]);
    ay[0] = ax[0] + ax[1] - abs0 - abs1;
    ay[0] >>= problem.y[0];
  }
  trace_off();
 
  problem.numSwitches = get_num_switches(problem.tapeId);
  std::cout << "s = " << problem.numSwitches << "\n";
}

ADOLC::AbsNormalForm owns its matrix storage. ADOLC::AbsNormalForm::fromTape() obtains the dimensions from the tape, and ADOLC::abs_normal() fills the \(L\), \(Z\), \(Y\), and \(J\) blocks.

  ADOLC::AbsNormalForm anf = ADOLC::AbsNormalForm::fromTape(problem.tapeId);
 
  const int rc = ADOLC::abs_normal(problem.tapeId, problem.x, anf);
 
  std::cout << "rc = " << rc << "\n";
 
  printMatrix("L (s x s):", anf.L.data(), anf.shape.s, anf.shape.s);
  printMatrix("Z (s x n):", anf.Z.data(), anf.shape.s, anf.shape.n);
  printMatrix("Y (m x n):", anf.Y.data(), anf.shape.m, anf.shape.n);
  printMatrix("J (m x s):", anf.J.data(), anf.shape.m, anf.shape.s);

For this two-switch example, the executable verifies the shape and every matrix coefficient:

s = 2
rc = 0
L (s x s): 
0 0 
0 0 
Z (s x n): 
1 0 
0 1 
Y (m x n): 
1 1 
J (m x s): 
-1 -1 
validation passed

7. Add checkpointing to repeated time steps

Goal: Differentiate 100 Euler steps without keeping the complete unrolled time-stepping tape.

Checkpointing trades recomputation for memory. The step function is available for both adouble and double; ADOLC::CP::Context describes how many steps to execute, how many checkpoints may be retained, and where the state enters and leaves the checkpointed region.

The state variables must occupy consecutive tape locations. Call ensureContiguousLocations() before constructing them, as shown below.

  trace_on(tapeIdPart);
  {
    // ensure that the adoubles stored in y occupy consecutive locations
    currentTape().ensureContiguousLocations(dim);
    std::array<adouble, dim> y;
 
    std::array<adouble, dim> con;
 
    for (size_t i = 0; i < con.size(); ++i) {
      con[i] <<= conp[i];
      y[i] = con[i];
    }
 
    // Define the active variant of the time-step function.
    ADOLC::CP::Context cpc(tapeIdPart, tapeIdCheck, euler_step_act<adouble>);
 
    // Provide the passive variant of the time-step function.
    cpc.setDoubleFct(euler_step_act<double>);
 
    cpc.setNumberOfSteps(steps);
 
    cpc.setNumberOfCheckpoints(num_cpts);
 
    cpc.setDimensionXY(dim);
    cpc.setInput(y.data());
    cpc.setOutput(y.data());
    // Reuse the recorded time-step tape when possible.
    cpc.setAlwaysRetaping(false);
 
    cpc.checkpointing(tapeIdPart);
 
    double f[] = {0.0};
    y[0] + y[1] >>= f[0];
  }
  trace_off();
 
  const int checkpointStatus =
      gradient(tapeIdPart, dim, conp.data(), gradCheckpoint.data());

The complete-loop tape provides a reference implementation. Both approaches must reproduce the analytic gradient \((1.01^{100},1.02^{100})\):

full taping gradient = [2.704814, 7.244646]
checkpoint gradient = [2.704814, 7.244646]
validation passed

8. Add an external differentiated function

Goal: Record each Euler step as an external operation with separate derivative callbacks.

Derive from EDFobject to supply primal and derivative callbacks. This example implements them in two ways:

• a nested implementation delegates the callbacks to a small inner ADOL-C tape;
• a manual implementation evaluates the primal and first-order formulas directly.

The object must outlive every evaluation of the outer tape that refers to it. Implement each callback requested by the application and check its result against full taping.

For the Euler step, the scalar forward callback first evaluates the primal function and then applies the Jacobian to the incoming direction:

  int fos_forward(short tapeId, int m, int n, int, double *x, double *X,
                  double *y, double *Y) override {
    const int rc = eulerStep(tapeId, m, n, x, y);
    if (rc != 0)
      return rc;
 
    Y[0] = eulerScale0() * X[0];
    Y[1] = eulerScale1() * X[1];
    return 0;
  }

The scalar reverse callback applies the transposed Jacobian to the incoming adjoints. It adds to z because operations visited earlier in the reverse sweep may already have contributed:

  int fos_reverse(short, int m, int n, double *u, double *z, double *,
                  double *) override {
    if (m != 2 || n != 2)
      return -1;
 
    // Accumulate: z may already contain adjoints from later operations.
    z[0] += eulerScale0() * u[0];
    z[1] += eulerScale1() * u[1];
    return 0;
  }

The outer tape can then use either the nested or manual implementation through the same object interface:

int main() {
  const PassiveState controls{1.0, 1.0};
  PassiveState gradFull{0.0, 0.0};
  PassiveState gradNested{0.0, 0.0};
  PassiveState gradManual{0.0, 0.0};
 
  const short fullTapeId = createNewTape();
  const short nestedOuterTapeId = createNewTape();
  const short nestedInnerTapeId = createNewTape();
  const short manualOuterTapeId = createNewTape();
 
  traceFullTape(fullTapeId, controls);
  gradient(fullTapeId, 2, controls.data(), gradFull.data());
  printGradient("full taping", gradFull);
 
  traceEulerInnerTape(nestedInnerTapeId, controls);
  std::unique_ptr<EulerStepNestedEdf> nestedEdf;
  trace_on(nestedOuterTapeId);
  nestedEdf = std::make_unique<EulerStepNestedEdf>(nestedOuterTapeId,
                                                   nestedInnerTapeId);
  traceWithObject(*nestedEdf, controls);
  trace_off();
  gradient(nestedOuterTapeId, 2, controls.data(), gradNested.data());
  printGradient("taping with EDFobject (nested ADOL-C)", gradNested);
 
  std::unique_ptr<EulerStepManualEdf> manualEdf;
  trace_on(manualOuterTapeId);
  manualEdf = std::make_unique<EulerStepManualEdf>(manualOuterTapeId);
  traceWithObject(*manualEdf, controls);
  trace_off();
  gradient(manualOuterTapeId, 2, controls.data(), gradManual.data());
  printGradient("taping with EDFobject (manual callbacks)", gradManual);

The executable compares full taping, nested ADOL-C, and manual callbacks with the same analytic time-step gradient:

full taping:
 gradient=( 2.704814, 7.244646)
 
taping with EDFobject (nested ADOL-C):
 gradient=( 2.704814, 7.244646)
 
taping with EDFobject (manual callbacks):
 gradient=( 2.704814, 7.244646)
 
validation passed

Build the examples

Configure with Doxygen and Graphviz installed, then run:

cmake -S . -B build-docs -DCMAKE_BUILD_TYPE=Release
cmake --build build-docs --target docs -j

CMake runs each example in a separate working directory and includes its output only after its result checks pass. Add new examples to ADOL-C/doc/doxygen/CMakeLists.txt with a snippet region, expected results, a nonzero exit status on failure, and stable output.