NAG Library Manual, Mark 28.7
Interfaces:  FL   CL   CPP   AD
```/* nag::ad::e04ur Tangent Over Adjoint Example Program.
*   this is a copy from the e04us example, but calls into e04ur to
*   set optional parameters
*/

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

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

// Driver with the tangent over adjoint calls.
// Solve the sum-of-squares minimization problem
//    F(x) = .5 sum_{i=0}^{44} (y_i - f_i(x))^2
// for
//    f_i(x) = x_1 + (0.49 - x_1) * e^{-x_2 * (a_i - 8)}
//
// and a_i is given via ruser parameter.
//
// Also, computes the sum of all Hessian entries d2x/druser2.
void driver(std::vector<double> const &ruser,
std::vector<double> &      x,
double &                   d2xdruser2);

int main()
{

std::cout << " nag::ad::e04ur Tangent Over Adjoint Example Program Results"
<< std::endl;

// (Initial) Solution
std::vector<double> x{0.4, 0.0};

// Parameters a
std::vector<double> ruser{
8.0,  8.0,  10.0, 10.0, 10.0, 10.0, 12.0, 12.0, 12.0, 12.0, 14.0,
14.0, 14.0, 16.0, 16.0, 16.0, 18.0, 18.0, 20.0, 20.0, 20.0, 22.0,
22.0, 22.0, 24.0, 24.0, 24.0, 26.0, 26.0, 26.0, 28.0, 28.0, 30.0,
30.0, 30.0, 32.0, 32.0, 34.0, 36.0, 36.0, 38.0, 38.0, 40.0, 42.0};

// Sum of Hessian entries dx/druser
double d2xdruser2;

driver(ruser, x, d2xdruser2);

// Print outputs
std::cout << "\n Derivatives calculated: First order tangent over adjoints\n";
std::cout << " Computational mode    : algorithmic\n";

// Print derivatives
std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(12);
std::cout << "\n Solution:\n";
std::cout << " x = [" << x[0] << ", " << x[1] << "]" << std::endl;
std::cout << "\n Sum of all Hessian entries:\n";
std::cout << " sum_ijk [d2x/druser2]_ijk = " << d2xdruser2 << std::endl;

return 0;
}

// Driver with the tangent over adjoint calls.
// Solve the sum-of-squares minimization problem
//    F(x) = .5 sum_{i=0}^{44} (y_i - f_i(x))^2
// for
//    f_i(x) = x_1 + (0.49 - x_1) * e^{-x_2 * (a_i - 8)}
//
// and a_i is given via ruser parameter.
//
// Also, computes the sum of all Hessian entries d2x/druser2.
void driver(std::vector<double> const &ruserv,
std::vector<double> &      xv,
double &                   d2xdruser2)
{
using mode = dco::ga1s<dco::gt1s<double>::type>;
using T    = mode::type;

// Create AD tape
mode::global_tape = mode::tape_t::create();

std::vector<T> x(begin(xv), end(xv));
std::vector<T> ruser(begin(ruserv), end(ruserv));

dco::derivative(dco::value(ruser)) = std::vector<double>(ruser.size(), 1.0);
mode::global_tape->register_variable(ruser);

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

// Get output
xv = dco::passive_value(x);

// Set up evaluation of derivatives
mode::global_tape->register_output_variable(x);
dco::value(dco::derivative(x)) = std::vector<double>(x.size(), 1.0);

// Get sum of Hessian elements.
d2xdruser2 = 0;
for (T const &r : ruser)
{
d2xdruser2 += dco::derivative(dco::derivative(r));
}

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

template <typename T> void func(std::vector<T> &ruser, std::vector<T> &x)
{
// Problem sizes
Integer m = 44, n = 2, nclin = 1, ncnln = 1, lb = n + nclin + ncnln,
liwork = 3 * n + nclin + 2 * ncnln, lda = nclin, sda = n,
ldcj = ncnln, ldfj = m, ldr = n, lwork = 20 * n + m * (n + 3);

lwork += 2 * n * n + 11 * nclin;
lwork += n * nclin + (2 * n + 21) * ncnln;

std::vector<T> c(ncnln), cjac(ncnln * n), f(m), fjac(m * n), clamda(lb),
r(ldr * n), rwsav(475), work(lwork);
std::vector<Integer> iwsav(610), iwork(liwork), istate(lb);
std::vector<logical> lwsav(120);

// Problem data
std::vector<T> a{1.0, 1.0};
std::vector<T> y{0.49, 0.49, 0.48, 0.47, 0.48, 0.47, 0.46, 0.46, 0.45,
0.43, 0.45, 0.43, 0.43, 0.44, 0.43, 0.43, 0.46, 0.45,
0.42, 0.42, 0.43, 0.41, 0.41, 0.40, 0.42, 0.40, 0.40,
0.41, 0.40, 0.41, 0.41, 0.40, 0.40, 0.40, 0.38, 0.41,
0.40, 0.40, 0.41, 0.38, 0.40, 0.40, 0.39, 0.39};
std::vector<T> bl{0.4, -4.0, 1.0, 0.0};
std::vector<T> bu{1.0E+25, 1.0E+25, 1.0E+25, 1.0E+25};

// Initialize sav arrays
Integer ifail = 0;
char    cwsav;
nag::ad::e04wb("E04USA", &cwsav, 1, lwsav.data(), lwsav.size(), iwsav.data(),
iwsav.size(), rwsav.data(), rwsav.size(), ifail);

// Create AD configuration data object
ifail = 0;

Integer inform;
nag::ad::e04ur("Print level = 10", lwsav.data(), iwsav.data(), rwsav.data(),
inform);

Integer &          mode,
const Integer &    m,
const Integer &    n,
const Integer &    ldfj,
const Integer &    needfi,
const T            x[],
T                  f[],
T                  fjac[],
const Integer &    nstate)
{
T x1{x[0]}, x2{x[1]};

#define FJAC(i, j) fjac[(i) + (j)*m]

if (mode == 0 && needfi > 0)
{
// Update only one element f(i)
Integer i  = needfi - 1;
T       ai = ruser[i] - 8.0;
f[i]       = x1 + (0.49 - x1) * exp(-x2 * ai);
}
else
{
for (int i = 0; i < m; ++i)
{
T ai = ruser[i] - 8.0;
if (mode == 0 || mode == 2)
{
// Update the whole of f(i)
f[i] = x1 + (0.49 - x1) * exp(-x2 * ai);
}
if (mode == 1 || mode == 2)
{
// Update df/dx1
FJAC(i, 0) = 1.0 - exp(-x2 * ai);
// Update df/dx2
FJAC(i, 1) = -(0.49 - x1) * ai * exp(-x2 * ai);
}
}
}
#undef FJAC
};

Integer &          mode,
const Integer &    ncnln,
const Integer &    n,
const Integer &    ldcj,
const Integer      needc[],
const T            x[],
T                  c[],
T                  cjac[],
const Integer &    nstate)
{
// calling confun for the first time
if (nstate == 1)
{
for (int i = 0; i < ncnln * n; ++i)
{
cjac[i] = 0.0;
}
}
// evaluate c and cjac (depending on mode) at point x
if (needc[0] > 0)
{
// c required
if (mode == 0 || mode == 2)
{
c[0] = -0.09 - x[0] * x[1] + 0.49 * x[1];
}
// cjac required
if (mode == 1 || mode == 2)
{
cjac[0]     = -x[1];
cjac[ncnln] = -x[0] + 0.49;
}
}
};

// Solve the problem
Integer iter;
T       objf;
ifail = 0;
nag::ad::e04us(ad_handle, m, n, nclin, ncnln, lda, ldcj, ldfj, ldr, a.data(),
bl.data(), bu.data(), y.data(), confun, objfun, iter,
istate.data(), c.data(), cjac.data(), f.data(), fjac.data(),
clamda.data(), objf, r.data(), x.data(), iwork.data(), liwork,
work.data(), lwork, lwsav.data(), iwsav.data(), rwsav.data(), ifail);
}
```