NAG Library Manual, Mark 28.5
```/* nag::ad::f07ca Tangent over Adjoint Example Program.
*/

#include <dco.hpp>
#include <iostream>

// Function which calls NAG AD Library routines.
template <typename T>
void func(std::vector<T> &l,
std::vector<T> &d,
std::vector<T> &u,
std::vector<T> &x);

// Driver with the adjoint calls.
// Computes the solution to a system of linear equations Ax=b where A is a
// triagonal matrix of size n. Matrix A is stored in arrays l(n-1), d(n), u(n-1)
// that store the lower, the main and the upper diagonals. Also, computes the
// sum of the Hessian elements of output x w.r.t. all inputs l, d, u, b.
void driver(const std::vector<double> &lv,
const std::vector<double> &dv,
const std::vector<double> &uv,
const std::vector<double> &bv,
std::vector<double> &      xv,
double &                   dx2_dall2);

int main()
{

// Problem dimension
Integer n = 5;
// Matrix A stored in diagonals
std::vector<double> lv = {3.4, 3.6, 7.0, -6.0};
std::vector<double> dv = {3.0, 2.3, -5.0, -0.9, 7.1};
std::vector<double> uv = {2.1, -1.0, 1.9, 8.0};
// Right-hand-side vector b
std::vector<double> bv = {2.7, -0.5, 2.6, 0.6, 2.7};
// Computed solution to the system Ax=b
std::vector<double> xv(n);

double dx2_dall2;
// Call driver
driver(lv, dv, uv, bv, xv, dx2_dall2);

std::cout << "\n Solution point = ";
for (int i = 0; i < n; i++)
{
std::cout.width(5);
std::cout << xv[i];
}
std::cout << std::endl;

std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(12);
std::cout << "\n Derivatives calculated: Second order adjoints\n";
std::cout << " Computational mode    : algorithmic\n\n";

std::cout
<< "\n Sum of all Hessian elements of solution x w.r.t. l,d,u and b:\n";
std::cout << " sum_ij [d2x/dall2]_ij = " << dx2_dall2 << std::endl;

return 0;
}

// Driver with the adjoint calls.
// Computes the solution to a system of linear equations Ax=b where A is a
// triagonal matrix of size n. Matrix A is stored in arrays l(n-1), d(n), u(n-1)
// that store the lower, the main and the upper diagonals. Also, computes the
// sum of the Hessian elements of output x w.r.t. all inputs l, d, u, b.
void driver(const std::vector<double> &lv,
const std::vector<double> &dv,
const std::vector<double> &uv,
const std::vector<double> &bv,
std::vector<double> &      xv,
double &                   dx2_dall2)
{
using mode = dco::ga1s<dco::gt1s<double>::type>;
using T    = mode::type;

mode::global_tape = mode::tape_t::create();

Integer n  = xv.size();
Integer n1 = n - 1;
// Stores the lower diagonal of A
std::vector<T> l(n1), l1(n1);
dco::passive_value(l)          = lv;
dco::derivative(dco::value(l)) = std::vector<double>(n1, 1.0);
mode::global_tape->register_variable(l);
l1 = l;
// Stores the main diagonal of A
std::vector<T> d(n), d1(n);
dco::passive_value(d)          = dv;
dco::derivative(dco::value(d)) = std::vector<double>(n, 1.0);
mode::global_tape->register_variable(d);
d1 = d;
// Stores the upper diagonal of A
std::vector<T> u(n1), u1(n1);
dco::passive_value(u)          = uv;
dco::derivative(dco::value(u)) = std::vector<double>(n1, 1.0);
mode::global_tape->register_variable(u);
u1 = u;
// Stores right-hand-side vector b, variable to differentiate w.r.t.
std::vector<T> b(n);
dco::passive_value(b)          = bv;
dco::derivative(dco::value(b)) = std::vector<double>(n, 1.0);
mode::global_tape->register_variable(b);

// Variable to differentiate
std::vector<T> x(n);

// nag::ad::f07ca modifies rhs b and returns solution into the same array.
// Since we seek dx/db, and b is overwritten, we must save these input nodes
// in the DAG.  Hence we save b and overwrite a copy of it.
x = b;
// Call the NAG AD Lib functions
func(l1, d1, u1, x);

// Solution point
xv = dco::passive_value(x);

mode::global_tape->register_output_variable(x);
dco::derivative(x) = std::vector<dco::gt1s<double>::type>(n, 1.0);

dx2_dall2 = 0;
// Get sum of Hessian elements of solution x w.r.t. d and b
for (int i = 0; i < n; i++)
{
dx2_dall2 += dco::derivative(dco::derivative(d[i]));
dx2_dall2 += dco::derivative(dco::derivative(b[i]));
}
// Get sum of Hessian elements of solution x w.r.t. l and u
for (int i = 0; i < n1; i++)
{
dx2_dall2 += dco::derivative(dco::derivative(l[i]));
dx2_dall2 += dco::derivative(dco::derivative(u[i]));
}

// Remove tape
mode::tape_t::remove(mode::global_tape);
}

// Function which calls NAG AD Library routines.
template <typename T>
void func(std::vector<T> &l,
std::vector<T> &d,
std::vector<T> &u,
std::vector<T> &x)
{
Integer n = x.size(), nrhs = 1;

// Create AD configuration data object
Integer           ifail = 0;