```/* E04GB_P0W_F C++ Header Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
* Mark 27, 2019.
*/

#include <nag.h>
#include <nagad.h>
#include <stdio.h>
#include <math.h>
#include <nagx02.h>
#include <iostream>
#include <string>
using namespace std;

extern "C"
{
static void NAG_CALL lsqgrd(const Integer &m, const Integer &n,
const double fvec[], const double fjac[],
const Integer &ldfjac, double g[]);
static void NAG_CALL lsqfun(void* &ad_handle, Integer &iflag,
const Integer &m, const Integer &n,
const double xc[], double fvec[],
double fjac[], const Integer &ldfjac,
Integer iuser[], double ruser[]);
static void NAG_CALL lsqmon(void* &ad_handle, const Integer &m,
const Integer &n, const double xc[],
const double fvec[], const double fjac[],
const Integer &ldfjac, const double s[],
const Integer &igrade, const Integer &niter,
const Integer &nf, Integer iuser[],
double ruser[]);
}

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

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

Integer       ifail = 0;
void          *ad_handle = 0;

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

// Read computational mode
Integer       mode;
cin >> mode;

// Set the mode

ifail = 0;
// Read problem parameters and register for differentiation

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

// passive routine arguments
Integer       iprint, maxcal, niter, nf, selct;
Integer       iuser[1];
double        inc = 1.0;
double        eta, stepmx, xtol, fsumsq;
double        x[n], fvec[m], fjac[m*n], g[n], s[n], v[ldv*n];

iprint = -1;
selct = 2;
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;
}

// Call the passive routine
ifail = -1;
e04gb_p0w_f_(ad_handle,m,n,selct,lsqfun,lsqmon,iprint,maxcal,eta,xtol,
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.width(12); cout << fsumsq;
cout << "\n Solution point = ";
for (int i=0; i<n; i++) {
cout.width(12); cout << x[i];
}
cout << endl;

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

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

cout.precision(1);
cout << "\n Residuals :\n";
for (int i=0; i<m; i++) {
cout.width(12); cout << fvec[i];
if (i%3==2) {
cout << endl;
}
}
} else {
cout << "e04gb_p0w_f_ failed with ifail = "<< ifail << endl;
}

return exit_status;
}

static void NAG_CALL lsqgrd(const Integer &m, const Integer &n,
const double fvec[], const double fjac[],
const Integer &ldfjac, double g[]){
Integer k;
for (int i=0; i<n; i++) {
g[i] = 0.0;
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 lsqfun(void* &ad_handle, Integer &iflag,
const Integer &m, const Integer &n,
const double xc[], double fvec[],
double fjac[], const Integer &ldfjac,
Integer iuser[], double ruser[])
{
for (int i=0;i<m; i++) {
double d1, d2, denom;
d1 = xc[1]*ruser[2*m+i];
d2 = xc[2]*ruser[3*m+i];
denom = d1 + d2;
d1 = xc[0]- ruser[i];
d2 = ruser[m+i]/denom;
fvec[i] = d1 + d2;
if (iflag>0) {
fjac[i] = 1.0;
d1 = 1.0/denom;
d2 = -d1*d1;
d1 = ruser[m+i]*d2;
fjac[m+i] = d1*ruser[2*m+i];
fjac[2*m+i] = d1*ruser[3*m+i];
}
}
return;
}

static void NAG_CALL lsqmon(void* &ad_handle, const Integer &m,
const Integer &n, const double xc[],
const double fvec[], const double fjac[],
const Integer &ldfjac, const double s[],
const Integer &igrade, const Integer &niter,
const Integer &nf, Integer iuser[],
double ruser[])
{
double fsumsq = 0.0;
for (int i=0;i<m; i++) {
double f = fvec[i];
fsumsq = fsumsq + f*f;
}
double g[3];
lsqgrd(m,n,fvec,fjac,ldfjac,g);

double gtg = 0.0;
for (int i=0;i<n; i++) {
double gg = g[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++) {
cout.width(13); cout << xc[i];
cout.width(19); cout << g[i];
cout.width(19); cout << s[i];
cout << endl;
}
return;
}
```