Revision 1099

branches/threaded2/examples/threaded/mforest/qromb.c (revision 1099)
1
#include <math.h>
2
#define EPS 1.0e-6
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[j-K],&s[j-K],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/mforest-nothreads.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++)
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff