1 #define bard_N 15
2 #define bard_P 3
3
4 static double bard_x0[bard_P] = { 1.0, 1.0, 1.0 };
5 static double bard_epsrel = 1.0e-7;
6
7 static double bard_Y[bard_N] = {
8 0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37,
9 0.58, 0.73, 0.96, 1.34, 2.10, 4.39
10 };
11
12 static void
bard_checksol(const double x[],const double sumsq,const double epsrel,const char * sname,const char * pname)13 bard_checksol(const double x[], const double sumsq,
14 const double epsrel, const char *sname,
15 const char *pname)
16 {
17 size_t i;
18 const double sumsq_exact1 = 8.214877306578963e-03;
19 double bard_x1[bard_P] = { 8.241055975623580e-02,
20 1.133036092245175,
21 2.343695178435405 };
22 const double sumsq_exact2 = 17.42869333333333;
23 double bard_x2[bard_P] = { 8.406666666666666e-01,
24 0.0, /* -inf */
25 0.0 }; /* -inf */
26 double *bard_x;
27 double sumsq_exact;
28
29 bard_x2[1] = GSL_NAN;
30 bard_x2[2] = GSL_NAN;
31
32 if (fabs(x[1]) < 10.0 && fabs(x[2]) < 10.0)
33 {
34 bard_x = bard_x1;
35 sumsq_exact = sumsq_exact1;
36 }
37 else
38 {
39 bard_x = bard_x2;
40 sumsq_exact = sumsq_exact2;
41 }
42
43 gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
44 sname, pname);
45
46 for (i = 0; i < bard_P; ++i)
47 {
48 if (!gsl_finite(bard_x[i]))
49 continue;
50
51 gsl_test_rel(x[i], bard_x[i], epsrel, "%s/%s i=%zu",
52 sname, pname, i);
53 }
54 }
55
56 static int
bard_f(const gsl_vector * x,void * params,gsl_vector * f)57 bard_f (const gsl_vector * x, void *params, gsl_vector * f)
58 {
59 double x1 = gsl_vector_get(x, 0);
60 double x2 = gsl_vector_get(x, 1);
61 double x3 = gsl_vector_get(x, 2);
62 size_t i;
63
64 for (i = 0; i < bard_N; ++i)
65 {
66 double ui = i + 1.0;
67 double vi = 16.0 - i - 1.0;
68 double wi = GSL_MIN(ui, vi);
69 double yi = bard_Y[i];
70 double fi = yi - (x1 + (ui / (x2*vi + x3*wi)));
71
72 gsl_vector_set(f, i, fi);
73 }
74
75 (void)params; /* avoid unused parameter warning */
76
77 return GSL_SUCCESS;
78 }
79
80 static int
bard_df(const gsl_vector * x,void * params,gsl_matrix * J)81 bard_df (const gsl_vector * x, void *params, gsl_matrix * J)
82 {
83 double x2 = gsl_vector_get(x, 1);
84 double x3 = gsl_vector_get(x, 2);
85 size_t i;
86
87 for (i = 0; i < bard_N; ++i)
88 {
89 double ui = i + 1.0;
90 double vi = 16.0 - i - 1.0;
91 double wi = GSL_MIN(ui, vi);
92 double term = x2 * vi + x3 * wi;
93
94 gsl_matrix_set(J, i, 0, -1.0);
95 gsl_matrix_set(J, i, 1, ui * vi / (term * term));
96 gsl_matrix_set(J, i, 2, ui * wi / (term * term));
97 }
98
99 (void)params; /* avoid unused parameter warning */
100
101 return GSL_SUCCESS;
102 }
103
104 static int
bard_fvv(const gsl_vector * x,const gsl_vector * v,void * params,gsl_vector * fvv)105 bard_fvv (const gsl_vector * x, const gsl_vector * v,
106 void *params, gsl_vector * fvv)
107 {
108 double x2 = gsl_vector_get(x, 1);
109 double x3 = gsl_vector_get(x, 2);
110 double v2 = gsl_vector_get(v, 1);
111 double v3 = gsl_vector_get(v, 2);
112 size_t i;
113
114 for (i = 0; i < bard_N; ++i)
115 {
116 double ui = i + 1.0;
117 double vi = 16.0 - i - 1.0;
118 double wi = GSL_MIN(ui, vi);
119 double term1 = x2 * vi + x3 * wi;
120 double term2 = v2 * vi + v3 * wi;
121 double ratio = term2 / term1;
122
123 gsl_vector_set(fvv, i, -2.0 * ui * ratio * ratio / term1);
124 }
125
126 (void)params; /* avoid unused parameter warning */
127
128 return GSL_SUCCESS;
129 }
130
131 static gsl_multifit_nlinear_fdf bard_func =
132 {
133 bard_f,
134 bard_df,
135 bard_fvv,
136 bard_N,
137 bard_P,
138 NULL,
139 0,
140 0,
141 0
142 };
143
144 static test_fdf_problem bard_problem =
145 {
146 "bard",
147 bard_x0,
148 NULL,
149 NULL,
150 &bard_epsrel,
151 &bard_checksol,
152 &bard_func
153 };
154