Revision 1099
branches/threaded2/examples/threaded/mforest/qromb.c (revision 1099)  

1 
#include <math.h> 

2 
#define EPS 1.0e6 

3 
#define JMAX 20 

4 
#define JMAXP (JMAX+1) 

5 
#define K 5 

6  
7 
float qromb(float (*func)(float), float a, float b) 

8 
{ 

9 
void polint(float xa[], float ya[], int n, float x, float *y, float *dy); 

10 
float trapzd(float (*func)(float), float a, float b, int n); 

11 
void nrerror(char error_text[]); 

12 
float ss,dss; 

13 
float s[JMAXP],h[JMAXP+1]; 

14 
int j; 

15  
16 
h[1]=1.0; 

17 
for (j=1;j<=JMAX;j++) { 

18 
s[j]=trapzd(func,a,b,j); 

19 
if (j >= K) { 

20 
polint(&h[jK],&s[jK],K,0.0,&ss,&dss); 

21 
if (fabs(dss) <= EPS*fabs(ss)) return ss; 

22 
} 

23 
h[j+1]=0.25*h[j]; 

24 
} 

25 
nrerror("Too many steps in routine qromb"); 

26 
return 0.0; 

27 
} 

28 
#undef EPS 

29 
#undef JMAX 

30 
#undef JMAXP 

31 
#undef K 
branches/threaded2/examples/threaded/mforest/mforestnothreads.tpl (revision 1099)  

1 
GLOBALS_SECTION 

2 
#include "trace.h" 

3  
4 
DATA_SECTION 

5 
!! ad_comm::change_datafile_name("mforest.dat"); 

6 
init_int nsteps 

7 
init_int k 

8 
init_vector a(1,k+1) 

9 
init_vector freq(1,k) 

10 
int a_index; 

11 
number sum_freq 

12 
!! sum_freq=sum(freq); 

13 
PARAMETER_SECTION 

14 
init_bounded_number log_tau(14,15,2) 

15 
init_bounded_number log_nu(15,4) 

16 
init_bounded_number beta(.1,1.0,1) 

17 
init_bounded_number log_sigma(5,3) 

18 
likeprof_number tau 

19 
!! tau.set_stepnumber(25); 

20 
sdreport_number nu 

21 
sdreport_number sigma 

22 
vector S(1,k+1) 

23 
objective_function_value f 

24 
INITIALIZATION_SECTION 

25 
log_tau 0 

26 
beta 0.6666667 

27 
log_nu 0 

28 
log_sigma 2 

29 
PROCEDURE_SECTION 

30 
//cout << endl; 

31 
//cout << "ifn = " << ifn << endl; 

32 
//cout << "quit_flag = " << quit_flag << endl; 

33 
//cout << "ihflag = " << ihflag << endl; 

34 
//cout << "last_phase() " << last_phase() << endl; 

35 
//cout << "iexit = " << iexit << endl; 

36 
tau=exp(log_tau); 

37 
nu=exp(log_nu); 

38 
sigma=exp(log_sigma); 

39 
funnel_dvariable Integral; 

40 
int i; 

41 
for (i=1;i<=k+1;i++) 

42 
{ 

43 
a_index=i; 

44 
ad_begin_funnel(); 

45 
Integral=adromb(&model_parameters::h,3.0,3.0,nsteps); 

46 
//TRACE(Integral) 

47 
S(i)=Integral; 

48 
} 

49 
f=0.0; 

50 
for (i=1;i<=k;i++) 
Also available in: Unified diff