/* E02BC_T1W_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>
#include <string>
using namespace std;

// typedef DCO_TYPE nagad_t1w_w_rtype
typedef double                   DCO_BASE_TYPE;
typedef dco::gt1s<DCO_BASE_TYPE> DCO_MODE;
typedef DCO_MODE::type           DCO_TYPE;

int main(void)
{
  Integer exit_status = 0;

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

  // Create AD handle
  Integer ifail     = 0;
  void *  ad_handle = 0;
  nag::ad::x10aa(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++)
        {

          nag::ad::e02bc(ad_handle, ncap7, lamda, c, x, left, s, ifail);

          if (ifail != 0)
            {
              printf("\nError from e02bc_t1w_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;
            }

          Integer indx = (left - 1) * m * ncap7 + i * ncap7;
          for (int l = 0; l < ncap7; l++)
            {
              // lamda[l].tangent = 1.0;
              dco::derivative(lamda[l]) = 1.0;

              nag::ad::e02bc(ad_handle, ncap7, lamda, c, x, left, s, ifail);

              // tmp = s[0].tangent;
              tmp            = dco::derivative(s[0]);
              dsdl[indx + l] = tmp;
              // lamda[l].tangent = 0.0;
              dco::derivative(lamda[l]) = 0.0;
            }
          for (int l = 0; l < ncap7; l++)
            {
              // c[l].tangent = 1.0;
              dco::derivative(c[l]) = 1.0;

              nag::ad::e02bc(ad_handle, ncap7, lamda, c, x, left, s, ifail);

              // tmp = s[0].tangent;
              tmp            = dco::derivative(s[0]);
              dsdc[indx + l] = tmp;
              // c[l].tangent = 0.0;
              dco::derivative(c[l]) = 0.0;
            }
          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);

  nag::ad::x10ab(ad_handle, ifail);

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

  return exit_status;
}