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

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

int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;
  /* 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);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  fprintf(fpout, "nag_1d_cheb_interp (e01aec) Example Program Results\n");

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

  while (fscanf(fpin, "%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)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }

          /* Read p, x and y arrays */
          n = 0;
          pmax = 0;
          for (i = 1; i <= m; ++i)
            {
              fscanf(fpin, "%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)))
                {
                  fprintf(fpout, "Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }
              for (j = n + 1; j <= k; ++j)
                fscanf(fpin, "%lf", &y[j-1]);
              fscanf(fpin, "%*[^\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)))
            {
              fprintf(fpout, "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);

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

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

              fprintf(fpout, "\n");

              fprintf(fpout, "  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)
                        fprintf(fpout, "    %4ld%12.1f%17.6f\n",
                                j - 1, y[iy-1], perf[ires-1]);
                      else
                        fprintf(fpout, "%4.1f   0%12.1f%17.6f\n",
                                x[i-1], y[iy-1], perf[ires-1]);
                    }
                }
            }
          else
            {
              fprintf(fpout, "Error from nag_1d_cheb_interp (e01aec).\n%s\n",
                      fail.message);
              exit_status = 1;
            }
        }
    }

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (a) NAG_FREE(a);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (p) NAG_FREE(p);
  if (perf) NAG_FREE(perf);

  return exit_status;
}