[BACK]Return to call_gsl.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / ox_gsl

Diff for /OpenXM/src/ox_gsl/call_gsl.c between version 1.1 and 1.4

version 1.1, 2018/03/29 11:52:18 version 1.4, 2018/04/18 02:20:51
Line 1 
Line 1 
   /* $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_types.h>
 //#include <gsl/gsl_sys.h>  //#include <gsl/gsl_sys.h>
   #include <unistd.h>
 #include <gsl/gsl_sf_result.h>  #include <gsl/gsl_sf_result.h>
   #include <gsl/gsl_errno.h>
 #include <gsl/gsl_sf_gamma.h>  #include <gsl/gsl_sf_gamma.h>
   #include <gsl/gsl_integration.h>
 #include "ox_gsl.h"  #include "ox_gsl.h"
 extern int Debug;  extern int Debug;
   // local prototype declarations
   
 void  call_gsl_sf_lngamma_complex_e() {  void  call_gsl_sf_lngamma_complex_e() {
   cmo *c;    int argc;
   double zr;    double zr;
   double zi;    double zi;
   gsl_sf_result lnr;    gsl_sf_result lnr;
Line 13  void  call_gsl_sf_lngamma_complex_e() {
Line 20  void  call_gsl_sf_lngamma_complex_e() {
   int status;    int status;
   cmo *r[3];    cmo *r[3];
   cmo *ans;    cmo *ans;
   gsl_set_error_handler_off();    //  gsl_set_error_handler_off();
   // Todo, gsl_set_error_handler(my_handler);    gsl_set_error_handler((gsl_error_handler_t *)myhandler);
   c = pop(); // number of args    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();    zr=get_double();
   zi=get_double();    zi=get_double();
   if (Debug) printf("gsl_sf_lngamma_complex_e(zr=%lg,zi=%lg)\n",zr,zi);    if (Debug) printf("gsl_sf_lngamma_complex_e(zr=%lg,zi=%lg)\n",zr,zi);
Line 23  void  call_gsl_sf_lngamma_complex_e() {
Line 35  void  call_gsl_sf_lngamma_complex_e() {
   r[0] = (cmo *)new_cmo_double(lnr.val);    r[0] = (cmo *)new_cmo_double(lnr.val);
   r[1] = (cmo *)new_cmo_double(arg.val);    r[1] = (cmo *)new_cmo_double(arg.val);
   r[2] = (cmo *)new_cmo_int32(status);    r[2] = (cmo *)new_cmo_int32(status);
   ans = (cmo *)new_cmo_list_array((void **)r,3);    ans = (cmo *)new_cmo_list_array((void *)r,3);
   push(ans);    push(ans);
 }  }
   
   cmo *Func_x=NULL;
   double f_x(double x,void *params) {
     double d;
     if (Debug) ox_printf("f_x\n");
     replace(1,"x",x);
     if (Debug) ox_printf("f_x after replace x=%lg\n",x);
     if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_x",GSL_ETOL);
     if (Debug) ox_printf("f_x(%lg) -> d=%lg\n",x,d);
     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;
     F.params=NULL;
   
     status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,
                                  w, &result, &error);
   
     if (Debug) ox_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);
   }

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>