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

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

// typedef DCO_TYPE nagad_a1w_w_rtype
typedef double DCO_BASE_TYPE;
typedef dco::ga1s<DCO_BASE_TYPE> DCO_MODE;
typedef DCO_MODE::type DCO_TYPE;
typedef DCO_MODE::tape_t DCO_TAPE_TYPE;

int main(void)
{
  Integer exit_status = 0;

  cout << "E02BC_A1W_F C++ Header Example Program\n";

  // Create AD handle
  Integer ifail = 0;
  void    *ad_handle = 0;
  x10aa_a1w_f_(ad_handle,ifail);

  // Skip heading in data file
  string mystr;
  getline (cin, mystr);
  
  Integer ncap, m;
  cin >> ncap; cin >> m;

  Integer ncap7 = ncap + 7;

  if (m <= 0) {
    printf("Invalid m.\n");
    exit_status = 1;
    return exit_status;
  }
  if (ncap <= 0) {
    printf("Invalid ncap.\n");
    exit_status = 1;
    return exit_status;
  }
  
  // Allocate arrays
  DCO_TYPE  *c = 0, *lamda = 0;
  double    *dsdl = 0, *dsdc = 0;
  c     = new DCO_TYPE [ncap7];
  lamda = new DCO_TYPE [ncap7];
  dsdc  = new double [2*m*ncap7];
  dsdl  = new double [2*m*ncap7];

  // Read knots and spline coefficients from file
  double tmp;
  for (int j = 0; j < ncap7; j++) {
    cin >> tmp;
    lamda[j] = tmp;
  }
  for (int j = 0; j < ncap + 3; j++) {
    cin >> tmp;
    c[j] = tmp;
  }

  cout << "\n        x             Spline    1st deriv";
  cout << "   2nd deriv   3rd deriv\n";
  cout.setf(ios::scientific,ios::floatfield);
  cout.setf(ios::right);
  cout.precision(4);

  DCO_TYPE s[4], x;
  for (int i = 0; i < m; i++) {
    cin >> tmp;
    x = tmp;

    for (Integer left = 1; left <= 2; left++) {
      // Create AD tape
      // nagad_a1w_ir_create();
      DCO_MODE::global_tape=DCO_TAPE_TYPE::create();
      for (int l = 0; l < ncap7; l++) {
        // nagad_a1w_ir_register_variable(&lamda[l]);
        DCO_MODE::global_tape->register_variable(lamda[l]);
      }
      for (int l = 0; l < ncap7; l++) {
        // nagad_a1w_ir_register_variable(&c[l]);
        DCO_MODE::global_tape->register_variable(c[l]);
      }

      e02bc_a1w_f_(ad_handle, ncap7, lamda, c, x, left, s, ifail);
      
      if (ifail != 0) {
        printf("\nError from e02bc_a1w_f_.\n%" NAG_IFMT "\n", ifail);
        exit_status = 1;
        return exit_status;
      }

      // double xv = x.value;
      double xv = dco::value(x);
      cout.width(11); cout << xv;

      if (left == 1) {
        cout << "   Left";
      } else {
        cout << "  Right";
      }
      for (int l = 0; l < 4; l++) {
        // tmp = s[l].value;
        tmp = dco::value(s[l]);
        cout.width(12); cout << tmp;
      }
      // s[0].tangent = 1.0;
      dco::derivative(s[0]) = 1.0;
      // nagad_a1w_ir_interpret_adjoint(ifail);
      DCO_MODE::global_tape->interpret_adjoint();

      Integer indx = (left-1)*m*ncap7 + i*ncap7;
      for (int l = 0; l < ncap7; l++) {
        // tmp = lamda[l].tangent;
        tmp = dco::derivative(lamda[l]);
        dsdl[indx + l] = tmp;
        // tmp = c[l].tangent;
        tmp = dco::derivative(c[l]);
        dsdc[indx + l] = tmp;
      }
      // nagad_a1w_ir_remove();
      DCO_TAPE_TYPE::remove(DCO_MODE::global_tape);
      cout << "\n";
    }
    cout << "\n";
  }
  
  cout << "\n Derivatives of spline values with respect to knots\n";
  NagError  fail;
  INIT_FAIL(fail);
  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,dsdl,ncap7,
         "Left:  ds/dlamda",0,&fail);
  cout << "\n";
  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,&dsdl[m*ncap7],
         ncap7, "Right: ds/dlamda",0,&fail);
  cout << "\n\n Derivatives with respect to spline coefficients\n";

  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,dsdc,ncap7,
         "Left:  ds/dC",0,&fail);
  cout << "\n";
  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,&dsdc[m*ncap7],
         ncap7, "Right: ds/dC",0,&fail);

  x10ab_a1w_f_(ad_handle,ifail);

  delete [] c;
  delete [] lamda;
  delete [] dsdc;
  delete [] dsdl;

  return exit_status;
}