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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
/*
Divonne.c
Multidimensional integration by partitioning
originally by J.H. Friedman and M.H. Wright
(CERNLIB subroutine D151)
this version by Thomas Hahn
last modified 15 Nov 11 th
*/
#include "decl.h"
#define Print(s) puts(s); fflush(stdout)
/*********************************************************************/
#define DIVONNE
static void DoSample(This *t, number n, creal *x, real *f, ccount ldx);
static int ExploreParent(This *t, cint iregion);
/*********************************************************************/
static inline count SampleExtra(This *t, cBounds *b)
{
number n = t->nextra;
t->peakfinder(&t->ndim, b, &n, t->xextra);
DoSample(t, n, t->xextra, t->fextra, t->ldxgiven);
return n;
}
/*********************************************************************/
static inline void AllocGiven(This *t, creal *xgiven)
{
if( t->ngiven | t->nextra ) {
cnumber nxgiven = t->ngiven*(t->ldxgiven = IMax(t->ldxgiven, t->ndim));
cnumber nxextra = t->nextra*t->ldxgiven;
cnumber nfgiven = t->ngiven*t->ncomp;
cnumber nfextra = t->nextra*t->ncomp;
Alloc(t->xgiven, nxgiven + nxextra + nfgiven + nfextra);
t->xextra = t->xgiven + nxgiven;
t->fgiven = t->xextra + nxextra;
t->fextra = t->fgiven + nfgiven;
if( nxgiven ) {
t->phase = 0;
Copy(t->xgiven, xgiven, nxgiven);
DoSample(t, t->ngiven, t->xgiven, t->fgiven, t->ldxgiven);
}
}
}
/*********************************************************************/
#include "common.c"
#include "DoSample.c"
Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
Integrand integrand, void *userdata,
creal epsrel, creal epsabs,
cint flags, cint seed,
cnumber mineval, cnumber maxeval,
cint key1, cint key2, cint key3, ccount maxpass,
creal border, creal maxchisq, creal mindeviation,
cnumber ngiven, ccount ldxgiven, creal *xgiven,
cnumber nextra, PeakFinder peakfinder,
int *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
{
This t;
t.ndim = ndim;
t.ncomp = ncomp;
t.integrand = integrand;
t.userdata = userdata;
t.epsrel = epsrel;
t.epsabs = epsabs;
t.flags = flags;
t.seed = seed;
t.mineval = mineval;
t.maxeval = maxeval;
t.key1 = key1;
t.key2 = key2;
t.key3 = key3;
t.maxpass = maxpass;
t.border.upper = 1 - (t.border.lower = border);
t.maxchisq = maxchisq;
t.mindeviation = mindeviation;
t.ngiven = ngiven;
t.xgiven = NULL;
t.ldxgiven = ldxgiven;
t.nextra = nextra;
t.peakfinder = peakfinder;
t.nregions = 0;
t.neval = 0;
ForkCores(&t);
AllocGiven(&t, xgiven);
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;
*pneval = t.neval;
free(t.xgiven);
WaitCores(&t);
}
/*********************************************************************/
Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
Integrand integrand, void *userdata,
creal *pepsrel, creal *pepsabs,
cint *pflags, cint *pseed,
cnumber *pmineval, cnumber *pmaxeval,
cint *pkey1, cint *pkey2, cint *pkey3, ccount *pmaxpass,
creal *pborder, creal *pmaxchisq, creal *pmindeviation,
cnumber *pngiven, ccount *pldxgiven, creal *xgiven,
cnumber *pnextra, PeakFinder peakfinder,
int *pnregions, number *pneval, int *pfail,
real *integral, real *error, real *prob)
{
This t;
t.ndim = *pndim;
t.ncomp = *pncomp;
t.integrand = integrand;
t.userdata = userdata;
t.epsrel = *pepsrel;
t.epsabs = *pepsabs;
t.flags = *pflags;
t.seed = *pseed;
t.mineval = *pmineval;
t.maxeval = *pmaxeval;
t.key1 = *pkey1;
t.key2 = *pkey2;
t.key3 = *pkey3;
t.maxpass = *pmaxpass;
t.border.upper = 1 - (t.border.lower = *pborder);
t.maxchisq = *pmaxchisq;
t.mindeviation = *pmindeviation;
t.ngiven = *pngiven;
t.xgiven = NULL;
t.ldxgiven = *pldxgiven;
t.nextra = *pnextra;
t.peakfinder = peakfinder;
t.nregions = 0;
t.neval = 0;
ForkCores(&t);
AllocGiven(&t, xgiven);
*pfail = Integrate(&t, integral, error, prob);
*pnregions = t.nregions;
*pneval = t.neval;
free(t.xgiven);
WaitCores(&t);
}
|
