#define thurber_N 37
#define thurber_P 7
static double thurber_x0a[thurber_P] = { 1000.0, 1000.0, 400.0, 40.0,
0.7, 0.3, 0.03 };
static double thurber_x0b[thurber_P] = { 1300.0, 1500.0, 500.0, 75.0,
1.0, 0.4, 0.05 };
static double thurber_epsrel = 1.0e-6;
static double thurber_sigma[thurber_P] = {
4.6647963344E+00, 3.9571156086E+01, 2.8698696102E+01,
5.5675370270E+00, 3.1333340687E-02, 1.4984928198E-02,
6.5842344623E-03
};
static double thurber_X[thurber_N] = {
-3.067, -2.981, -2.921, -2.912, -2.840, -2.797, -2.702,
-2.699, -2.633, -2.481, -2.363, -2.322, -1.501, -1.460,
-1.274, -1.212, -1.100, -1.046, -0.915, -0.714, -0.566,
-0.545, -0.400, -0.309, -0.109, -0.103, 0.010, 0.119,
0.377, 0.790, 0.963, 1.006, 1.115, 1.572, 1.841,
2.047, 2.200
};
static double thurber_F[thurber_N] = {
80.574, 84.248, 87.264, 87.195, 89.076, 89.608, 89.868, 90.101,
92.405, 95.854, 100.696, 101.060, 401.672, 390.724, 567.534,
635.316, 733.054, 759.087, 894.206, 990.785, 1090.109, 1080.914,
1122.643, 1178.351, 1260.531, 1273.514, 1288.339, 1327.543, 1353.863,
1414.509, 1425.208, 1421.384, 1442.962, 1464.350, 1468.705,
1447.894, 1457.628
};
static void
thurber_checksol(const double x[], const double sumsq,
const double epsrel, const char *sname,
const char *pname)
{
size_t i;
const double sumsq_exact = 5.6427082397E+03;
const double thurber_x[thurber_P] = {
1.2881396800E+03, 1.4910792535E+03, 5.8323836877E+02,
7.5416644291E+01, 9.6629502864E-01, 3.9797285797E-01,
4.9727297349E-02 };
gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
sname, pname);
for (i = 0; i < thurber_P; ++i)
{
gsl_test_rel(x[i], thurber_x[i], epsrel, "%s/%s i=%zu",
sname, pname, i);
}
}
static int
thurber_f (const gsl_vector * x, void *params, gsl_vector * f)
{
double b[thurber_P];
size_t i;
for (i = 0; i < thurber_P; i++)
{
b[i] = gsl_vector_get(x, i);
}
for (i = 0; i < thurber_N; i++)
{
double xi = thurber_X[i];
double yi;
yi = b[0] + b[1]*xi + b[2]*xi*xi + b[3]*xi*xi*xi;
yi /= 1.0 + b[4]*xi + b[5]*xi*xi + b[6]*xi*xi*xi;
gsl_vector_set (f, i, yi - thurber_F[i]);
}
(void)params; /* avoid unused parameter warning */
return GSL_SUCCESS;
}
static int
thurber_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
double b[thurber_P];
size_t i;
for (i = 0; i < thurber_P; i++)
{
b[i] = gsl_vector_get(x, i);
}
for (i = 0; i < thurber_N; i++)
{
double xi = thurber_X[i];
double d, n, d_sq;
n = b[0] + b[1]*xi + b[2]*xi*xi + b[3]*xi*xi*xi;
d = 1.0 + b[4]*xi + b[5]*xi*xi + b[6]*xi*xi*xi;
d_sq = d * d;
gsl_matrix_set (df, i, 0, 1.0 / d);
gsl_matrix_set (df, i, 1, xi / d);
gsl_matrix_set (df, i, 2, (xi * xi) / d);
gsl_matrix_set (df, i, 3, (xi * xi * xi) / d);
gsl_matrix_set (df, i, 4, -xi * n / d_sq);
gsl_matrix_set (df, i, 5, -xi * xi * n / d_sq);
gsl_matrix_set (df, i, 6, -xi * xi * xi * n / d_sq);
}
(void)params; /* avoid unused parameter warning */
return GSL_SUCCESS;
}
static gsl_multifit_nlinear_fdf thurber_func =
{
thurber_f,
thurber_df,
NULL, /* analytic expression too complex */
thurber_N,
thurber_P,
NULL,
0,
0,
0
};
static test_fdf_problem thurbera_problem =
{
"nist-thurbera",
thurber_x0a,
NULL,
thurber_sigma,
&thurber_epsrel,
&thurber_checksol,
&thurber_func
};
static test_fdf_problem thurberb_problem =
{
"nist-thurberb",
thurber_x0b,
NULL,
thurber_sigma,
&thurber_epsrel,
&thurber_checksol,
&thurber_func
};