1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
for (I = problems ; I->f != 0; I++)
{
size_t i;
double res, err;
gsl_rng * r;
if (I->dim > 3)
{
continue ;
}
r = gsl_rng_alloc (gsl_rng_default);
for (i = 0; i < TRIALS ; i++)
{
MONTE_STATE *s = MONTE_ALLOC (I->dim);
#ifdef MONTE_PARAMS
MONTE_PARAMS params;
#endif
I->f->dim = I->dim;
MONTE_INTEGRATE (I->f, I->xl, I->xu,
I->dim, I->calls / MONTE_SPEEDUP, r, s,
&res, &err);
gsl_test_abs (res, I->expected_result,
5 * GSL_MAX(err, 1024*GSL_DBL_EPSILON),
NAME ", %s, result[%d]", I->description, i);
MONTE_ERROR_TEST (err, I->expected_error);
result[i] = res;
error[i] = err;
MONTE_FREE (s);
}
/* Check the results for consistency as an ensemble */
{
double mean = 0, sumd2 = 0, sd;
/* We need to compute the mean exactly when all terms are equal,
to get an exact zero for the standard deviation (this is a
common case when integrating a constant). */
for (i = 0; i < TRIALS; i++)
{
mean += (result[i] - mean) / (i + 1.0);
}
for (i = 0; i < TRIALS; i++)
{
sumd2 += pow(result[i] - mean, 2.0);
}
sd = sqrt(sumd2 / (TRIALS-1.0)) ;
for (i = 0; i < TRIALS; i++)
{
gsl_test_factor (error[i], sd, 5.0,
NAME ", %s, abserr[%d] vs sd", I->description, i);
}
}
gsl_rng_free (r);
}
|
