/* D01UA_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 <nagad.h>
#include <stdio.h>
using namespace std;

extern "C"
{
  static void NAG_CALL f_a1w(void *&                 ad_handle,
                             const nagad_a1w_w_rtype x[],
                             const Integer &         nx,
                             nagad_a1w_w_rtype       fv[],
                             Integer &               iflag,
                             Integer                 iuser[],
                             nagad_a1w_w_rtype       ruser[]);
}

int main(void)
{
  // Scalars
  int     exit_status = 0;
  Integer n           = 64;
  // Arrays
  nagad_a1w_w_rtype ruser[1];
  Integer           iuser[1];

  ruser[0] = 1.0;

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

  // 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;
  nag::ad::x10aa(ad_handle, ifail);

  // Register variables to differentiate w.r.t.
  dco::ga1s<double>::global_tape->register_variable(ruser[0]);

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

  nagad_a1w_w_rtype a, b;
  Integer           key;
  for (int funid = 1; funid <= 6; funid++)
    {
      if (funid == 1)
        {
          a   = 0.0;
          b   = 1.0;
          key = 0;
        }
      if (funid == 2)
        {
          a   = 2.0;
          b   = 0.0;
          key = -5;
        }
      if (funid == 3)
        {
          a   = 2.0;
          b   = 1.0;
          key = -3;
        }
      if (funid == 4)
        {
          a   = 2.0;
          b   = 1.0;
          key = 3;
        }
      if (funid == 5)
        {
          a   = -1.0;
          b   = 3.0;
          key = -4;
        }
      if (funid == 6)
        {
          a   = -1.0;
          b   = 3.0;
          key = 4;
        }

      iuser[0] = funid;

      // Call the AD routine
      ifail = 0;
      nagad_a1w_w_rtype dinest;
      nag::ad::d01ua(ad_handle, key, a, b, n, f_a1w, dinest, -1, iuser, -1,
                     ruser, ifail);

      dco::ga1s<double>::global_tape->zero_adjoints();
      double inc = 1.0;
      dco::derivative(dinest) += inc;
      dco::ga1s<double>::global_tape->sparse_interpret() = true;
      dco::ga1s<double>::global_tape->interpret_adjoint();

      // Get derivatives
      double dr;
      dr = dco::derivative(ruser[0]);

      cout.setf(ios::scientific, ios::floatfield);
      cout.setf(ios::right);
      cout.precision(4);
      cout << "\n dinest = ";
      double dinest_value = dco::value(dinest);
      cout.width(12);
      cout << dinest_value << "     ";
      cout << " d/druser[0] = ";
      cout.width(12);
      cout << dr << endl;
    }

  // Remove computational data object and tape
  nag::ad::x10ab(ad_handle, ifail);
  dco::ga1s<double>::tape_t::remove(dco::ga1s<double>::global_tape);

  return exit_status;
}

static void NAG_CALL f_a1w(void *&                 ad_handle,
                           const nagad_a1w_w_rtype x[],
                           const Integer &         nx,
                           nagad_a1w_w_rtype       fv[],
                           Integer &               iflag,
                           Integer                 iuser[],
                           nagad_a1w_w_rtype       ruser[])
{
  // dco/c++ overloading used here to perform AD
  Integer           i;
  nagad_a1w_w_rtype xx, fx, rux;
  for (i = 0; i < nx && iflag >= 0; i++)
    {
      xx  = x[i];
      rux = ruser[0] * xx;
      if (iuser[0] == 1)
        {
          fx = 4.0 / (1.0 + ruser[0] * xx * xx);
        }
      if (iuser[0] == 2)
        {
          fx = 1.0 / (xx * xx * log(ruser[0] * xx));
        }
      if (iuser[0] == 3)
        {
          fx = exp(-ruser[0] * xx) / xx;
        }
      if (iuser[0] == 4)
        {
          fx = ruser[0] / xx;
        }
      if (iuser[0] == 5)
        {
          fx = exp(-3.0 * ruser[0] * xx * xx - 4.0 * xx - 1.0);
        }
      if (iuser[0] == 6)
        {
          fx = exp(2.0 * ruser[0] * xx + 2.0);
        }
      if (iuser[0] < 0 || iuser[0] > 6)
        {
          iflag = -1;
        }
      fv[i] = fx;
    }

  return;
}