/* E04FC_A1W_F C++ Header Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
* Mark 27, 2019.
*/
#include <nag.h>
#include <dco.hpp>
#include <nagad.h>
#include <stdio.h>
#include <math.h>
#include <nag_stdlib.h>
#include <nagx02.h>
#include <iostream>
using namespace std;
extern "C"
{
static void NAG_CALL lsqgrd(const Integer &m,
const Integer &n,
const nagad_a1w_w_rtype fvec[],
const nagad_a1w_w_rtype fjac[],
const Integer &ldfjac,
nagad_a1w_w_rtype g[]);
static void NAG_CALL lsqfun(void* &ad_handle,
Integer &iflag,
const Integer &m,
const Integer &n,
const nagad_a1w_w_rtype xc[],
nagad_a1w_w_rtype fvec[],
Integer iuser[],
nagad_a1w_w_rtype ruser[]);
static void NAG_CALL lsqmon(void* &ad_handle,
const Integer &m,
const Integer &n,
const nagad_a1w_w_rtype xc[],
const nagad_a1w_w_rtype fvec[],
const nagad_a1w_w_rtype fjac[],
const Integer &ldfjac,
const nagad_a1w_w_rtype s[],
const Integer &igrade,
const Integer &niter,
const Integer &nf,
Integer iuser[],
nagad_a1w_w_rtype ruser[]);
}
int main(void)
{
// Scalars
int exit_status = 0;
const Integer m = 15, ldfjac = 15, n = 3, nt = 3, ldv = 3;
cout << "E04FC_A1W_F C++ Header Example Program Results\n\n";
// Create AD tape
nagad_a1w_ir_create();
// Create AD configuration data object
Integer ifail = 0;
void *ad_handle = 0;
x10aa_a1w_f_(ad_handle,ifail);
// Read problem parameters and register for differentiation
// Skip first line of data file
string mystr;
getline (cin, mystr);
double yr;
nagad_a1w_w_rtype ruser[m*nt+m];
for (int i=0; i<m; i++) {
cin >> yr;
ruser[i] = yr;
nagad_a1w_ir_register_variable(&ruser[i]);
for (int j=1; j<=nt; j++) {
cin >> yr;
Integer k = j*m + i;
ruser[k] = yr;
}
}
// AD routine arguments
Integer iprint, maxcal, niter, nf;
Integer iuser[1];
nagad_a1w_w_rtype eta, stepmx, xtol, fsumsq;
nagad_a1w_w_rtype x[n], fvec[m], fjac[m*n], g[n], s[n], v[ldv*n];
iprint = -1;
maxcal = 400*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 AD routine
ifail = -1;
e04fc_a1w_f_(ad_handle,m,n,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 << nagad_a1w_get_value(fsumsq);
cout << "\n Solution point = ";
for (int i=0; i<n; i++) {
cout.width(12); cout << nagad_a1w_get_value(x[i]);
}
lsqgrd(m,n,fvec,fjac,m,g);
cout.precision(3);
cout << "\n Estim gradient = ";
for (int i=0; i<n; i++) {
cout.width(12); cout << nagad_a1w_get_value(g[i]);
}
cout.precision(1);
cout << "\n Residuals :\n";
for (int i=0; i<m; i++) {
cout.width(12); cout << nagad_a1w_get_value(fvec[i]);
if (i%3==2) {
cout << endl;
}
}
}
cout << "\n Derivatives calculated: First order adjoints\n";
cout << " Computational mode : algorithmic\n\n";
cout << " Derivatives:\n\n";
// Setup evaluation of derivatives of fsumsq via adjoints.
double inc = 1.0;
nagad_a1w_inc_derivative(&fsumsq,inc);
ifail = 0;
nagad_a1w_ir_interpret_adjoint(ifail);
// Get derivatives of fsumsq
cout << " sum of squares w.r.t y:\n";
for (int i=0; i<m; i++) {
double d = nagad_a1w_get_derivative(ruser[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++) {
nagad_a1w_ir_zero_adjoints();
nagad_a1w_inc_derivative(&x[j],inc);
ifail = 0;
nagad_a1w_ir_interpret_adjoint(ifail);
cout << "\n dx(";
cout.width(1); cout << j+1;
cout << ")/dy :\n";
for (int i=0; i<m; i++) {
double d = nagad_a1w_get_derivative(ruser[i]);
cout.width(12); cout << d;
if (i%3==2) {
cout << endl;
}
}
}
// Remove computational data object and tape
x10ab_a1w_f_(ad_handle,ifail);
nagad_a1w_ir_remove();
return exit_status;
}
// dco/c++ used to perform AD of callbacks
static void NAG_CALL lsqgrd(const Integer &m,
const Integer &n,
const nagad_a1w_w_rtype fvec[],
const nagad_a1w_w_rtype fjac[],
const Integer &ldfjac,
nagad_a1w_w_rtype g[]){
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 lsqfun(void* &ad_handle,
Integer &iflag,
const Integer &m,
const Integer &n,
const nagad_a1w_w_rtype xc[],
nagad_a1w_w_rtype fvec[],
Integer iuser[],
nagad_a1w_w_rtype ruser[])
{
nagad_a1w_w_rtype x1, x2;
for (int i=0;i<m; i++) {
x1 = xc[1]*ruser[2*m+i];
x2 = xc[2]*ruser[3*m+i];
fvec[i] = ruser[m+i]/(x1+x2) + xc[0] - ruser[i];
}
return;
}
static void NAG_CALL lsqmon(void* &ad_handle,
const Integer &m,
const Integer &n,
const nagad_a1w_w_rtype xc[],
const nagad_a1w_w_rtype fvec[],
const nagad_a1w_w_rtype fjac[],
const Integer &ldfjac,
const nagad_a1w_w_rtype s[],
const Integer &igrade,
const Integer &niter,
const Integer &nf,
Integer iuser[],
nagad_a1w_w_rtype ruser[])
{
double fsumsq = 0.0;
for (int i=0;i<m; i++) {
double f = nagad_a1w_get_value(fvec[i]);
fsumsq = fsumsq + f*f;
}
nagad_a1w_w_rtype g[3];
lsqgrd(m,n,fvec,fjac,ldfjac,g);
double gtg = 0.0;
for (int i=0;i<n; i++) {
double gg = nagad_a1w_get_value(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 << nagad_a1w_get_value(xc[i]);
cout.width(19); cout << nagad_a1w_get_value(g[i]);
cout.width(19); cout << nagad_a1w_get_value(s[i]);
cout << endl;
}
return;
}