/* nag_1d_spline_fit_knots (e02bac) Example Program. * * Copyright 1996 Numerical Algorithms Group. * * Mark 4, 1996. * * Mark 6 revised, 2000. * Mark 8 revised, 2004. */ #include #include #include #include int main(void) { Integer exit_status = 0, j, m, ncap, ncap7, r, wght; NagError fail; Nag_Spline spline; double fit, ss, *weights = 0, *x = 0, xarg, *y = 0; INIT_FAIL(fail); /* Initialise spline */ spline.lamda = 0; spline.c = 0; printf("nag_1d_spline_fit_knots (e02bac) Example Program Results\n"); scanf("%*[^\n]"); /* Skip heading in data file */ while (scanf("%ld", &m) != EOF) { if (m >= 4) { if (!(weights = NAG_ALLOC(m, double)) || !(x = NAG_ALLOC(m, double)) || !(y = NAG_ALLOC(m, double)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid m.\n"); exit_status = 1; goto END; } scanf("%ld%ld", &ncap, &wght); if (ncap > 0) { ncap7 = ncap+7; spline.n = ncap7; if (!(spline.lamda = NAG_ALLOC(ncap7, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid ncap.\n"); exit_status = 1; goto END; } for (j = 4; j < ncap+3; ++j) scanf("%lf", &(spline.lamda[j])); for (r = 0; r < m; ++r) { if (wght == 1) { scanf("%lf%lf", &x[r], &y[r]); weights[r] = 1.0; } else scanf("%lf%lf%lf", &x[r], &y[r], &weights[r]); } /* nag_1d_spline_fit_knots (e02bac). * Least-squares curve cubic spline fit (including * interpolation), one variable */ nag_1d_spline_fit_knots(m, x, y, weights, &ss, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_1d_spline_fit_knots (e02bac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\nNumber of distinct knots = %ld\n\n", ncap+1); printf("Distinct knots located at \n\n"); for (j = 3; j < ncap+4; j++) printf("%8.4f%s", spline.lamda[j], (j-3)%6 == 5 || j == ncap+3?"\n":" "); printf("\n\n J B-spline coeff c\n\n"); for (j = 0; j < ncap+3; ++j) printf(" %ld %13.4f\n", j+1, spline.c[j]); printf("\nResidual sum of squares = "); printf("%11.2e\n\n", ss); printf("Cubic spline approximation and residuals\n"); printf(" r Abscissa Weight Ordinate" " Spline Residual\n\n"); for (r = 0; r < m; ++r) { /* nag_1d_spline_evaluate (e02bbc). * Evaluation of fitted cubic spline, function only */ nag_1d_spline_evaluate(x[r], &fit, &spline, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_1d_spline_evaluate (e02bbc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("%3ld %11.4f %11.4f %11.4f %11.4f" " %10.1e\n", r+1, x[r], weights[r], y[r], fit, fit-y[r]); if (r < m-1) { xarg = (x[r] + x[r+1]) * 0.5; /* nag_1d_spline_evaluate (e02bbc), see above. */ nag_1d_spline_evaluate(xarg, &fit, &spline, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_1d_spline_evaluate (e02bbc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" %14.4f %33.4f\n", xarg, fit); } } /* Free memory used by spline */ NAG_FREE(spline.lamda); NAG_FREE(spline.c); END: NAG_FREE(weights); NAG_FREE(x); NAG_FREE(y); } return exit_status; }