/* nag_1d_cheb_fit_constr (e02agc) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 */

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

int main(void)
{
  /* Scalars */
  double fiti, xmax, xmin;
  Integer exit_status, i, iy, j, k, h, m, mf, n, pda, stride;
  NagError fail;
  Nag_OrderType order;

  /* Arrays */
  double *a = 0, *s = 0, *w = 0, *resid = 0,
         *x = 0, *xf = 0, *y = 0, *yf = 0;
  Integer *p = 0;

#ifdef NAG_COLUMN_MAJOR
#define A(I,J) a[(J-1)*pda + I - 1]
  order = Nag_ColMajor;
#else
#define A(I,J) a[(I-1)*pda + J - 1]
  order = Nag_RowMajor;
#endif

  INIT_FAIL(fail);
  exit_status = 0;
  Vprintf("e02agc Example Program Results\n");

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

  while (scanf("%ld%*[^\n] ", &mf) != EOF)
    {
      if (mf > 0)
        {
          /* Allocate memory for p and xf. */
          if (!(p = NAG_ALLOC(mf, Integer)) ||
              !(xf = NAG_ALLOC(mf, double)))
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

          /* Read p, xf and yf arrays */
          iy = 1;
          n = mf;
          for (i = 0; i < mf; ++i)
            {
              Vscanf("%ld%lf", &p[i], &xf[i]);
              h = iy + p[i] + 1;
              /* We need to extend array yf as we go along */
              if (!(yf = NAG_REALLOC(yf, h - 1, double)))
                {
                  Vprintf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }
              for (j = iy-1; j < h - 1; ++j)
                Vscanf("%lf", &yf[j]);
              Vscanf("%*[^\n] ");
              n += p[i];
              iy = h;
            }
          Vscanf("%ld%*[^\n] ", &m);

          if (m > 0)
            {
              /* Allocate memory for x, y and w. */
              if (!(x = NAG_ALLOC(m, double)) ||
                  !(y = NAG_ALLOC(m, double)) ||
                  !(w = NAG_ALLOC(m, double)))
                {
                  Vprintf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }
              for (i = 0; i < m; ++i)
                Vscanf("%lf%lf%lf", &x[i], &y[i], &w[i]);
              Vscanf("%*[^\n] ");

              Vscanf("%ld%lf%lf%*[^\n] ", &k, &xmin, &xmax);
              pda = k + 1;

              /* Allocate arrays a, s and resid */
              if ( !(a = NAG_ALLOC((k + 1) * (k + 1), double)) ||
                   !(s = NAG_ALLOC((k + 1), double)) ||
                   !(resid = NAG_ALLOC(m, double))  )
                {
                  Vprintf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }

              e02agc(order, m, k, xmin, xmax, x, y, w, mf, xf, yf,
                     p, a, s, &n, resid, &fail);
              if (fail.code != NE_NOERROR)
                {
                  Vprintf("Error from e02agc.\n%s\n", fail.message);
                  exit_status = 1;
                  goto END;
                }

              Vprintf("\n");
              Vprintf("Degree  RMS residual\n");
              for (i = n; i <= k; ++i)
                Vprintf("%4ld%15.2e\n", i, s[i]);
              Vprintf("\n");

              Vprintf("Details of the fit of degree %2ld\n", k);
              Vprintf("\n");
              Vprintf("  Index   Coefficient\n");
              for (i = 0; i < k + 1; ++i)
                Vprintf("%6ld%11.4f\n", i, A(k+1, i+1));
              Vprintf("\n");

              Vprintf("     i      x(i)       y(i)       Fit     Residual\n");
              for (i = 0; i < m; ++i)
                {
                 /* Note that the coefficients of polynomial are stored in the 
                  * rows of A hence when the storage is in Nag_ColMajor order 
                  * then stride is the first dimension of A, k + 1.
                  * When the storage is in Nag_RowMajor order then stride is 1.
                  */
#ifdef NAG_COLUMN_MAJOR
                  stride = k + 1;
#else
                  stride = 1;
#endif
                  e02akc(k, xmin, xmax, &A(k+1, 1), stride, x[i], &fiti, &fail);
                  if (fail.code != NE_NOERROR)
                    {
                      Vprintf("Error from e02akc.\n%s\n", fail.message);
                      exit_status = 1;
                      goto END;
                    }
                  Vprintf("%6ld%11.4f%11.4f%11.4f", i, x[i], y[i], fiti);
                  Vprintf("%11.2e\n", fiti - y[i]);
                }
            }
        }
    }

 END:
  if (a) NAG_FREE(a);
  if (s) NAG_FREE(s);
  if (w) NAG_FREE(w);
  if (resid) NAG_FREE(resid);
  if (x) NAG_FREE(x);
  if (xf) NAG_FREE(xf);
  if (y) NAG_FREE(y);
  if (yf) NAG_FREE(yf);
  if (p) NAG_FREE(p);

  return exit_status;
}