/* nag_2d_spline_fit_panel (e02dac) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 8, 2004.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage02.h>
int main(void)
{
/* Initialized data */
char label[] = "xy";
/* Scalars */
double d, eps, sigma, sum, temp;
Integer exit_status=0, i, iadres, itemp, j, m, nc, np, npoint, px, py,
rank;
/* Arrays */
double *dl = 0, *f = 0, *ff = 0, *lamda = 0, *mu = 0, *w = 0, *x = 0;
double *y = 0;
Integer *point = 0;
/* Nag Types */
Nag_2dSpline spline;
NagError fail;
exit_status = 0;
INIT_FAIL(fail);
/* Initialise spline */
spline.lamda = 0;
spline.mu = 0;
spline.c = 0;
printf("nag_2d_spline_fit_panel (e02dac) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
while (scanf("%lf", &eps) != EOF && exit_status==0)
{
/* Read data, interchanging X and Y axes if PX.LT.PY */
scanf("%ld%*[^\n] ", &m);
if (m > 1)
{
/* Allocate memory */
if (!(f = NAG_ALLOC(m, double)) ||
!(ff = NAG_ALLOC(m, double)) ||
!(w = 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;
goto END;
}
scanf("%ld%ld%*[^\n] ", &px, &py);
if (px < 8 && py < 8) {
printf("px or py is too small.\n");
exit_status = 1;
goto END;
}
nc = (px - 4) * (py - 4);
np = (px - 7) * (py - 7);
npoint = m+(px-7)*(py-7);
/* Allocate memory */
if (!(dl = NAG_ALLOC(nc, double)) ||
!(point = NAG_ALLOC(npoint, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
if (px < py)
{
itemp = px;
px = py;
py = itemp;
itemp = 1;
/* Allocate memory */
if (!(lamda = NAG_ALLOC(px, double)) ||
!(mu = NAG_ALLOC(py, double))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < m; ++i)
scanf("%lf%lf%lf%lf", &y[i], &x[i], &f[i], &w[i]);
scanf("%*[^\n] ");
if (py > 8)
{
for (j = 4; j < py - 4; ++j) scanf("%lf", &mu[j]);
scanf("%*[^\n] ");
}
if (px > 8)
{
for (j = 4; j < px - 4; ++j) scanf("%lf", &lamda[j]);
scanf("%*[^\n] ");
}
}
else
{
/* Allocate memory */
if (!(lamda = NAG_ALLOC(px, double)) ||
!(mu = NAG_ALLOC(py, double))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
itemp = 0;
for (i = 0; i < m; ++i)
scanf("%lf%lf%lf%lf", &x[i], &y[i], &f[i], &w[i]);
scanf("%*[^\n] ");
if (px > 8)
{
for (j = 4; j < px - 4; ++j) scanf("%lf", &lamda[j]);
scanf("%*[^\n] ");
}
if (py > 8)
{
for (j = 4; j < py - 4; ++j) scanf("%lf", &mu[j]);
scanf("%*[^\n] ");
}
}
printf("\nInterior %1.1s -knots\n", label + itemp);
for (j = 4; j < px - 4; ++j) printf("%11.4f\n", lamda[j]);
if (px == 8) printf("None\n");
printf("\nInterior %1.1s -knots\n", label + (2 - itemp - 1));
for (j = 4; j < py - 4; ++j) printf("%1s%11.4f\n", "", mu[j]);
if (py == 8) printf("None\n");
/* nag_2d_panel_sort (e02zac).
* Sort two-dimensional data into panels for fitting bicubic
* splines
*/
nag_2d_panel_sort(px, py, lamda, mu, m, x, y, point, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_2d_panel_sort (e02zac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Fit bicubic spline to data points */
spline.nx = px;
spline.ny = py;
if (!(spline.c = NAG_ALLOC((spline.nx-4)*(spline.ny-4), double)) ||
!(spline.lamda = NAG_ALLOC(spline.nx, double)) ||
!(spline.mu = NAG_ALLOC(spline.ny, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < spline.nx; i++) spline.lamda[i] = lamda[i];
for (i = 0; i < spline.ny; i++) spline.mu[i] = mu[i];
/* nag_2d_spline_fit_panel (e02dac).
* Least-squares surface fit, bicubic splines
*/
nag_2d_spline_fit_panel(m, x, y, f, w, point, dl, eps, &sigma, &rank,
&spline, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_2d_spline_fit_panel (e02dac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\nSum of squares of residual RHS%16.3e\n", sigma);
printf("\nRank%5ld\n", rank);
/* nag_2d_spline_eval (e02dec).
* Evaluation of bicubic spline, at a set of points
*/
nag_2d_spline_eval(m, x, y, ff, &spline, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_2d_spline_eval (e02dec).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
sum = 0.;
if (itemp == 1) printf("\nx and y have been interchanged\n\n");
/*Output data points, fitted values and residuals */
printf(" X Y Data Fit Residual\n");
for (i = 0; i < np; ++i)
{
iadres = i + m;
while ((iadres = point[iadres]-1) >= 0)
{
temp = ff[iadres] - f[iadres];
printf("%11.4f%11.4f%11.4f%11.4f%12.3e\n", x[iadres],
y[iadres], f[iadres], ff[iadres], temp);
/* Computing 2nd power */
d = temp * w[iadres];
sum += d * d;
}
}
printf("\nSum of squared residuals%16.3e\n", sum);
printf("\nSpline coefficients\n");
for (i = 0; i < px - 4; ++i)
{
for (j = 0; j < py - 4; ++j)
printf("%11.4f", spline.c[i * (py - 4) + j]);
printf("\n");
}
END:
NAG_FREE(dl);
NAG_FREE(f);
NAG_FREE(ff);
NAG_FREE(lamda);
NAG_FREE(mu);
NAG_FREE(w);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(point);
NAG_FREE(spline.lamda);
NAG_FREE(spline.mu);
NAG_FREE(spline.c);
}
return exit_status;
}