/* nag_2d_spline_deriv_rect (e02dhc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 23, 2011.
*/
#include <math.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage02.h>
#include <nagx04.h>
#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);
}