NAG Library Manual, Mark 29.2
```/* G02DA_T1W_F C++ Header Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
* Mark 29.2, 2023.
*/

#include <dco.hpp>
#include <iostream>
#include <math.h>
#include <nag.h>
#include <nagx04.h>
#include <stdio.h>
#include <string>

typedef double                   DCO_BASE_TYPE;
typedef dco::gt1s<DCO_BASE_TYPE> DCO_MODE;
typedef DCO_MODE::type           DCO_TYPE;

using namespace std;

int main()
{
int               exit_status = 0;
DCO_TYPE          tol;
Integer           ifail = 0;

cout << "G02DA_T1W_F C++ Header Example Program Results\n\n";

tol = 0.000001;

// Skip heading in data file
string mystr;
getline(cin, mystr);

Integer n, m;
cin >> n;
cin >> m;

// Allocate arrays depending on m and n.
DCO_TYPE *x = 0, *y = 0, *y_in = 0;
Integer * isx = 0;
double *  dy  = 0;
x             = new DCO_TYPE[m * n];
y             = new DCO_TYPE[n];
y_in          = new DCO_TYPE[n];
dy            = new double[n];
isx           = new Integer[m];

// Create AD configuration data object
ifail = 0;

double dd;
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < m; ++j)
{
cin >> dd;
int k = i + j * n;
x[k]  = dd;
}
cin >> dd;
y_in[i] = dd;
y[i]    = y_in[i];
}

// Calculate ip
Integer ip = 0;
for (int j = 0; j < m; ++j)
{
cin >> isx[j];
if (isx[j] > 0)
{
ip++;
}
}
// Mean = 'M'
ip++;

// Allocate arrays depending on ip
DCO_TYPE *b = 0, *cov = 0, *h = 0, *p = 0, *q = 0, *res = 0, *se = 0, *wk = 0;
double *  dbdy;
Integer   lcov = (ip * ip + ip) / 2, lp = ip * (ip + 2), lq = n * (ip + 1);
Integer   lwk = ip * ip + 5 * (ip - 1);
b             = new DCO_TYPE[ip];
cov           = new DCO_TYPE[lcov];
h             = new DCO_TYPE[n];
p             = new DCO_TYPE[lp];
q             = new DCO_TYPE[lq];
res           = new DCO_TYPE[n];
se            = new DCO_TYPE[ip];
wk            = new DCO_TYPE[lwk];
dbdy          = new double[n * ip];

// Perform Regression
Integer  idf, irank;
logical  svd;
ifail = 0;
se, cov, res, h, q, n, svd, irank, p, tol, wk, ifail);

// Display results
if (svd)
{
cout << "Model is not of full rank, rank = " << irank << endl;
}
cout << "Residual sum of squares         = " << dco::value(rss);
cout << "\nDegrees of freedom              = " << idf << endl;
cout << "\nVariable   Parameter estimate   Standard error\n\n";
cout.setf(ios::scientific, ios::floatfield);
cout.precision(2);
for (int i = 0; i < ip; ++i)
{
cout.width(5);
cout << i << "       ";
cout.width(12);
cout << dco::value(b[i]) << "       ";
cout.width(12);
cout << dco::value(se[i]) << endl;
}

cout << "\n\n Derivatives calculated: First order tangents\n";
cout << " Computational mode    : algorithmic\n";
cout << "\n Derivatives:\n\n";

// Obtain derivatives
for (int i = 0; i < n; i++)
{
y[i] = y_in[i];
}
for (int i = 0; i < n; i++)
{
dco::derivative(y[i]) = 1.0;
b, se, cov, res, h, q, n, svd, irank, p, tol, wk, ifail);

for (int j = 0; j < ip; j++)
{
int k   = i + j * n;
dbdy[k] = dco::derivative(b[j]);
}
for (int j = 0; j < n; ++j)
{
y[j] = y_in[j];
}
}
cout.precision(4);
for (int i = 0; i < n; ++i)
{
cout.width(5);
cout << i << "       ";
cout.width(12);
cout << dy[i] << endl;
}

// Print matrix routine
cout << endl;

NagError fail;
INIT_FAIL(fail);
x04cac(Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, n, ip, dbdy, n,
"db/dy", 0, &fail);

ifail = 0;

delete[] b;
delete[] cov;
delete[] h;
delete[] p;
delete[] q;
delete[] res;
delete[] se;
delete[] wk;
delete[] dbdy;
delete[] x;
delete[] y;
delete[] y_in;
delete[] dy;
delete[] isx;
return exit_status;
}
```