/* 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;
  printf("nag_1d_cheb_fit_constr (e02agc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\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)))
            {
              printf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

          /* Read p, xf and yf arrays */
          iy = 1;
          n = mf;
          for (i = 0; i < mf; ++i)
            {
              scanf("%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)))
                {
                  printf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }
              for (j = iy-1; j < h - 1; ++j)
                scanf("%lf", &yf[j]);
              scanf("%*[^\n] ");
              n += p[i];
              iy = h;
            }
          scanf("%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)))
                {
                  printf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }
              for (i = 0; i < m; ++i)
                scanf("%lf%lf%lf", &x[i], &y[i], &w[i]);
              scanf("%*[^\n] ");

              scanf("%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)))
                {
                  printf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }

              /* nag_1d_cheb_fit_constr (e02agc).
               * Least-squares polynomial fit, values and derivatives may
               * be constrained, arbitrary data points
               */
              nag_1d_cheb_fit_constr(order, m, k, xmin, xmax, x, y, w, mf, xf,
                                     yf, p, a, s, &n, resid, &fail);
              if (fail.code != NE_NOERROR)
                {
                  printf(
                          "Error from nag_1d_cheb_fit_constr (e02agc).\n%s\n",
                          fail.message);
                  exit_status = 1;
                  goto END;
                }

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

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

              printf(
                      "     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
                  /* nag_1d_cheb_eval2 (e02akc).
                   * Evaluation of fitted polynomial in one variable from
                   * Chebyshev series form
                   */
                  nag_1d_cheb_eval2(k, xmin, xmax, &A(k+1, 1), stride, x[i],
                                    &fiti, &fail);
                  if (fail.code != NE_NOERROR)
                    {
                      printf(
                              "Error from nag_1d_cheb_eval2 (e02akc).\n%s\n",
                              fail.message);
                      exit_status = 1;
                      goto END;
                    }
                  printf("%6ld%11.4f%11.4f%11.4f", i, x[i], y[i],
                          fiti);
                  printf("%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;
}