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