/* nag_fit_1dspline_deriv_vector (e02bfc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage02.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_1dspline_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("%ld", &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_1d_spline_fit (e02bec).
   * Least squares cubic spline curve fit, automatic knot placement,
   * one variable.
   */
  nag_1d_spline_fit(start_e02bec, m, xdata, ydata, wdata, sfac, nest,
                    &fp, &warmstartinf, &spline, &fail);

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

  /* Read in the number of sample points requested.*/
  scanf("%ld", &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_1dspline_deriv_vector (e02bfc).
   * Evaluation of fitted cubic spline, function and optionally derivatives
   * at a vector of points.
   */
  nag_fit_1dspline_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 %7ld ",x[j],ixloc[j]);
                for (d=0; d<4; d++)
                  printf("%12.4e ", S(j+1,d+1));
                printf("\n");
              }
            else
              printf("%f %ld\n", x[j],ixloc[j]);
          }
        break;
      }
    default:
      {
        printf("Error from nag_fit_1dspline_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;
}