/* nag_2d_spline_fit_grid (e02dcc) 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, mx, my, npx, npy, nx, ny; Nag_2dSpline spline; Nag_Comm warmstartinf; Nag_Start start; double delta, *f = 0, *fg = 0, fp, *px = 0, *py = 0, s, *x = 0, xhi, xlo, *y = 0, yhi; double ylo; NagError fail; INIT_FAIL(fail); /* Initialise spline */ spline.lamda = 0; spline.mu = 0; spline.c = 0; warmstartinf.nag_w = 0; warmstartinf.nag_iw = 0; printf("nag_2d_spline_fit_grid (e02dcc) Example Program Results\n"); scanf("%*[^\n]"); /* Skip heading in data file */ /* Input the number of x, y co-ordinates mx, my. */ scanf("%ld%ld", &mx, &my); if (mx >= 4 && my >= 4) { if (!(f = NAG_ALLOC(mx*my, double)) || !(x = NAG_ALLOC(mx, double)) || !(y = NAG_ALLOC(my, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid mx or my.\n"); exit_status = 1; return exit_status; } /* Input the x co-ordinates followed by the y co-ordinates. */ for (i = 0; i < mx; i++) scanf("%lf", &x[i]); for (i = 0; i < my; i++) scanf("%lf", &y[i]); /* Input the mx*my function values f at the grid points. */ for (i = 0; i < mx*my; i++) scanf("%lf", &f[i]); start = Nag_Cold; scanf("%lf", &s); /* Determine the spline approximation. */ /* nag_2d_spline_fit_grid (e02dcc). * Least-squares bicubic spline fit with automatic knot * placement, two variables (rectangular grid) */ nag_2d_spline_fit_grid(start, mx, x, my, y, f, s, mx+4, my+4, &fp, &warmstartinf, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_2d_spline_fit_grid (e02dcc).\n%s\n", fail.message); exit_status = 1; goto END; } nx = spline.nx; ny = spline.ny; printf("\nCalling with smoothing factor s = %13.4e:" " spline.nx = %2ld, spline.ny = %2ld.\n\n", s, nx, ny); /* 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.4f", spline.c[i+j*(ny-4)]); printf("\n"); } printf("\nSum of squared residuals fp = %13.4e\n", fp); if (fp == 0.0) printf("\nThe spline is an interpolating spline\n"); else if (nx == 8 && ny == 8) printf("\nThe 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(xlo+i*delta, xhi); for (i = 0; i < npy; i++) py[i] = MIN(ylo+i*delta, yhi); /* 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"); printf("\n x"); for (i = 0; i < npx; i++) printf("%7.2f ", px[i]); printf("\n y\n"); for (i = npy-1; i >= 0; i--) { printf("%7.2f ", py[i]); for (j = 0; j < npx; ++j) printf("%8.2f ", fg[npy*j+i]); printf("\n"); } END: NAG_FREE(spline.lamda); NAG_FREE(spline.mu); NAG_FREE(spline.c); NAG_FREE(warmstartinf.nag_w); NAG_FREE(warmstartinf.nag_iw); NAG_FREE(f); NAG_FREE(x); NAG_FREE(y); NAG_FREE(fg); NAG_FREE(px); NAG_FREE(py); return exit_status; }