/* nag_1d_cheb_interp (e01aec) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include 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.3f\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: 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; }