Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
test_cube.c
Go to the documentation of this file.
1
2/*
3 Usage: ./test_cube <dim> <tol> <integrand> <maxeval>
4
5 where <dim> = # dimensions, <tol> = relative tolerance,
6 <integrand> is either 0/1/2 for the three test integrands (see below),
7 and <maxeval> is the maximum # function evaluations (0 for none).
8*/
9
10/* author Steven G. Johnson */
11
12#include "cubature.h"
13#include <math.h>
14#include <gsl/gsl_math.h> /* M_PI */
15
16int count = 0;
18const double radius = 0.50124145262344534123412; /* random */
19
20double f_test(unsigned dim, const double *x, void *data)
21{
22 double val;
23 unsigned i;
24 ++count;
25 switch (which_integrand) {
26 case 0: /* discontinuous objective: volume of hypersphere */
27 val = 0;
28 for (i = 0; i < dim; ++i)
29 val += x[i] * x[i];
30 val = val < radius * radius;
31 break;
32 case 1: /* simple smooth (separable) objective: prod. cos(x[i]). */
33 val = 1;
34 for (i = 0; i < dim; ++i)
35 val *= cos(x[i]);
36 break;
37 case 2: { /* integral of exp(-x^2), rescaled to (0,infinity) limits */
38 double scale = 1.0;
39 val = 0;
40 for (i = 0; i < dim; ++i) {
41 double z = (1 - x[i]) / x[i];
42 val += z * z;
43 scale *= M_2_SQRTPI / (x[i] * x[i]);
44 }
45 val = exp(-val) * scale;
46 break;
47 }
48
49 default:
50 fprintf(stderr, "unknown integrand %d\n", which_integrand);
51 exit(EXIT_FAILURE);
52 }
53 /* if (count < 100) printf("%d: f(%g, ...) = %g\n", count, x[0], val); */
54 return val;
55}
56
57/* surface area of n-dimensional unit hypersphere */
58static double S(unsigned n)
59{
60 double val;
61 int fact = 1;
62 if (n % 2 == 0) { /* n even */
63 val = 2 * pow(M_PI, n * 0.5);
64 n = n / 2;
65 while (n > 1) fact *= (n -= 1);
66 val /= fact;
67 }
68 else { /* n odd */
69 val = (1 << (n/2 + 1)) * pow(M_PI, n/2);
70 while (n > 2) fact *= (n -= 2);
71 val /= fact;
72 }
73 return val;
74}
75
76static double exact_integral(unsigned dim, const double *xmax) {
77 unsigned i;
78 double val;
79 switch(which_integrand) {
80 case 0:
81 val = dim == 0 ? 1 : S(dim) * pow(radius * 0.5, dim) / dim;
82 break;
83 case 1:
84 val = 1;
85 for (i = 0; i < dim; ++i)
86 val *= sin(xmax[i]);
87 break;
88 case 2:
89 val = 1;
90 break;
91 default:
92 fprintf(stderr, "unknown integrand %d\n", which_integrand);
93 exit(EXIT_FAILURE);
94 }
95 return val;
96}
97
98int main(int argc, char **argv)
99{
100 double *xmin, *xmax;
101 double tol, val, err;
102 unsigned i, dim, maxEval;
103
104 dim = argc > 1 ? atoi(argv[1]) : 2;
105 tol = argc > 2 ? atof(argv[2]) : 1e-2;
106 which_integrand = argc > 3 ? atoi(argv[3]) : 1;
107 maxEval = argc > 4 ? atoi(argv[4]) : 0;
108
109 xmin = (double *) malloc(dim * sizeof(double));
110 xmax = (double *) malloc(dim * sizeof(double));
111 for (i = 0; i < dim; ++i) {
112 xmin[i] = 0;
113 xmax[i] = 1 + (which_integrand == 2 ? 0 : 0.4 * sin(i));
114 }
115
116 printf("%u-dim integral, tolerance = %g, integrand = %d\n",
117 dim, tol, which_integrand);
118 adapt_integrate(f_test, 0, dim, xmin, xmax, maxEval, 0, tol, &val, &err);
119 printf("integration val = %g, est. err = %g, true err = %g\n",
120 val, err, fabs(val - exact_integral(dim, xmax)));
121 printf("#evals = %d\n", count);
122
123 free(xmax);
124 free(xmin);
125
126 return 0;
127}
128
int adapt_integrate(integrand f, void *fdata, unsigned dim, const double *xmin, const double *xmax, unsigned maxEval, double reqAbsError, double reqRelError, double *val, double *estimated_error)
Definition cubature.c:701
#define n
Definition exactTestM.c:37
static double S(unsigned n)
Definition test_cube.c:58
const double radius
Definition test_cube.c:18
int main(int argc, char **argv)
Definition test_cube.c:98
static double exact_integral(unsigned dim, const double *xmax)
Definition test_cube.c:76
int which_integrand
Definition test_cube.c:17
int count
Definition test_cube.c:16
double f_test(unsigned dim, const double *x, void *data)
Definition test_cube.c:20