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

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

std::stringstream filecontent("4 \
2.0   -1.0    0.0    0.0 \
-1.0    2.0   -1.0    0.0 \
0.0   -1.0    2.0   -1.0 \
0.0    0.0   -1.0    2.0 ");

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

// Computes the Nearest Correlation Matrix X for an input matrix G.
// Also, computes the sum of all Hessian elements of d2X/dG2.
void driver(const std::vector<double> &gv,
std::vector<double> &      xv,
double &                   d2x_dg2);

int main()
{

Integer n;
filecontent >> n;

// Input Matrix G, whose NCM we want to compute
std::vector<double> gv(n * n);
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
filecontent >> gv[i + j * n];
}
}

// Output NCM X
std::vector<double> xv(n * n);

// Sum of all Hessian elements d2X/dG2
double d2x_dg2;

driver(gv, xv, d2x_dg2);

std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(5);
std::cout << "\n Nearest Correlation Matrix:\n\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
std::cout.width(13);
std::cout << xv[i + j * n];
}
std::cout << std::endl;
}

std::cout << "\n Derivatives calculated: Second order adjoint\n";
std::cout << " Computational mode    : algorithmic\n\n";

// Print derivatives of NCM X w.r.t. input matrix G
std::cout << "\n Sum of Hessian elements of NCM X w.r.t. input matrix G:\n";
std::cout << " sum_ij [d2X/dG2]_ij = " << d2x_dg2 << std::endl;

return 0;
}

// Computes the Nearest Correlation Matrix X for an input matrix G.
// Also, computes the sum of all Hessian elements of d2X/dG2.
void driver(const std::vector<double> &gv,
std::vector<double> &      xv,
double &                   d2x_dg2)
{
using mode = dco::ga1s<dco::gt1s<double>::type>;
using T    = mode::type;
mode::global_tape = mode::tape_t::create();

std::vector<T> g(gv.size()), g1(gv.size());
dco::passive_value(g)          = gv;
dco::derivative(dco::value(g)) = std::vector<double>(gv.size(), 1.0);
mode::global_tape->register_variable(g);
// nag::ad::g02aa modifies input matrix g.  Since we seek dx/dg, and g is
// overwritten, we must save these input nodes in the DAG.  Hence we save g
// and overwrite a copy of it.
g1 = g;

// Variable to differentiate
std::vector<T> x(xv.size());

// Call the NAG AD Lib functions
func(g1, x);

// Extract the computed NCM
xv = dco::passive_value(x);

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

// Sum of the adjoints of g
d2x_dg2 = 0.0;
for (int i = 0; i < g.size(); i++)
{
d2x_dg2 += dco::derivative(dco::derivative(g[i]));
}

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

// Function which calls NAG AD routines.
template <typename T> void func(std::vector<T> &g, std::vector<T> &x)
{
Integer n   = sqrt(g.size());
Integer pdg = n;
Integer pdx = n;

// Set up method parameters
T       errtol = 1.00e-7;
Integer maxits = 200;
Integer maxit  = 10;
// Output variables
Integer iter, feval;
T       nrmgrd;
// Create AD configuration data object
Integer ifail     = 0;