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