/* nag_2d_spline_deriv_rect (e02dhc) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL print_spline(Integer *ngx, double *gridx, Integer *ngy, double *gridy, double *z, double *zder, Integer *exit_status); #ifdef __cplusplus } #endif #define F(I, J) f[my*(I)+(J)] int main(void) { /* Scalars */ Integer exit_status = 0; Integer i, j, mx, my, ngx, ngy, nux, nuy, nxest, nyest; double delta, fp, s, xhi, xlo, yhi, ylo; /* Arrays */ double *f = 0, *gridx = 0, *gridy = 0, *x = 0, *y = 0, *z = 0, *zder = 0; /* NAG types */ Nag_2dSpline spline; Nag_Comm warmstartinf; Nag_Start startc; NagError fail; INIT_FAIL(fail); printf("nag_2d_spline_deriv_rect (e02dhc) Example Program Results\n"); fflush(stdout); /* Skip heading in data file*/ scanf("%*[^\n] "); /* Input the number of X, Y co-ordinates MX, MY.*/ scanf("%ld %ld%*[^\n]", &mx, &my); nxest = mx + 4; nyest = my + 4; spline.nx = 4; spline.ny = 4; /* Alocations for spline fitting */ if (!(x = NAG_ALLOC((mx), double)) || !(y = NAG_ALLOC((my), double)) || !(f = NAG_ALLOC((mx * my), double)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 0; i < mx; i++) scanf("%lf", &x[i]); scanf("%*[^\n]"); for (i = 0; i < my; i++) scanf("%lf", &y[i]); scanf("%*[^\n]"); /* Input the MX*MY function values F at grid points and smoothing factor.*/ for (i = 0; i < mx; i++) for (j = 0; j < my; j++) scanf("%lf", &F(i, j)); scanf("%*[^\n]"); scanf("%lf%*[^\n]", &s); /* nag_2d_spline_fit_grid (e02dcc). * Least squares bicubic spline fit with automatic knot placement, * two variables (rectangular grid) */ startc = Nag_Cold; nag_2d_spline_fit_grid(startc, mx, x, my, y, f, s, nxest, nyest, &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; } printf("\nSpline fit used smoothing factor s = %13.4e.\n", s); printf("Number of knots in each direction = %5ld,%5ld.\n\n", spline.nx, spline.ny); printf("Sum of squared residuals = %13.4e.\n", fp); fflush(stdout); /* Spline and its derivative to be evaluated on rectangular grid with * ngx*ngy points on the domain [xlo,xhi] by [ylo,yhi]. */ scanf("%ld%lf%lf%*[^\n]", &ngx, &xlo, &xhi); scanf("%ld%lf%lf%*[^\n]", &ngy, &ylo, &yhi); /* Allocations for spline evaluation.*/ if (!(gridx = NAG_ALLOC((ngx), double)) || !(gridy = NAG_ALLOC((ngy), double)) || !(z = NAG_ALLOC((ngx * ngy), double)) || !(zder = NAG_ALLOC((ngx * ngy), double)) ) { printf("Allocation failure\n"); exit_status = -2; goto END; } delta = (xhi - xlo)/(double) (ngx - 1); gridx[0] = xlo; for (i = 1; i < ngx - 1; i++) gridx[i] = gridx[i-1] + delta; gridx[ngx-1] = xhi; delta = (yhi - ylo)/(double) (ngy - 1); gridy[0] = ylo; for (i = 1; i < ngy - 1; i++) gridy[i] = gridy[i-1] + delta; gridy[ngy-1] = yhi; /* Evaluate spline (nux=nuy=0) using * nag_2d_spline_deriv_rect (e02dhc). * Evaluation of spline surface at mesh of points with derivatives */ nux = 0; nuy = 0; nag_2d_spline_deriv_rect(ngx, ngy, gridx, gridy, nux, nuy, z, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_2d_spline_deriv_rect (e02dhc))\n%s\n", fail.message); exit_status = 2; goto END; } /* Evaluate spline partial derivative of order (nux,nuy)*/ scanf("%ld%ld%*[^\n]", &nux, &nuy); printf("\nDerivative of spline has order nux, nuy =%5ld, %5ld.\n", nux, nuy); fflush(stdout); /* nag_2d_spline_deriv_rect (e02dhc). * Evaluation of spline surface at mesh of points with derivatives */ nag_2d_spline_deriv_rect(ngx, ngy, gridx, gridy, nux, nuy, zder, &spline, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_2d_spline_deriv_rect (e02dhc)\n%s\n", fail.message); exit_status = 3; goto END; } fflush(stdout); /* Print tabulated spline and derivative evaluations.*/ print_spline(&ngx, gridx, &ngy, gridy, z, zder, &exit_status); END: NAG_FREE(f); NAG_FREE(gridx); NAG_FREE(gridy); NAG_FREE(x); NAG_FREE(y); NAG_FREE(z); NAG_FREE(zder); NAG_FREE(spline.lamda); NAG_FREE(spline.mu); NAG_FREE(spline.c); NAG_FREE(warmstartinf.nag_w); NAG_FREE(warmstartinf.nag_iw); return exit_status; } static void NAG_CALL print_spline(Integer *ngx, double *gridx, Integer *ngy, double *gridy, double *z, double *zder, Integer *exit_status) { /* Print spline function and spline derivative evaluation*/ Integer indent = 0, ncols = 80; char formc[] = "%8.3f"; Integer i; char title[49]; char *outfile = 0; char **clabsc = 0, **rlabsc = 0; Nag_OrderType order; Nag_MatrixType matrixc; Nag_DiagType diagc; Nag_LabelType chlabelc; NagError fail; INIT_FAIL(fail); /* Allocate for row and column label*/ if ( !(clabsc = NAG_ALLOC(*ngx, char *)) || !(rlabsc = NAG_ALLOC(*ngy, char *)) ) { printf("Allocation failure\n"); *exit_status = -3; goto END; } /* Allocate memory to clabsc and rlabsc elements and generate * column and row labels to print the results with. */ for (i = 0; i < *ngx; i++) { clabsc[i] = NAG_ALLOC(11, char); sprintf(clabsc[i], "%5.2f%5s", gridx[i], ""); } for (i = 0; i < *ngy; i++) { rlabsc[i] = NAG_ALLOC(11, char); sprintf(rlabsc[i], "%5.2f%5s", gridy[i], ""); } order = Nag_ColMajor; matrixc = Nag_GeneralMatrix; diagc = Nag_NonUnitDiag; chlabelc = Nag_CharacterLabels; /* Print the spline evaluations, z. */ strcpy(title, "Spline evaluated on X-Y grid (X across, Y down):"); printf("\n"); fflush(stdout); /* nag_gen_real_mat_print_comp (x04cbc). * Print real general matrix (comprehensive) */ nag_gen_real_mat_print_comp(order, matrixc, diagc, *ngy, *ngx, z, *ngy, formc, title, chlabelc, (const char **) rlabsc, chlabelc, (const char **) clabsc, ncols, indent, outfile, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print_comp (x04cbc)\n%s\n", fail.message); *exit_status = 4; goto END; } /* Print the spline derivative evaluations, zder. */ strcpy(title, "Spline derivative evaluated on X-Y grid:"); printf("\n"); fflush(stdout); nag_gen_real_mat_print_comp(order, matrixc, diagc, *ngy, *ngx, zder, *ngy, formc, title, chlabelc, (const char **) rlabsc, chlabelc, (const char **) clabsc, ncols, indent, outfile, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print_comp (x04cbc)\n%s\n", fail.message); *exit_status = 5; goto END; } END: for (i = 0; i < *ngy; i++) NAG_FREE(rlabsc[i]); NAG_FREE(rlabsc); for (i = 0; i < *ngx; i++) NAG_FREE(clabsc[i]); NAG_FREE(clabsc); }