/* 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; 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; }