version 1.5, 2018/06/06 07:40:32 |
version 1.7, 2018/06/07 01:53:33 |
|
|
/* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.4 2018/04/18 02:20:51 takayama Exp $ |
/* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.6 2018/06/06 10:46:00 takayama Exp $ |
*/ |
*/ |
//#include <gsl/gsl_types.h> |
//#include <gsl/gsl_types.h> |
//#include <gsl/gsl_sys.h> |
//#include <gsl/gsl_sys.h> |
Line 43 void call_gsl_sf_lngamma_complex_e() { |
|
Line 43 void call_gsl_sf_lngamma_complex_e() { |
|
push(ans); |
push(ans); |
} |
} |
|
|
cmo *Func_x=NULL; |
cmo *Func_x=NULL; // function to evalute. |
|
cmo *Vlist=NULL; // list of variables (not implemented yet). |
double f_x(double x,void *params) { |
double f_x(double x,void *params) { |
double d; |
double d; |
if (Debug) ox_printf("f_x\n"); |
if (Debug) ox_printf("f_x\n"); |
Line 100 void call_gsl_integration_qags() { |
|
Line 101 void call_gsl_integration_qags() { |
|
} |
} |
|
|
double f_n(double x[],size_t dim,void *params) { |
double f_n(double x[],size_t dim,void *params) { |
|
int debug_f_n=0; |
double d; |
double d; |
int i; |
int i; |
char *xx[30]={"x0","x1","x2","x3","x4","x5","x6","x7","x8","x9", |
char *xx[30]={"x0","x1","x2","x3","x4","x5","x6","x7","x8","x9", |
"x10","x11","x12","x13","x14","x15","x16","x17","x18","x19", |
"x10","x11","x12","x13","x14","x15","x16","x17","x18","x19", |
"x20","x21","x22","x23","x24","x25","x26","x27","x28","x29" |
"x20","x21","x22","x23","x24","x25","x26","x27","x28","x29" |
}; |
}; |
|
// debug_f_n = Debug; |
if (dim > 30) GSL_ERROR("f_n supports functions with args <= 30",GSL_ETOL); |
if (dim > 30) GSL_ERROR("f_n supports functions with args <= 30",GSL_ETOL); |
if (Debug) ox_printf("f_n\n"); |
init_dic(); |
|
if (debug_f_n) ox_printf("f_n, dim=%d\n",dim); |
for (i=0; i<dim; i++) { |
for (i=0; i<dim; i++) { |
replace(1,xx[i],x[i]); |
if (debug_f_n) ox_printf("x%d=%lg, ",i,x[i]); |
|
register_entry(xx[i],x[i]); |
} |
} |
|
if (debug_f_n) { ox_printf("\n"); fflush(NULL); } |
if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_n",GSL_ETOL); |
if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_n",GSL_ETOL); |
if (Debug) ox_printf("f_x(...) -> d=%lg\n",d); |
if (debug_f_n) ox_printf("\nf_x(...) -> d=%lg\n",d); |
return(d); |
return(d); |
} |
} |
void call_gsl_monte_plain_integrate() { |
void call_gsl_monte_plain_miser_vegas_integrate(int type) { |
int argc; |
int argc; |
int dim; |
int dim; |
int dim2; |
int dim2; |
Line 128 void call_gsl_monte_plain_integrate() { |
|
Line 134 void call_gsl_monte_plain_integrate() { |
|
double res, err; |
double res, err; |
const gsl_rng_type *T; |
const gsl_rng_type *T; |
gsl_rng *r; |
gsl_rng *r; |
size_t calls = 500000; |
// size_t calls = 500000; |
|
size_t calls = 100000; // for test 1000 |
|
gsl_monte_plain_state *s_p; |
|
gsl_monte_miser_state *s_m; |
|
gsl_monte_vegas_state *s_v; |
gsl_rng_env_setup (); |
gsl_rng_env_setup (); |
T = gsl_rng_default; |
T = gsl_rng_default; |
r = gsl_rng_alloc (T); |
r = gsl_rng_alloc (T); |
Line 136 void call_gsl_monte_plain_integrate() { |
|
Line 146 void call_gsl_monte_plain_integrate() { |
|
// gsl_set_error_handler_off(); |
// gsl_set_error_handler_off(); |
gsl_set_error_handler((gsl_error_handler_t *)myhandler); |
gsl_set_error_handler((gsl_error_handler_t *)myhandler); |
argc = get_i(); // number of args |
argc = get_i(); // number of args |
if (argc != 3) { |
if (argc < 3) { |
pops(argc); |
pops(argc); |
push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1)); |
push(make_error2("The argc must be >=3 for gsl_integration_qags.",NULL,0,-1)); |
return; |
return; |
} |
} |
Func_x = pop(); |
Func_x = pop(); |
xl = get_double_list(&dim); |
xl = get_double_list(&dim); |
xu = get_double_list(&dim2); |
xu = get_double_list(&dim2); |
|
if (argc >= 4) { |
|
calls = (int) get_double(); |
|
} |
|
if (argc >= 5) { |
|
Vlist = pop(); // list of variables. Not implemented yet. |
|
} |
if (dim != dim2) { |
if (dim != dim2) { |
push(make_error2("gsl_monte_plain: dim of interval differs.",NULL,0,-1)); |
push(make_error2("gsl_monte_plain: dim of interval differs.",NULL,0,-1)); |
return; |
return; |
} |
} |
gsl_monte_plain_state *s = gsl_monte_plain_alloc (dim); |
G.dim=dim; |
gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s, |
switch(type) { |
&res, &err); |
case 0: // plain |
gsl_monte_plain_free (s); |
s_p = gsl_monte_plain_alloc (dim); |
|
gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s_p, |
|
&res, &err); |
|
gsl_monte_plain_free (s_p); |
|
break; |
|
case 1: // miser |
|
s_m = gsl_monte_miser_alloc (dim); |
|
gsl_monte_miser_integrate (&G, xl, xu, dim, calls, r, s_m, |
|
&res, &err); |
|
gsl_monte_miser_free (s_m); |
|
break; |
|
default: // vegas; |
|
s_v = gsl_monte_vegas_alloc (dim); |
|
gsl_monte_vegas_integrate (&G, xl, xu, dim, calls, r, s_v, |
|
&res, &err); |
|
gsl_monte_vegas_free (s_v); |
|
break; |
|
} |
if (Debug) ox_printf("result = %lg\n",res); |
if (Debug) ox_printf("result = %lg\n",res); |
cr[0] = (cmo *)new_cmo_double(res); |
cr[0] = (cmo *)new_cmo_double(res); |
cr[1] = (cmo *)new_cmo_double(err); |
cr[1] = (cmo *)new_cmo_double(err); |