Return to call_gsl.c CVS log

Up to [local] / OpenXM / src / ox_gsl

Annotation of OpenXM/src/ox_gsl/call_gsl.c, Revision 1.4

1.4     ! takayama    1: /* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.3 2018/04/17 00:56:38 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.1       takayama   10: #include "ox_gsl.h"
                     11: extern int Debug;
1.2       takayama   12: // local prototype declarations
                     13:
1.1       takayama   14: void  call_gsl_sf_lngamma_complex_e() {
1.3       takayama   15:   int argc;
1.1       takayama   16:   double zr;
                     17:   double zi;
                     18:   gsl_sf_result lnr;
                     19:   gsl_sf_result arg;
                     20:   int status;
                     21:   cmo *r[3];
                     22:   cmo *ans;
1.2       takayama   23:   //  gsl_set_error_handler_off();
                     24:   gsl_set_error_handler((gsl_error_handler_t *)myhandler);
1.3       takayama   25:   argc = get_i(); // number of args
                     26:   if (argc != 2) {
                     27:     pops(argc);
                     28:     push(make_error2("The argc must be 2 for gsl_sf_lngamma_complex_e.",NULL,0,-1));
                     29:     return;
                     30:   }
1.1       takayama   31:   zr=get_double();
                     32:   zi=get_double();
                     33:   if (Debug) printf("gsl_sf_lngamma_complex_e(zr=%lg,zi=%lg)\n",zr,zi);
                     34:   status = gsl_sf_lngamma_complex_e(zr,zi,&lnr,&arg);
                     35:   r[0] = (cmo *)new_cmo_double(lnr.val);
                     36:   r[1] = (cmo *)new_cmo_double(arg.val);
                     37:   r[2] = (cmo *)new_cmo_int32(status);
1.2       takayama   38:   ans = (cmo *)new_cmo_list_array((void *)r,3);
1.1       takayama   39:   push(ans);
                     40: }
1.3       takayama   41:
                     42: cmo *Func_x=NULL;
                     43: double f_x(double x,void *params) {
                     44:   double d;
1.4     ! takayama   45:   if (Debug) ox_printf("f_x\n");
1.3       takayama   46:   replace(1,"x",x);
1.4     ! takayama   47:   if (Debug) ox_printf("f_x after replace x=%lg\n",x);
1.3       takayama   48:   if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_x",GSL_ETOL);
1.4     ! takayama   49:   if (Debug) ox_printf("f_x(%lg) -> d=%lg\n",x,d);
1.3       takayama   50:   return(d);
                     51: }
                     52: void  call_gsl_integration_qags() {
                     53:   int argc;
                     54:   double a;
                     55:   double b;
                     56:   int status;
                     57:   cmo *r[3];
                     58:   cmo *ans;
                     59:   gsl_integration_workspace * w
                     60:     = gsl_integration_workspace_alloc (1000);
                     61:   double result, error;
                     62:   gsl_function F;
                     63:   double epsabs=0;
                     64:   double epsrel=1e-7;
                     65:   int limit=1000;
                     66:
                     67:   //  gsl_set_error_handler_off();
                     68:   gsl_set_error_handler((gsl_error_handler_t *)myhandler);
                     69:   argc = get_i(); // number of args
                     70:   if (argc != 3) {
                     71:     pops(argc);
                     72:     push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1));
                     73:     return;
                     74:   }
                     75:   Func_x = pop();
                     76:   a = get_double();
                     77:   b = get_double();
                     78:
                     79:   F.function = &f_x;
1.4     ! takayama   80:   F.params=NULL;
1.3       takayama   81:
                     82:   status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,
                     83:                                w, &result, &error);
                     84:
1.4     ! takayama   85:   if (Debug) ox_printf ("result          = % .18f\n", result);
1.3       takayama   86: //  printf ("estimated error = % .18f\n", error);
                     87: //  printf ("intervals       = %zu\n", w->size);
                     88:
                     89:   gsl_integration_workspace_free(w);
                     90:
                     91:   r[0] = (cmo *)new_cmo_double(result);
                     92:   r[1] = (cmo *)new_cmo_double(error);
                     93:   r[2] = (cmo *)new_cmo_int32(status);
                     94:   ans = (cmo *)new_cmo_list_array((void *)r,3);
                     95:   push(ans);
                     96: }