/* nag_1d_spline_fit (e02bec) Example Program. * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * * Mark 6 revised, 2000. * Mark 8 revised, 2004. */ #include #include #include #include int main(void) { Integer exit_status = 0, j, m, nest, r; NagError fail; Nag_Comm warmstartinf; Nag_Spline spline; Nag_Start start; double fp, s, *sp = 0, txr, *weights = 0, *x = 0, *y = 0; INIT_FAIL(fail); /* Initialise spline */ spline.lamda = 0; spline.c = 0; warmstartinf.nag_w = 0; warmstartinf.nag_iw = 0; printf("nag_1d_spline_fit (e02bec) Example Program Results\n"); scanf("%*[^\n]"); /* Skip heading in data file */ /* Input the number of data points, followed by the data * points x, the function values y and the weights w. */ scanf("%ld", &m); nest = m + 4; if (m >= 4) { if (!(weights = NAG_ALLOC(m, double)) || !(x = NAG_ALLOC(m, double)) || !(y = NAG_ALLOC(m, double)) || !(sp = NAG_ALLOC(2*m-1, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid m.\n"); exit_status = 1; return exit_status; } start = Nag_Cold; for (r = 0; r < m; r++) scanf("%lf%lf%lf", &x[r], &y[r], &weights[r]); /* Read in successive values of s until end of data file. */ while (scanf("%lf", &s) != EOF) { /* Determine the spline approximation. */ /* nag_1d_spline_fit (e02bec). * Least-squares cubic spline curve fit, automatic knot * placement, one variable */ nag_1d_spline_fit(start, m, x, y, weights, s, nest, &fp, &warmstartinf, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_1d_spline_fit (e02bec).\n%s\n", fail.message); exit_status = 1; goto END; } /* Evaluate the spline at each x point and midway * between x points, saving the results in sp. */ for (r = 0; r < m; r++) { /* nag_1d_spline_evaluate (e02bbc). * Evaluation of fitted cubic spline, function only */ nag_1d_spline_evaluate(x[r], &sp[(r-1)*2+2], &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_1d_spline_fit (e02bec).\n%s\n", fail.message); exit_status = 1; goto END; } } for (r = 0; r < m-1; r++) { txr = (x[r] + x[r+1]) / 2; /* nag_1d_spline_evaluate (e02bbc), see above. */ nag_1d_spline_evaluate(txr, &sp[r*2+1], &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; } } /* Output the results. */ printf("\nCalling with smoothing factor s = %12.3e\n", s); printf("\nNumber of distinct knots = %ld\n\n", spline.n-6); printf("Distinct knots located at \n\n"); for (j = 3; j < spline.n-3; j++) printf("%8.4f%s", spline.lamda[j], (j-3)%6 == 5 || j == spline.n-4?"\n":" "); printf("\n\n J B-spline coeff c\n\n"); for (j = 0; j < spline.n-4; ++j) printf(" %3ld %13.4f\n", j+1, spline.c[j]); printf("\nWeighted sum of squared residuals fp = %12.3e\n", fp); if (fp == 0.0) printf("The spline is an interpolating spline\n"); else if (spline.n == 8) printf("The spline is the weighted least-squares cubic" "polynomial\n"); start = Nag_Warm; } /* Free memory allocated in spline and warmstartinf */ END: NAG_FREE(spline.lamda); NAG_FREE(spline.c); NAG_FREE(warmstartinf.nag_w); NAG_FREE(warmstartinf.nag_iw); NAG_FREE(weights); NAG_FREE(x); NAG_FREE(y); NAG_FREE(sp); return exit_status; }