/* nag_1d_cheb_fit_constr (e02agc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include 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; }