/* nag_2d_spline_fit_panel (e02dac) Example Program. * * Copyright 2005 Numerical Algorithms Group. * * Mark 8, 2004. */ #include #include #include #include #include int main(int argc, char *argv[]) { FILE *fpin, *fpout; /* Initialized data */ char label[] = "XY"; /* Scalars */ double d__1, eps, sigma, sum, temp; Integer exit_status, 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; /* Check for command-line IO options */ fpin = nag_example_file_io(argc, argv, "-data", NULL); fpout = nag_example_file_io(argc, argv, "-results", NULL); fprintf(fpout, "nag_2d_spline_fit_panel (e02dac) Example Program Results\n"); /* Skip heading in data file */ fscanf(fpin, "%*[^\n] "); while (fscanf(fpin, "%lf", &eps) != EOF) { /* Read data, interchanging X and Y axes if PX.LT.PY */ fscanf(fpin, "%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)) ) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } } else { fprintf(fpout, "Invalid m.\n"); exit_status = 1; return exit_status; } fscanf(fpin, "%ld%ld%*[^\n] ", &px, &py); if (px >= 8 && py >= 8) { 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))) { fprintf(fpout, "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)) ) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } for (i = 1; i <= m; ++i) { fscanf(fpin, "%lf%lf%lf%lf", &y[i - 1], &x[i - 1], &f[i - 1], &w[i - 1]); } fscanf(fpin, "%*[^\n] "); if (py > 8) { for (j = 5; j <= py - 4; ++j) { fscanf(fpin, "%lf", &mu[j - 1]); } fscanf(fpin, "%*[^\n] "); } if (px > 8) { for (j = 5; j <= px - 4; ++j) { fscanf(fpin, "%lf", &lamda[j - 1]); } fscanf(fpin, "%*[^\n] "); } } else { /* Allocate memory */ if (!(lamda = NAG_ALLOC(px, double)) || !(mu = NAG_ALLOC(py, double)) ) { fprintf(fpout, "Allocation failure\n"); exit_status = -1; goto END; } itemp = 0; for (i = 1; i <= m; ++i) { fscanf(fpin, "%lf%lf%lf%lf", &x[i - 1], &y[i - 1], &f[i - 1], &w[i - 1]); } fscanf(fpin, "%*[^\n] "); if (px > 8) { for (j = 5; j <= px - 4; ++j) { fscanf(fpin, "%lf", &lamda[j - 1]); } fscanf(fpin, "%*[^\n] "); } if (py > 8) { for (j = 5; j <= py - 4; ++j) { fscanf(fpin, "%lf", &mu[j - 1]); } fscanf(fpin, "%*[^\n] "); } } fprintf(fpout, "\n%s%1.1s %s\n", "Interior ", label + itemp, "-knots"); for (j = 5; j <= px - 4; ++j) { fprintf(fpout, "%11.4f\n", lamda[j - 1]); } if (px == 8) { fprintf(fpout, "%s\n", "None"); } fprintf(fpout, "\n%s%1.1s %s\n", "Interior ", label + (2 - itemp - 1), "-knots"); for (j = 5; j <= py - 4; ++j) { fprintf(fpout, "%1s%11.4f\n", "", mu[j - 1]); } if (py == 8) { fprintf(fpout, "%s\n", "None"); } /* Sort points into panel order */ /* 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) { fprintf(fpout, "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))) { fprintf(fpout, "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) { fprintf(fpout, "Error from nag_2d_spline_fit_panel (e02dac).\n%s\n", fail.message); exit_status = 1; goto END; } fprintf(fpout, "\n%s%16.3e\n", "Sum of squares of residual RHS", sigma); fprintf(fpout, "\n%s%5ld\n", "Rank", rank); /* Evaluate spline at the data points */ /* 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) { fprintf(fpout, "Error from nag_2d_spline_eval (e02dec).\n%s\n", fail.message); exit_status = 1; goto END; } sum = 0.; if (itemp == 1) { fprintf(fpout, "\n%s\n", "X and Y have been interchanged"); } /*Output data points, fitted values and residuals */ fprintf(fpout, "\n%s", " X Y Data Fit" " Residual"); fprintf(fpout, "\n"); for (i = 1; i <= np; ++i) { iadres = i + m; while ((iadres = point[iadres - 1]) > 0) { temp = ff[iadres - 1] - f[iadres - 1]; fprintf(fpout, "%11.4f%11.4f%11.4f%11.4f%12.3e\n", x[iadres - 1], y[iadres - 1], f[iadres - 1], ff[iadres - 1], temp); /* Computing 2nd power */ d__1 = temp * w[iadres - 1]; sum += d__1 * d__1; } } fprintf(fpout, "\n%s%16.3e\n", "Sum of squared residuals", sum); fprintf(fpout, "\n%s\n", "Spline coefficients"); for (i = 1; i <= px - 4; ++i) { for (j = 1; j <= py - 4; ++j) { fprintf(fpout, "%11.4f", spline.c[(i - 1) * (py - 4) + j - 1]); } fprintf(fpout, "\n"); } } END: if (fpin != stdin) fclose(fpin); if (fpout != stdout) fclose(fpout); if (dl) NAG_FREE(dl); if (f) NAG_FREE(f); if (ff) NAG_FREE(ff); if (lamda) NAG_FREE(lamda); if (mu) NAG_FREE(mu); if (w) NAG_FREE(w); if (x) NAG_FREE(x); if (y) NAG_FREE(y); if (point) NAG_FREE(point); if (spline.lamda) NAG_FREE(spline.lamda); if (spline.mu) NAG_FREE(spline.mu); if (spline.c) NAG_FREE(spline.c); } return exit_status; }