From 197e88c3b56d28840217010a2871c6ea3a4dd1a4 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 21 May 2014 10:30:55 +0200 Subject: update dependencies, move examples etc --- packages/hmatrix/src/Numeric/GSL/gsl-ode.c | 182 ----------------------------- 1 file changed, 182 deletions(-) delete mode 100644 packages/hmatrix/src/Numeric/GSL/gsl-ode.c (limited to 'packages/hmatrix/src/Numeric/GSL/gsl-ode.c') diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c b/packages/hmatrix/src/Numeric/GSL/gsl-ode.c deleted file mode 100644 index 3f2771b..0000000 --- a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c +++ /dev/null @@ -1,182 +0,0 @@ - -#ifdef GSLODE1 - -////////////////////////////// ODE V1 ////////////////////////////////////////// - -#include - -typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; - -int odefunc (double t, const double y[], double f[], void *params) { - Tode * P = (Tode*) params; - (P->f)(t,P->n,y,P->n,f); - return GSL_SUCCESS; -} - -int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { - Tode * P = ((Tode*) params); - (P->j)(t,P->n,y,P->n,P->n,dfdy); - int j; - for (j=0; j< P->n; j++) - dfdt[j] = 0.0; - return GSL_SUCCESS; -} - - -int ode(int method, double h, double eps_abs, double eps_rel, - int f(double, int, const double*, int, double*), - int jac(double, int, const double*, int, int, double*), - KRVEC(xi), KRVEC(ts), RMAT(sol)) { - - const gsl_odeiv_step_type * T; - - switch(method) { - case 0 : {T = gsl_odeiv_step_rk2; break; } - case 1 : {T = gsl_odeiv_step_rk4; break; } - case 2 : {T = gsl_odeiv_step_rkf45; break; } - case 3 : {T = gsl_odeiv_step_rkck; break; } - case 4 : {T = gsl_odeiv_step_rk8pd; break; } - case 5 : {T = gsl_odeiv_step_rk2imp; break; } - case 6 : {T = gsl_odeiv_step_rk4imp; break; } - case 7 : {T = gsl_odeiv_step_bsimp; break; } - case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); } - case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); } - case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); } - default: ERROR(BAD_CODE); - } - - gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); - gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); - gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); - - Tode P; - P.f = f; - P.j = jac; - P.n = xin; - - gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; - - double t = tsp[0]; - - double* y = (double*)calloc(xin,sizeof(double)); - int i,j; - for(i=0; i< xin; i++) { - y[i] = xip[i]; - solp[i] = xip[i]; - } - - for (i = 1; i < tsn ; i++) - { - double ti = tsp[i]; - while (t < ti) - { - gsl_odeiv_evolve_apply (e, c, s, - &sys, - &t, ti, &h, - y); - // if (h < hmin) h = hmin; - } - for(j=0; j - -typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; - -int odefunc (double t, const double y[], double f[], void *params) { - Tode * P = (Tode*) params; - (P->f)(t,P->n,y,P->n,f); - return GSL_SUCCESS; -} - -int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { - Tode * P = ((Tode*) params); - (P->j)(t,P->n,y,P->n,P->n,dfdy); - int j; - for (j=0; j< P->n; j++) - dfdt[j] = 0.0; - return GSL_SUCCESS; -} - - -int ode(int method, double h, double eps_abs, double eps_rel, - int f(double, int, const double*, int, double*), - int jac(double, int, const double*, int, int, double*), - KRVEC(xi), KRVEC(ts), RMAT(sol)) { - - const gsl_odeiv2_step_type * T; - - switch(method) { - case 0 : {T = gsl_odeiv2_step_rk2; break; } - case 1 : {T = gsl_odeiv2_step_rk4; break; } - case 2 : {T = gsl_odeiv2_step_rkf45; break; } - case 3 : {T = gsl_odeiv2_step_rkck; break; } - case 4 : {T = gsl_odeiv2_step_rk8pd; break; } - case 5 : {T = gsl_odeiv2_step_rk2imp; break; } - case 6 : {T = gsl_odeiv2_step_rk4imp; break; } - case 7 : {T = gsl_odeiv2_step_bsimp; break; } - case 8 : {T = gsl_odeiv2_step_rk1imp; break; } - case 9 : {T = gsl_odeiv2_step_msadams; break; } - case 10: {T = gsl_odeiv2_step_msbdf; break; } - default: ERROR(BAD_CODE); - } - - Tode P; - P.f = f; - P.j = jac; - P.n = xin; - - gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; - - gsl_odeiv2_driver * d = - gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); - - double t = tsp[0]; - - double* y = (double*)calloc(xin,sizeof(double)); - int i,j; - int status=0; - for(i=0; i< xin; i++) { - y[i] = xip[i]; - solp[i] = xip[i]; - } - - for (i = 1; i < tsn ; i++) - { - double ti = tsp[i]; - - status = gsl_odeiv2_driver_apply (d, &t, ti, y); - - if (status != GSL_SUCCESS) { - printf ("error in ode, return value=%d\n", status); - break; - } - -// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); - - for(j=0; j