/* nag_2d_spline_fit_scat (e02ddc) 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, i, j, m, npx, npy, nx, nxest, ny, nyest, rank; NagError fail; Nag_2dSpline spline; Nag_Start start; double delta, *f = 0, *fg = 0, fp, *px = 0, *py = 0, s, warmstartinf; double *weights = 0, *x = 0, xhi, xlo, *y = 0, yhi, ylo; INIT_FAIL(fail); /* Initialise spline */ spline.lamda = 0; spline.mu = 0; spline.c = 0; nxest = 14; nyest = 14; printf("nag_2d_spline_fit_scat (e02ddc) Example Program Results\n"); scanf("%*[^\n]"); /* Skip heading in data file */ /* Input the number of data-points m. */ scanf("%ld", &m); if (m >= 16) { if (!(f = NAG_ALLOC(m, double)) || !(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; return exit_status; } /* Input the data-points and the weights. */ for (i = 0; i < m; i++) scanf("%lf%lf%lf%lf", &x[i], &y[i], &f[i], &weights[i]); start = Nag_Cold; if (scanf("%lf", &s) != EOF) { /* Determine the spline approximation. */ /* nag_2d_spline_fit_scat (e02ddc). * Least-squares bicubic spline fit with automatic knot * placement, two variables (scattered data) */ nag_2d_spline_fit_scat(start, m, x, y, f, weights, s, nxest, nyest, &fp, &rank, &warmstartinf, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_2d_spline_fit_scat (e02ddc).\n%s\n", fail.message); exit_status = 1; goto END; } nx = spline.nx; ny = spline.ny; printf("\nCalling with smoothing factor s = %13.4e, nx = %1ld," " ny = %1ld\n", s, nx, ny); printf("rank deficiency = %1ld\n\n", (nx-4)*(ny-4)-rank); /* Print the knot sets, lamda and mu. */ printf("Distinct knots in x direction located at\n"); for (j = 3; j < spline.nx-3; j++) printf("%12.4f%s", spline.lamda[j], ((j-3)%5 == 4 || j == spline.nx-4)?"\n":" "); printf("\nDistinct knots in y direction located at\n"); for (j = 3; j < spline.ny-3; j++) printf("%12.4f%s", spline.mu[j], ((j-3)%5 == 4 || j == spline.ny-4)?"\n":" "); printf("\nThe B-spline coefficients:\n\n"); for (i = 0; i < ny-4; i++) { for (j = 0; j < nx-4; j++) printf("%9.2f", spline.c[i+j*(ny-4)]); printf("\n"); } printf("\n Sum of squared residuals fp = %13.4e\n", fp); if (nx == 8 && ny == 8) printf("The spline is the least-squares bi-cubic polynomial\n"); /* Evaluate the spline on a rectangular grid at npx*npy points * over the domain (xlo to xhi) x (ylo to yhi). */ scanf("%ld%lf%lf", &npx, &xlo, &xhi); scanf("%ld%lf%lf", &npy, &ylo, &yhi); if (npx >= 1 && npy >= 1) { if (!(fg = NAG_ALLOC(npx*npy, double)) || !(px = NAG_ALLOC(npx, double)) || !(py = NAG_ALLOC(npy, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid npx or npy.\n"); exit_status = 1; return exit_status; } delta = (xhi-xlo)/(npx-1); for (i = 0; i < npx; i++) px[i] = MIN(xhi, xlo+i*delta); for (i = 0; i < npy; i++) py[i] = MIN(yhi, ylo+i*delta); /* nag_2d_spline_eval_rect (e02dfc). * Evaluation of bicubic spline, at a mesh of points */ nag_2d_spline_eval_rect(npx, npy, px, py, fg, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_2d_spline_eval_rect (e02dfc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\nValues of computed spline:\n\n"); printf(" x"); for (i = 0; i < npx; i++) printf("%8.2f ", px[i]); printf("\n y\n"); for (i = npy-1; i >= 0; i--) { printf("%8.2f ", py[i]); for (j = 0; j < npx; j++) printf("%8.2f ", fg[npy*j+i]); printf("\n"); } /* Free memory used by spline */ NAG_FREE(spline.lamda); NAG_FREE(spline.mu); NAG_FREE(spline.c); if (fg) NAG_FREE(fg); if (px) NAG_FREE(px); if (py) NAG_FREE(py); } END: if (f) NAG_FREE(f); if (weights) NAG_FREE(weights); if (x) NAG_FREE(x); if (y) NAG_FREE(y); return exit_status; }