[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.3 and 1.5

version 1.3, 2018/04/17 00:56:38 version 1.5, 2018/06/06 07:40:32
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.2 2018/03/30 04:43:16 takayama Exp $  /* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.4 2018/04/18 02:20:51 takayama Exp $
 */  */
 //#include <gsl/gsl_types.h>  //#include <gsl/gsl_types.h>
 //#include <gsl/gsl_sys.h>  //#include <gsl/gsl_sys.h>
Line 7 
Line 7 
 #include <gsl/gsl_errno.h>  #include <gsl/gsl_errno.h>
 #include <gsl/gsl_sf_gamma.h>  #include <gsl/gsl_sf_gamma.h>
 #include <gsl/gsl_integration.h>  #include <gsl/gsl_integration.h>
   #include <gsl/gsl_monte.h>
   #include <gsl/gsl_monte_plain.h>
   #include <gsl/gsl_monte_miser.h>
   #include <gsl/gsl_monte_vegas.h>
 #include "ox_gsl.h"  #include "ox_gsl.h"
 extern int Debug;  extern int Debug;
 // local prototype declarations  // local prototype declarations
Line 42  void  call_gsl_sf_lngamma_complex_e() {
Line 46  void  call_gsl_sf_lngamma_complex_e() {
 cmo *Func_x=NULL;  cmo *Func_x=NULL;
 double f_x(double x,void *params) {  double f_x(double x,void *params) {
   double d;    double d;
     if (Debug) ox_printf("f_x\n");
   replace(1,"x",x);    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 (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);    return(d);
 }  }
 void  call_gsl_integration_qags() {  void  call_gsl_integration_qags() {
Line 74  void  call_gsl_integration_qags() {
Line 81  void  call_gsl_integration_qags() {
   b = get_double();    b = get_double();
   
   F.function = &f_x;    F.function = &f_x;
     F.params=NULL;
   
   status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,    status=gsl_integration_qags (&F, a, b, epsabs, epsrel, limit,
                                w, &result, &error);                                 w, &result, &error);
   
 //  printf ("result          = % .18f\n", result);    if (Debug) ox_printf ("result          = % .18f\n", result);
 //  printf ("estimated error = % .18f\n", error);  //  printf ("estimated error = % .18f\n", error);
 //  printf ("intervals       = %zu\n", w->size);  //  printf ("intervals       = %zu\n", w->size);
   
Line 90  void  call_gsl_integration_qags() {
Line 98  void  call_gsl_integration_qags() {
   ans = (cmo *)new_cmo_list_array((void *)r,3);    ans = (cmo *)new_cmo_list_array((void *)r,3);
   push(ans);    push(ans);
 }  }
   
   double f_n(double x[],size_t dim,void *params) {
     double d;
     int i;
     char *xx[30]={"x0","x1","x2","x3","x4","x5","x6","x7","x8","x9",
                   "x10","x11","x12","x13","x14","x15","x16","x17","x18","x19",
                   "x20","x21","x22","x23","x24","x25","x26","x27","x28","x29"
                  };
     if (dim > 30) GSL_ERROR("f_n supports functions with args <= 30",GSL_ETOL);
     if (Debug) ox_printf("f_n\n");
     for (i=0; i<dim; i++) {
       replace(1,xx[i],x[i]);
     }
     if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_n",GSL_ETOL);
     if (Debug) ox_printf("f_x(...) -> d=%lg\n",d);
     return(d);
   }
   void  call_gsl_monte_plain_integrate() {
     int argc;
     int dim;
     int dim2;
     cmo *cr[2];
     cmo *ans;
     double *xl;
     double *xu;
     gsl_monte_function G = {&f_n,0,0};
   
     double res, err;
     const gsl_rng_type *T;
     gsl_rng *r;
     size_t calls = 500000;
     gsl_rng_env_setup ();
     T = gsl_rng_default;
     r = gsl_rng_alloc (T);
   
     //  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();
     xl = get_double_list(&dim);
     xu = get_double_list(&dim2);
     if (dim != dim2) {
       push(make_error2("gsl_monte_plain: dim of interval differs.",NULL,0,-1));
       return;
     }
     gsl_monte_plain_state *s = gsl_monte_plain_alloc (dim);
     gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s,
                                &res, &err);
     gsl_monte_plain_free (s);
   
     if (Debug) ox_printf("result = %lg\n",res);
     cr[0] = (cmo *)new_cmo_double(res);
     cr[1] = (cmo *)new_cmo_double(err);
     ans = (cmo *)new_cmo_list_array((void *)cr,2);
     push(ans);
   }
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>