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

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

int main(void)
{
  int exit_status = 0;

  cout << "G01GC_A1W_F C++ Header Example Program Results\n\n";
  // Skip heading in data file
  string mystr;
  getline (cin, mystr);

  // Read number of x values
  Integer n;
  cin >> n;

  cout << "    x       df      rlamda     p";
  cout << "          p'         df'     rlamda'\n";
  cout.setf(ios::scientific,ios::floatfield);
  cout.setf(ios::right);
  cout.precision(2);

  nagad_a1w_w_rtype tol;
  tol.value = 0.0;
  tol.id = 0;
  Integer maxit = 100;

  nagad_a1w_w_rtype *x      = new nagad_a1w_w_rtype[n];
  nagad_a1w_w_rtype *df     = new nagad_a1w_w_rtype[n];
  nagad_a1w_w_rtype *rlamda = new nagad_a1w_w_rtype[n];
      
  // Create AD tape
  nagad_a1w_ir_create();

  // Create AD configuration data object
  Integer ifail = -1;
  void    *ad_handle = 0;
  x10aa_a1w_f_(ad_handle,ifail);
  if (ifail != 0) {
    exit_status = 1;
    exit(exit_status);
  }

  // Loop over x values
  for (Integer i = 0; i < n; ++i) {
    // Read next x, df, rlamda
    double xr, dfr, rr;
    cin >> xr;     cin >> dfr;      cin >> rr;

    x[i].value = xr;  df[i].value = dfr;  rlamda[i].value = rr;
    x[i].id = 0;      df[i].id = 0;       rlamda[i].id = 0;


    // Register variables, i.e. differentiate w.r.t. x, df and rlamda
    nagad_a1w_ir_register_variable(&x[i]);
    nagad_a1w_ir_register_variable(&df[i]);
    nagad_a1w_ir_register_variable(&rlamda[i]);

    // Call NAG AD Routine
    //    nagad_a1w_ir_zero_adjoints();
    nagad_a1w_w_rtype p;
    ifail = -1;
    g01gc_a1w_f_(ad_handle,x[i],df[i],rlamda[i],tol,maxit,p,ifail);
    if (ifail != 0) {
      exit_status = 2;
      goto END;
    }
    
    // Reset adjoints, increment p, and evaluate adjoints
    nagad_a1w_ir_zero_adjoints();
    double inc = 1.0;
    nagad_a1w_inc_derivative(&p,inc);
    ifail = 0;
    nagad_a1w_ir_interpret_adjoint_sparse(ifail);

    // Get derivatives
    double dpdx, dpdd, dpdr;
    dpdx = nagad_a1w_get_derivative(x[i]);
    dpdd = nagad_a1w_get_derivative(df[i]);
    dpdr = nagad_a1w_get_derivative(rlamda[i]);

    // Output results
    cout.width(9);  cout << x[i].value;
    cout.width(9);  cout << df[i].value;
    cout.width(9);  cout << rlamda[i].value;
    cout.width(9);  cout << p.value;
    cout.width(11);  cout << dpdx;
    cout.width(11);  cout << dpdd;
    cout.width(11);  cout << dpdr << endl;

    // Remove computational data object and tape
  }
  x10ab_a1w_f_(ad_handle,ifail);
  nagad_a1w_ir_remove();

 END:
  return exit_status;
}