Annotation of OpenXM/src/ox_gsl/call_gsl.c, Revision 1.8
1.8 ! takayama 1: /* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.7 2018/06/07 01:53:33 takayama Exp $
1.2 takayama 2: */
1.1 takayama 3: //#include <gsl/gsl_types.h>
4: //#include <gsl/gsl_sys.h>
1.2 takayama 5: #include <unistd.h>
1.1 takayama 6: #include <gsl/gsl_sf_result.h>
1.2 takayama 7: #include <gsl/gsl_errno.h>
1.1 takayama 8: #include <gsl/gsl_sf_gamma.h>
1.3 takayama 9: #include <gsl/gsl_integration.h>
1.5 takayama 10: #include <gsl/gsl_monte.h>
11: #include <gsl/gsl_monte_plain.h>
12: #include <gsl/gsl_monte_miser.h>
13: #include <gsl/gsl_monte_vegas.h>
1.1 takayama 14: #include "ox_gsl.h"
15: extern int Debug;
1.2 takayama 16: // local prototype declarations
17:
1.1 takayama 18: void call_gsl_sf_lngamma_complex_e() {
1.3 takayama 19: int argc;
1.1 takayama 20: double zr;
21: double zi;
22: gsl_sf_result lnr;
23: gsl_sf_result arg;
24: int status;
25: cmo *r[3];
26: cmo *ans;
1.2 takayama 27: // gsl_set_error_handler_off();
28: gsl_set_error_handler((gsl_error_handler_t *)myhandler);
1.3 takayama 29: argc = get_i(); // number of args
30: if (argc != 2) {
31: pops(argc);
32: push(make_error2("The argc must be 2 for gsl_sf_lngamma_complex_e.",NULL,0,-1));
33: return;
34: }
1.1 takayama 35: zr=get_double();
36: zi=get_double();
37: if (Debug) printf("gsl_sf_lngamma_complex_e(zr=%lg,zi=%lg)\n",zr,zi);
38: status = gsl_sf_lngamma_complex_e(zr,zi,&lnr,&arg);
39: r[0] = (cmo *)new_cmo_double(lnr.val);
40: r[1] = (cmo *)new_cmo_double(arg.val);
41: r[2] = (cmo *)new_cmo_int32(status);
1.2 takayama 42: ans = (cmo *)new_cmo_list_array((void *)r,3);
1.1 takayama 43: push(ans);
44: }
1.3 takayama 45:
1.7 takayama 46: cmo *Func_x=NULL; // function to evalute.
47: cmo *Vlist=NULL; // list of variables (not implemented yet).
1.3 takayama 48: double f_x(double x,void *params) {
49: double d;
1.4 takayama 50: if (Debug) ox_printf("f_x\n");
1.3 takayama 51: replace(1,"x",x);
1.4 takayama 52: if (Debug) ox_printf("f_x after replace x=%lg\n",x);
1.3 takayama 53: if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_x",GSL_ETOL);
1.4 takayama 54: if (Debug) ox_printf("f_x(%lg) -> d=%lg\n",x,d);
1.3 takayama 55: return(d);
56: }
57: void call_gsl_integration_qags() {
58: int argc;
59: double a;
60: double b;
61: int status;
62: cmo *r[3];
63: cmo *ans;
64: gsl_integration_workspace * w
65: = gsl_integration_workspace_alloc (1000);
66: double result, error;
67: gsl_function F;
68: double epsabs=0;
69: double epsrel=1e-7;
70: int limit=1000;
71:
72: // gsl_set_error_handler_off();
73: gsl_set_error_handler((gsl_error_handler_t *)myhandler);
74: argc = get_i(); // number of args
75: if (argc != 3) {
76: pops(argc);
77: push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1));
78: return;
79: }
80: Func_x = pop();
81: a = get_double();
82: b = get_double();
83:
84: F.function = &f_x;
1.4 takayama 85: F.params=NULL;
1.3 takayama 86:
87: status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,
88: w, &result, &error);
89:
1.4 takayama 90: if (Debug) ox_printf ("result = % .18f\n", result);
1.3 takayama 91: // printf ("estimated error = % .18f\n", error);
92: // printf ("intervals = %zu\n", w->size);
93:
94: gsl_integration_workspace_free(w);
95:
96: r[0] = (cmo *)new_cmo_double(result);
97: r[1] = (cmo *)new_cmo_double(error);
98: r[2] = (cmo *)new_cmo_int32(status);
99: ans = (cmo *)new_cmo_list_array((void *)r,3);
100: push(ans);
101: }
1.5 takayama 102:
1.8 ! takayama 103: #define VEC_MAX 30
1.5 takayama 104: double f_n(double x[],size_t dim,void *params) {
1.7 takayama 105: int debug_f_n=0;
1.5 takayama 106: double d;
107: int i;
1.8 ! takayama 108: char *xx[VEC_MAX]={"x0","x1","x2","x3","x4","x5","x6","x7","x8","x9",
1.5 takayama 109: "x10","x11","x12","x13","x14","x15","x16","x17","x18","x19",
110: "x20","x21","x22","x23","x24","x25","x26","x27","x28","x29"
111: };
1.7 takayama 112: // debug_f_n = Debug;
1.8 ! takayama 113: if (dim > VEC_MAX) GSL_ERROR("f_n supports functions with args <= VEC_MAX",GSL_ETOL);
1.6 takayama 114: init_dic();
1.7 takayama 115: if (debug_f_n) ox_printf("f_n, dim=%d\n",dim);
1.5 takayama 116: for (i=0; i<dim; i++) {
1.7 takayama 117: if (debug_f_n) ox_printf("x%d=%lg, ",i,x[i]);
1.6 takayama 118: register_entry(xx[i],x[i]);
1.5 takayama 119: }
1.7 takayama 120: if (debug_f_n) { ox_printf("\n"); fflush(NULL); }
1.5 takayama 121: if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_n",GSL_ETOL);
1.7 takayama 122: if (debug_f_n) ox_printf("\nf_x(...) -> d=%lg\n",d);
1.5 takayama 123: return(d);
124: }
1.7 takayama 125: void call_gsl_monte_plain_miser_vegas_integrate(int type) {
1.5 takayama 126: int argc;
127: int dim;
128: int dim2;
129: cmo *cr[2];
130: cmo *ans;
131: double *xl;
132: double *xu;
133: gsl_monte_function G = {&f_n,0,0};
134:
135: double res, err;
136: const gsl_rng_type *T;
137: gsl_rng *r;
1.6 takayama 138: // size_t calls = 500000;
1.7 takayama 139: size_t calls = 100000; // for test 1000
140: gsl_monte_plain_state *s_p;
141: gsl_monte_miser_state *s_m;
142: gsl_monte_vegas_state *s_v;
1.5 takayama 143: gsl_rng_env_setup ();
144: T = gsl_rng_default;
145: r = gsl_rng_alloc (T);
146:
147: // gsl_set_error_handler_off();
148: gsl_set_error_handler((gsl_error_handler_t *)myhandler);
149: argc = get_i(); // number of args
1.7 takayama 150: if (argc < 3) {
1.5 takayama 151: pops(argc);
1.7 takayama 152: push(make_error2("The argc must be >=3 for gsl_integration_qags.",NULL,0,-1));
1.5 takayama 153: return;
154: }
155: Func_x = pop();
156: xl = get_double_list(&dim);
157: xu = get_double_list(&dim2);
1.7 takayama 158: if (argc >= 4) {
159: calls = (int) get_double();
160: }
161: if (argc >= 5) {
162: Vlist = pop(); // list of variables. Not implemented yet.
163: }
1.5 takayama 164: if (dim != dim2) {
165: push(make_error2("gsl_monte_plain: dim of interval differs.",NULL,0,-1));
166: return;
167: }
1.6 takayama 168: G.dim=dim;
1.7 takayama 169: switch(type) {
170: case 0: // plain
171: s_p = gsl_monte_plain_alloc (dim);
172: gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s_p,
173: &res, &err);
174: gsl_monte_plain_free (s_p);
175: break;
176: case 1: // miser
177: s_m = gsl_monte_miser_alloc (dim);
178: gsl_monte_miser_integrate (&G, xl, xu, dim, calls, r, s_m,
179: &res, &err);
180: gsl_monte_miser_free (s_m);
181: break;
182: default: // vegas;
183: s_v = gsl_monte_vegas_alloc (dim);
184: gsl_monte_vegas_integrate (&G, xl, xu, dim, calls, r, s_v,
185: &res, &err);
186: gsl_monte_vegas_free (s_v);
187: break;
188: }
1.5 takayama 189: if (Debug) ox_printf("result = %lg\n",res);
190: cr[0] = (cmo *)new_cmo_double(res);
191: cr[1] = (cmo *)new_cmo_double(err);
192: ans = (cmo *)new_cmo_list_array((void *)cr,2);
193: push(ans);
194: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to call_gsl.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / ox_gsl