File: [local] / OpenXM / src / ox_gsl / call_gsl.c (download)
Revision 1.3, Tue Apr 17 00:56:38 2018 UTC (6 years, 2 months ago) by takayama
Branch: MAIN
Changes since 1.2: +61 -3
lines
gsl_integration_qags (adaptive numerical integration) is added.
ox_cmo_rpc(Pid,"gsl_integration_qags", f(x),a,b) returns an approximate value of
int_a^b f(x) dx.
Example:
ox_launch(0,getenv("OpenXM_HOME")+"/src/ox_gsl/ox_gsl");
0
[1817] ox_cmo_rpc(0,"gsl_integration_qags",quote(x^2*log(x)),0.0,1.0);
0
[1818] ox_pop_cmo(0);
[-0.111111,9.27977e-09,0]
|
/* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.3 2018/04/17 00:56:38 takayama Exp $
*/
//#include <gsl/gsl_types.h>
//#include <gsl/gsl_sys.h>
#include <unistd.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_integration.h>
#include "ox_gsl.h"
extern int Debug;
// local prototype declarations
void call_gsl_sf_lngamma_complex_e() {
int argc;
double zr;
double zi;
gsl_sf_result lnr;
gsl_sf_result arg;
int status;
cmo *r[3];
cmo *ans;
// gsl_set_error_handler_off();
gsl_set_error_handler((gsl_error_handler_t *)myhandler);
argc = get_i(); // number of args
if (argc != 2) {
pops(argc);
push(make_error2("The argc must be 2 for gsl_sf_lngamma_complex_e.",NULL,0,-1));
return;
}
zr=get_double();
zi=get_double();
if (Debug) printf("gsl_sf_lngamma_complex_e(zr=%lg,zi=%lg)\n",zr,zi);
status = gsl_sf_lngamma_complex_e(zr,zi,&lnr,&arg);
r[0] = (cmo *)new_cmo_double(lnr.val);
r[1] = (cmo *)new_cmo_double(arg.val);
r[2] = (cmo *)new_cmo_int32(status);
ans = (cmo *)new_cmo_list_array((void *)r,3);
push(ans);
}
cmo *Func_x=NULL;
double f_x(double x,void *params) {
double d;
replace(1,"x",x);
if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_x",GSL_ETOL);
return(d);
}
void call_gsl_integration_qags() {
int argc;
double a;
double b;
int status;
cmo *r[3];
cmo *ans;
gsl_integration_workspace * w
= gsl_integration_workspace_alloc (1000);
double result, error;
gsl_function F;
double epsabs=0;
double epsrel=1e-7;
int limit=1000;
// gsl_set_error_handler_off();
gsl_set_error_handler((gsl_error_handler_t *)myhandler);
argc = get_i(); // number of args
if (argc != 3) {
pops(argc);
push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1));
return;
}
Func_x = pop();
a = get_double();
b = get_double();
F.function = &f_x;
status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,
w, &result, &error);
// printf ("result = % .18f\n", result);
// printf ("estimated error = % .18f\n", error);
// printf ("intervals = %zu\n", w->size);
gsl_integration_workspace_free(w);
r[0] = (cmo *)new_cmo_double(result);
r[1] = (cmo *)new_cmo_double(error);
r[2] = (cmo *)new_cmo_int32(status);
ans = (cmo *)new_cmo_list_array((void *)r,3);
push(ans);
}