```/* E04GB_A1W_F C++ Header Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
* Mark 26.2, 2017.
*/

#include <nag.h>
#include <dco.hpp>
#include <stdio.h>
#include <math.h>
#include <nag_stdlib.h>
#include <nagx02.h>
#include <iostream>
using namespace std;

extern "C"
{
static void NAG_CALL lsqlin(Integer &selct);
static void NAG_CALL lsqgrd(const Integer &m,
const Integer &n,
const Integer &ldfjac,
Integer &iflag,
const Integer &m,
const Integer &n,
const Integer &ldfjac,
Integer iuser[],
const Integer &m,
const Integer &n,
const Integer &ldfjac,
const Integer &niter,
const Integer &nf,
Integer iuser[],
}

int main(void)
{
// Scalars
int               exit_status = 0;
const Integer     m = 15, ldfjac = 15, n = 3, nt = 3, ldv = 3;

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

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

// Skip first line of data file
string mystr;
getline (cin, mystr);

Integer mode;
cin >> mode;

// Set the mode
ifail = 0;

// Read problem parameters and register for differentiation

double            yr;
for (int i=0; i<m; i++) {
cin >> yr;
ruser[i] = yr;
for (int j=1; j<=nt; j++) {
cin >> yr;
Integer k = j*m + i;
ruser[k] = yr;
}
}

Integer           iprint, maxcal, niter, nf;
Integer           iuser[1];
nagad_a1w_w_rtype x[n], fvec[m], fjac[m*n], g[n], s[n], v[ldv*n];

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

// Starting points
x[0] = 0.5;
for (int i=1; i<n; i++) {
x[i] = x[i-1] + 0.5;
}

ifail = -1;
stepmx,x,fsumsq,fvec,fjac,ldfjac,s,v,ldv,niter,nf,
iuser,ruser,ifail);

// Primal results
cout.setf(ios::scientific,ios::floatfield);
if (ifail==0 || ifail>1) {
cout.precision(4);
cout << "\n Sum of squares = ";
cout << "\n Solution point = ";
for (int i=0; i<n; i++) {
}
cout << endl;

lsqgrd(m,n,fvec,fjac,m,g);

cout.precision(3);
cout << " Estim gradient = ";
for (int i=0; i<n; i++) {
}

cout.precision(1);
cout << "\n Residuals :\n";
for (int i=0; i<m; i++) {
if (i%3==2) {
cout << endl;
}
}
}

cout << "\n Derivatives calculated: First order adjoints\n";
cout << " Computational mode    : algorithmic\n\n";
} else {
cout << " Computational mode    : symbolic\n\n";
}
cout << " Derivatives:\n\n";

// Setup evaluation of derivatives of fsumsq via adjoints.
double inc = 1.0;
ifail = 0;

// Get derivatives of fsumsq
cout << "  sum of squares w.r.t y:\n";
for (int i=0; i<m; i++) {
cout.width(12); cout << d;
if (i%3==2) {
cout << endl;
}
}

// Get derivatives of solution points w.r.t. y
for (int j=0; j<n; j++) {
ifail = 0;
cout << "\n dx(";
cout.width(1); cout << j+1;
cout << ")/dy :\n";
for (int i=0; i<m; i++) {
cout.width(12); cout << d;
if (i%3==2) {
cout << endl;
}
}
}

// Remove computational data object and tape

return exit_status;
}

// dco/c++ used to perform AD of callbacks

static void NAG_CALL lsqgrd(const Integer &m,
const Integer &n,
const Integer &ldfjac,
Integer k;
for (int i=0; i<n; i++) {
k = i*ldfjac;
for (int j=0; j<m; j++) {
g[i] = g[i] + fjac[k]*fvec[j];
k++;
}
g[i] = 2.0*g[i];
}

return;
}

static void NAG_CALL lsqlin(Integer &selct) {
selct = 2;
return;
}

Integer &iflag,
const Integer &m,
const Integer &n,
const Integer &ldfjac,
Integer iuser[],
{
for (int i=0;i<m; i++) {
fvec[i] = xc[0] + ruser[m+i]/denom - ruser[i];
if (iflag>0) {
fjac[i] = 1.0;
fjac[m+i] = ruser[m+i]*ruser[2*m+i]*denom2;
fjac[2*m+i] = ruser[m+i]*ruser[3*m+i]*denom2;
}
}
return;
}

const Integer &m,
const Integer &n,
const Integer &ldfjac,
const Integer &niter,
const Integer &nf,
Integer iuser[],
{
double fsumsq = 0.0;
for (int i=0;i<m; i++) {
fsumsq = fsumsq + f*f;
}
lsqgrd(m,n,fvec,fjac,ldfjac,g);

double gtg = 0.0;
for (int i=0;i<n; i++) {
gtg += gg*gg;
}

cout.setf(ios::scientific,ios::floatfield);
cout.precision(4);
cout << "  Itn      F evals        SUMSQ             GTG        Grade\n";
cout.width(4); cout << niter;
cout.width(11); cout << nf;
cout.width(19); cout << fsumsq;
cout.width(15); cout << gtg;
cout.width(9); cout << igrade; cout << endl;
cout << "\n       x                    g           Singular values\n";
for (int i=0;i<n;i++) {