Annotation of OpenXM/src/ox_gsl/call_gsl.c, Revision 1.5
1.5 ! takayama 1: /* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.4 2018/04/18 02:20:51 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:
46: cmo *Func_x=NULL;
47: double f_x(double x,void *params) {
48: double d;
1.4 takayama 49: if (Debug) ox_printf("f_x\n");
1.3 takayama 50: replace(1,"x",x);
1.4 takayama 51: if (Debug) ox_printf("f_x after replace x=%lg\n",x);
1.3 takayama 52: if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_x",GSL_ETOL);
1.4 takayama 53: if (Debug) ox_printf("f_x(%lg) -> d=%lg\n",x,d);
1.3 takayama 54: return(d);
55: }
56: void call_gsl_integration_qags() {
57: int argc;
58: double a;
59: double b;
60: int status;
61: cmo *r[3];
62: cmo *ans;
63: gsl_integration_workspace * w
64: = gsl_integration_workspace_alloc (1000);
65: double result, error;
66: gsl_function F;
67: double epsabs=0;
68: double epsrel=1e-7;
69: int limit=1000;
70:
71: // gsl_set_error_handler_off();
72: gsl_set_error_handler((gsl_error_handler_t *)myhandler);
73: argc = get_i(); // number of args
74: if (argc != 3) {
75: pops(argc);
76: push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1));
77: return;
78: }
79: Func_x = pop();
80: a = get_double();
81: b = get_double();
82:
83: F.function = &f_x;
1.4 takayama 84: F.params=NULL;
1.3 takayama 85:
86: status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,
87: w, &result, &error);
88:
1.4 takayama 89: if (Debug) ox_printf ("result = % .18f\n", result);
1.3 takayama 90: // printf ("estimated error = % .18f\n", error);
91: // printf ("intervals = %zu\n", w->size);
92:
93: gsl_integration_workspace_free(w);
94:
95: r[0] = (cmo *)new_cmo_double(result);
96: r[1] = (cmo *)new_cmo_double(error);
97: r[2] = (cmo *)new_cmo_int32(status);
98: ans = (cmo *)new_cmo_list_array((void *)r,3);
99: push(ans);
100: }
1.5 ! takayama 101:
! 102: double f_n(double x[],size_t dim,void *params) {
! 103: double d;
! 104: int i;
! 105: char *xx[30]={"x0","x1","x2","x3","x4","x5","x6","x7","x8","x9",
! 106: "x10","x11","x12","x13","x14","x15","x16","x17","x18","x19",
! 107: "x20","x21","x22","x23","x24","x25","x26","x27","x28","x29"
! 108: };
! 109: if (dim > 30) GSL_ERROR("f_n supports functions with args <= 30",GSL_ETOL);
! 110: if (Debug) ox_printf("f_n\n");
! 111: for (i=0; i<dim; i++) {
! 112: replace(1,xx[i],x[i]);
! 113: }
! 114: if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_n",GSL_ETOL);
! 115: if (Debug) ox_printf("f_x(...) -> d=%lg\n",d);
! 116: return(d);
! 117: }
! 118: void call_gsl_monte_plain_integrate() {
! 119: int argc;
! 120: int dim;
! 121: int dim2;
! 122: cmo *cr[2];
! 123: cmo *ans;
! 124: double *xl;
! 125: double *xu;
! 126: gsl_monte_function G = {&f_n,0,0};
! 127:
! 128: double res, err;
! 129: const gsl_rng_type *T;
! 130: gsl_rng *r;
! 131: size_t calls = 500000;
! 132: gsl_rng_env_setup ();
! 133: T = gsl_rng_default;
! 134: r = gsl_rng_alloc (T);
! 135:
! 136: // gsl_set_error_handler_off();
! 137: gsl_set_error_handler((gsl_error_handler_t *)myhandler);
! 138: argc = get_i(); // number of args
! 139: if (argc != 3) {
! 140: pops(argc);
! 141: push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1));
! 142: return;
! 143: }
! 144: Func_x = pop();
! 145: xl = get_double_list(&dim);
! 146: xu = get_double_list(&dim2);
! 147: if (dim != dim2) {
! 148: push(make_error2("gsl_monte_plain: dim of interval differs.",NULL,0,-1));
! 149: return;
! 150: }
! 151: gsl_monte_plain_state *s = gsl_monte_plain_alloc (dim);
! 152: gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s,
! 153: &res, &err);
! 154: gsl_monte_plain_free (s);
! 155:
! 156: if (Debug) ox_printf("result = %lg\n",res);
! 157: cr[0] = (cmo *)new_cmo_double(res);
! 158: cr[1] = (cmo *)new_cmo_double(err);
! 159: ans = (cmo *)new_cmo_list_array((void *)cr,2);
! 160: push(ans);
! 161: }
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