NAG Library Manual, Mark 27.2
Interfaces:  FL   CL   CPP   AD
```/* E01AA_A1W_F C++ Header Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
* Mark 27.2, 2021.
*/

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

int main(void)
{
// Scalars
int exit_status = 0;

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

// Skip first line of data file
string mystr;
getline(cin, mystr);
// Read number of data points
Integer n;
cin >> n;

// Allocate arrays for data and interpolant
nagad_a1w_w_rtype x, *a = 0, *b = 0, *c = 0;
double            xr, ar;
Integer           n1, n2;
n1 = n + 1;
n2 = n * (n + 1) / 2;
if (!(a = NAG_ALLOC(n1, nagad_a1w_w_rtype)) ||
!(b = NAG_ALLOC(n1, nagad_a1w_w_rtype)) ||
!(c = NAG_ALLOC(n2, nagad_a1w_w_rtype)))
{
printf("Allocation failure\n");
exit_status = -1;
}
if (exit_status == 0)
{
// Create AD tape
dco::ga1s<double>::global_tape = dco::ga1s<double>::tape_t::create();

cin >> xr;
x = xr;

for (int i = 0; i < n1; i++)
{
cin >> ar;
a[i] = ar;
}
for (int i = 0; i < n1; i++)
{
cin >> ar;
b[i] = ar;
}

// Create AD tape
dco::ga1s<double>::global_tape = dco::ga1s<double>::tape_t::create();

// Create AD configuration data object
Integer ifail     = 0;
void *  ad_handle = 0;

// Register variables to differentiate w.r.t.
for (int i = 0; i < n1; i++)
{
dco::ga1s<double>::global_tape->register_variable(a[i]);
}

// Call the AD routine
ifail = 0;
nag::ad::e01aa(ad_handle, a, b, c, n, x, ifail);

cout.setf(ios::scientific, ios::floatfield);
cout.precision(5);
cout << "\n Interpolation point = ";
cout.width(12);
cout << xr << endl;
cout << "\n Function value at interpolation point = ";
cout.width(12);
cout << dco::value(c[n2 - 1]) << endl;

// Setup evaluation of derivatives via adjoints.
double inc = 1.0;
dco::derivative(c[n2 - 1]) += inc;

ifail                                              = 0;
dco::ga1s<double>::global_tape->sparse_interpret() = true;

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

// Get derivatives
cout << "\n Derivatives of fitted value w.r.t. x values:\n";
cout << "  j         d/da(j)\n";
for (int j = 0; j < n1; j++)
{
double dx = dco::derivative(a[j]);
cout.width(3);
cout << j + 1 << "     ";
cout.width(12);
cout << dx << endl;
}

// Remove computational data object and tape