/* nag_fit_dim1_spline_deriv_vector (e02bfc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */

#include <nag.h>
#include <stdio.h>

int main(void)
{
#define S(I,J) s[(J-1)*pds + I-1]
  Integer exit_status = 0;
  double fp, sfac;
  Integer pds, liwrk, m, nest, nx, d, j;
  double *s = 0, *wdata = 0, *x = 0, *xdata = 0, *ydata = 0;
  Integer *iwrk = 0, *ixloc = 0;
  Nag_Comm warmstartinf;
  Nag_Spline spline;
  Nag_Start start_e02bec;
  Nag_SplineVectorSort start;
  Nag_Boolean xord;
  Nag_DerivType deriv;
  NagError fail;

  printf("nag_fit_dim1_spline_deriv_vector (e02bfc) Example Program Results\n");

  INIT_FAIL(fail);

  /* Initialize spline */
  spline.lamda = 0;
  spline.c = 0;
  warmstartinf.nag_w = 0;
  warmstartinf.nag_iw = 0;

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Input the number of data points for the spline, */
  /* followed by the data points (xdata), the function values (ydata) */
  /* and the weights (wdata). */
  scanf("%" NAG_IFMT "", &m);
  scanf("%*[^\n] ");
  nest = m + 4;
  if (m >= 4) {
    if (!(wdata = NAG_ALLOC(m, double)) ||
        !(xdata = NAG_ALLOC(m, double)) || !(ydata = NAG_ALLOC(m, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    printf("Invalid m.\n");
    exit_status = 1;
    return exit_status;
  }
  start_e02bec = Nag_Cold;

  for (j = 0; j < m; j++) {
    scanf("%lf", &xdata[j]);
    scanf("%lf", &ydata[j]);
    scanf("%lf", &wdata[j]);
  }
  scanf("%*[^\n] ");

  /* Read in the requested smoothing factor. */
  scanf("%lf", &sfac);
  scanf("%*[^\n] ");

  /* Determine the spline approximation.
   * nag_fit_dim1_spline_auto (e02bec).
   * Least squares cubic spline curve fit, automatic knot placement,
   * one variable.
   */
  nag_fit_dim1_spline_auto(start_e02bec, m, xdata, ydata, wdata, sfac, nest,
                    &fp, &warmstartinf, &spline, &fail);

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_dim1_spline_auto (e02bec).\n%s\n", fail.message);
    exit_status = 2;
    goto END;
  }

  /* Read in the number of sample points requested. */
  scanf("%" NAG_IFMT "", &nx);
  scanf("%*[^\n] ");

  /* Allocate memory for sample point locations and */
  /* function and derivative approximations. */
  pds = nx;
  liwrk = 3 + 3 * nx;
  if (!(x = NAG_ALLOC(nx, double)) ||
      !(s = NAG_ALLOC(pds * 4, double)) ||
      !(ixloc = NAG_ALLOC(nx, Integer)) || !(iwrk = NAG_ALLOC(liwrk, Integer))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in sample points. */
  for (j = 0; j < nx; j++)
    scanf("%lf", &x[j]);
  scanf("%*[^\n] ");

  xord = Nag_FALSE;
  start = Nag_SplineVectorSort_Sorted;
  deriv = Nag_RightDerivs_3;
  /*
   * nag_fit_dim1_spline_deriv_vector (e02bfc).
   * Evaluation of fitted cubic spline, function and optionally derivatives
   * at a vector of points.
   */
  nag_fit_dim1_spline_deriv_vector(start, &spline, deriv, xord, x, ixloc, nx,
                                s, pds, iwrk, liwrk, &fail);
  switch (fail.code) {
  case NE_NOERROR:
  case NW_SOME_SOLUTIONS:
    {
      /* Output the results. */
      printf("\n");
      printf("       x   ixloc         s(x) ");
      printf("       ds/dx      d2s/dx2      d3s/dx3\n");
      for (j = 0; j < nx; j++) {
        if (ixloc[j] >= 4 && ixloc[j] <= spline.n - 3) {
          printf("%8.4f %7" NAG_IFMT " ", x[j], ixloc[j]);
          for (d = 0; d < 4; d++)
            printf("%12.4e ", S(j + 1, d + 1));
          printf("\n");
        }
        else
          printf("%f %" NAG_IFMT "\n", x[j], ixloc[j]);
      }
      break;
    }
  default:
    {
      printf("Error from nag_fit_dim1_spline_deriv_vector (e02bfc).\n%s\n",
             fail.message);
      exit_status = 3;
      goto END;
    }
  }
END:

  NAG_FREE(xdata);
  NAG_FREE(ydata);
  NAG_FREE(wdata);
  NAG_FREE(warmstartinf.nag_w);
  NAG_FREE(warmstartinf.nag_iw);
  NAG_FREE(spline.lamda);
  NAG_FREE(spline.c);
  NAG_FREE(x);
  NAG_FREE(ixloc);
  NAG_FREE(s);
  NAG_FREE(iwrk);

  return exit_status;
}