/* nag_1d_cheb_interp (e01aec) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 */

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

int main(void)
{
  /* Scalars */
  double   xmax, xmin;
  Integer  exit_status, i, pmax, ires, iy, j, k, m, n,
           itmin, itmax, num_iter;
  NagError fail;

  /* Arrays */
  double   *a = 0, *perf = 0, *x = 0, *y = 0;
  Integer  *p = 0;

  exit_status = 0;

  INIT_FAIL(fail);

  printf("nag_1d_cheb_interp (e01aec) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  while (scanf("%ld%lf%lf%*[^\n] ", &m, &xmin, &xmax) != EOF)
    {
      if (m > 0)
        {
          /* Allocate memory for p and x. */
          if (!(p = NAG_ALLOC(m, Integer)) ||
              !(x = NAG_ALLOC(m, double)))
            {
              printf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

          /* Read p, x and y arrays */
          n = 0;
          pmax = 0;
          for (i = 1; i <= m; ++i)
            {
              scanf("%ld%lf", &p[i-1], &x[i-1]);
              k = n + p[i-1] + 1;
              /* We need to extend array y as we go along */
              if (!(y = NAG_REALLOC(y, k, double)))
                {
                  printf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }
              for (j = n + 1; j <= k; ++j)
                scanf("%lf", &y[j-1]);
              scanf("%*[^\n] ");
              if (p[i-1] > pmax)
                pmax = p[i-1];
              n = k;
            }

          /* Allocate array a */
          if (!(a = NAG_ALLOC(n, double)) ||
              !(perf = NAG_ALLOC(pmax+n+1, double)))
            {
              printf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

          itmin = -1;
          itmax = -1;
          /* nag_1d_cheb_interp (e01aec).
           * Interpolating functions, polynomial interpolant, data may
           * include derivative values, one variable
           */
          nag_1d_cheb_interp(m, xmin, xmax, x, y, p, itmin, itmax, a, perf,
                             &num_iter, &fail);

          printf("\n");
          if (fail.code == NE_NOERROR)
            {
              printf("Total number of interpolating conditions ="
                      " %ld\n", n);
              printf("\n");
              printf("Interpolating polynomial\n");
              printf("\n");
              printf("   i    Chebyshev Coefficient a(i+1)\n");

              for (i = 1; i <= n; ++i)
                printf("%4ld%20.4f\n", i - 1, a[i-1]);

              printf("\n");

              printf("  x    R   Rth derivative    Residual\n");
              iy = 0;
              ires = pmax + 1;
              for (i = 1; i <= m; ++i)
                {
                  for (j = 1; j <= p[i-1] + 1; ++j)
                    {
                      ++iy;
                      ++ires;
                      if (j - 1 != 0)
                        printf("    %4ld%12.1f%17.6f\n",
                                j - 1, y[iy-1], perf[ires-1]);
                      else
                        printf("%4.1f   0%12.1f%17.6f\n",
                                x[i-1], y[iy-1], perf[ires-1]);
                    }
                }
            }
          else
            {
              printf("Error from nag_1d_cheb_interp (e01aec).\n%s\n",
                      fail.message);
              exit_status = 1;
            }
        }
    }

 END:
  NAG_FREE(a);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(p);
  NAG_FREE(perf);

  return exit_status;
}