NAG Library Manual, Mark 27.2
```/* nag::ad::e04gb Passive Example Program.
*/

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

std::stringstream filecontent("0.14  1.0 15.0  1.0\
0.18  2.0 14.0  2.0\
0.22  3.0 13.0  3.0\
0.25  4.0 12.0  4.0\
0.29  5.0 11.0  5.0\
0.32  6.0 10.0  6.0\
0.35  7.0  9.0  7.0\
0.39  8.0  8.0  8.0\
0.37  9.0  7.0  7.0\
0.58 10.0  6.0  6.0\
0.73 11.0  5.0  5.0\
0.96 12.0  4.0  4.0\
1.34 13.0  3.0  3.0\
2.10 14.0  2.0  2.0\
4.39 15.0  1.0  1.0");

// Function which calls NAG AD Library routines.
// Computes the minimum of the sum of squares of m nonlinear functions,
// the solution point and the corresponding residuals.
template <typename T>
void func(std::vector<T> &y,
std::vector<T> &t,
std::vector<T> &x,
std::vector<T> &fvec,
T &             fsumsq);

// Evaluates the residuals and their 1st derivatives.
template <typename T>
Integer &      iflag,
const Integer &m,
const Integer &n,
const T        xc[],
T              fvec[],
T              fjac[],
const Integer &ldfjac,
Integer        iuser[],
T              ruser[]);

int main(void)
{
std::cout << " nag::ad::e04gb Passive Example Program Results\n";

// Problem dimensions
const Integer       m = 15, n = 3, nt = 3;
std::vector<double> yv(m), tv(m * nt);
for (int i = 0; i < m; i++)
{
filecontent >> yv[i];
for (int j = 0; j < nt; j++)
{
Integer k = j * m + i;
filecontent >> tv[k];
}
}

std::vector<double> xv(n), fvecv(m);

// Sum of squares of the residuals at the computed point xv
double fsumsqv;

// Call the NAG AD Lib functions
func(yv, tv, xv, fvecv, fsumsqv);

// Primal results
std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(12);
std::cout << "\n Sum of squares = ";
std::cout.width(20);
std::cout << fsumsqv;
std::cout << "\n Solution point = ";
for (int i = 0; i < n; i++)
{
std::cout.width(20);
std::cout << xv[i];
}
std::cout << std::endl;

std::cout << "\n Residuals :\n";
for (int i = 0; i < m; i++)
{
std::cout.width(20);
std::cout << fvecv[i] << std::endl;
}
std::cout << std::endl;

return 0;
}

// Function which calls NAG AD Library routines.
template <typename T>
void func(std::vector<T> &y,
std::vector<T> &t,
std::vector<T> &x,
std::vector<T> &fvec,
T &             fsumsq)
{
// Problem dimensions
const Integer m = fvec.size(), n = x.size();
const Integer ldfjac = m, nt = n, ldv = n;

// All additional data accessed in the callback MUST be in ruser.
// Pack the parameters [y, t1, t2, t3] into the columns of ruser
std::vector<T> ruser(m * nt + m);
T *            _y = &ruser[0];
T *            _t = &ruser[m];

for (int i = 0; i < m; i++)
{
_y[i] = y[i];
}

for (int i = 0; i < m * nt; i++)
{
_t[i] = t[i];
}

// Initial guess of the position of the minimum
dco::passive_value(x[0]) = 0.5;
for (int i = 1; i < n; i++)
{
dco::passive_value(x[i]) = dco::passive_value(x[i - 1]) + 0.5;
}

std::vector<T> s(n), v(ldv * n), fjac(m * n);
Integer        niter, nf;

Integer iprint = -1;
Integer selct  = 2;
Integer maxcal = 200 * n;
T       eta    = 0.5;
T       xtol   = 10.0 * sqrt(X02AJC);
T       stepmx = 10.0;

// Create AD configuration data object
Integer ifail     = 0;
// Routine for computing the minimum of the sum of squares of m nonlinear
// functions.
ifail = 0;
xtol, stepmx, x.data(), fsumsq, fvec.data(), fjac.data(),
ldfjac, s.data(), v.data(), ldv, niter, nf, 0, nullptr,
ruser.size(), ruser.data(), ifail);

// Remove computational data object
ifail = 0;
}

// Evaluates the residuals and their 1st derivatives.
template <typename T>
Integer &      iflag,
const Integer &m,
const Integer &n,
const T        xc[],
T              fvec[],
T              fjac[],
const Integer &ldfjac,
Integer        iuser[],
T              ruser[])
{
const T *y  = ruser;
const T *t1 = ruser + m;
const T *t2 = ruser + 2 * m;
const T *t3 = ruser + 3 * m;
const T  x1 = xc[0], x2 = xc[1], x3 = xc[2];

for (int i = 0; i < m; i++)
{
T di = x2 * t2[i] + x3 * t3[i];
T yi = x1 + t1[i] / di;
// The values of the residuals
fvec[i] = yi - y[i];
// Evaluate the Jacobian
if (iflag > 0)
{
// dF/dx1
fjac[i] = 1.0;
// dF/dx2
fjac[i + m] = -(t1[i] * t2[i]) / (di * di);
// dF/dx3
fjac[i + 2 * m] = -(t1[i] * t3[i]) / (di * di);
}
}
}
```