summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-08 08:48:12 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-08 08:48:12 +0200
commit1925c123d7d8184a1d2ddc0a413e0fd2776e1083 (patch)
treefad79f909d9c3be53d68e6ebd67202650536d387 /packages/hmatrix/src/Numeric/GSL/gsl-ode.c
parenteb3f702d065a4a967bb754977233e6eec408fd1f (diff)
empty hmatrix-base
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/gsl-ode.c')
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-ode.c182
1 files changed, 182 insertions, 0 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c b/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
new file mode 100644
index 0000000..3f2771b
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
@@ -0,0 +1,182 @@
1
2#ifdef GSLODE1
3
4////////////////////////////// ODE V1 //////////////////////////////////////////
5
6#include <gsl/gsl_odeiv.h>
7
8typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
9
10int odefunc (double t, const double y[], double f[], void *params) {
11 Tode * P = (Tode*) params;
12 (P->f)(t,P->n,y,P->n,f);
13 return GSL_SUCCESS;
14}
15
16int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
17 Tode * P = ((Tode*) params);
18 (P->j)(t,P->n,y,P->n,P->n,dfdy);
19 int j;
20 for (j=0; j< P->n; j++)
21 dfdt[j] = 0.0;
22 return GSL_SUCCESS;
23}
24
25
26int ode(int method, double h, double eps_abs, double eps_rel,
27 int f(double, int, const double*, int, double*),
28 int jac(double, int, const double*, int, int, double*),
29 KRVEC(xi), KRVEC(ts), RMAT(sol)) {
30
31 const gsl_odeiv_step_type * T;
32
33 switch(method) {
34 case 0 : {T = gsl_odeiv_step_rk2; break; }
35 case 1 : {T = gsl_odeiv_step_rk4; break; }
36 case 2 : {T = gsl_odeiv_step_rkf45; break; }
37 case 3 : {T = gsl_odeiv_step_rkck; break; }
38 case 4 : {T = gsl_odeiv_step_rk8pd; break; }
39 case 5 : {T = gsl_odeiv_step_rk2imp; break; }
40 case 6 : {T = gsl_odeiv_step_rk4imp; break; }
41 case 7 : {T = gsl_odeiv_step_bsimp; break; }
42 case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); }
43 case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); }
44 case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); }
45 default: ERROR(BAD_CODE);
46 }
47
48 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
49 gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
50 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
51
52 Tode P;
53 P.f = f;
54 P.j = jac;
55 P.n = xin;
56
57 gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
58
59 double t = tsp[0];
60
61 double* y = (double*)calloc(xin,sizeof(double));
62 int i,j;
63 for(i=0; i< xin; i++) {
64 y[i] = xip[i];
65 solp[i] = xip[i];
66 }
67
68 for (i = 1; i < tsn ; i++)
69 {
70 double ti = tsp[i];
71 while (t < ti)
72 {
73 gsl_odeiv_evolve_apply (e, c, s,
74 &sys,
75 &t, ti, &h,
76 y);
77 // if (h < hmin) h = hmin;
78 }
79 for(j=0; j<xin; j++) {
80 solp[i*xin + j] = y[j];
81 }
82 }
83
84 free(y);
85 gsl_odeiv_evolve_free (e);
86 gsl_odeiv_control_free (c);
87 gsl_odeiv_step_free (s);
88 return 0;
89}
90
91#else
92
93///////////////////// ODE V2 ///////////////////////////////////////////////////
94
95#include <gsl/gsl_odeiv2.h>
96
97typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
98
99int odefunc (double t, const double y[], double f[], void *params) {
100 Tode * P = (Tode*) params;
101 (P->f)(t,P->n,y,P->n,f);
102 return GSL_SUCCESS;
103}
104
105int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
106 Tode * P = ((Tode*) params);
107 (P->j)(t,P->n,y,P->n,P->n,dfdy);
108 int j;
109 for (j=0; j< P->n; j++)
110 dfdt[j] = 0.0;
111 return GSL_SUCCESS;
112}
113
114
115int ode(int method, double h, double eps_abs, double eps_rel,
116 int f(double, int, const double*, int, double*),
117 int jac(double, int, const double*, int, int, double*),
118 KRVEC(xi), KRVEC(ts), RMAT(sol)) {
119
120 const gsl_odeiv2_step_type * T;
121
122 switch(method) {
123 case 0 : {T = gsl_odeiv2_step_rk2; break; }
124 case 1 : {T = gsl_odeiv2_step_rk4; break; }
125 case 2 : {T = gsl_odeiv2_step_rkf45; break; }
126 case 3 : {T = gsl_odeiv2_step_rkck; break; }
127 case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
128 case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
129 case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
130 case 7 : {T = gsl_odeiv2_step_bsimp; break; }
131 case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
132 case 9 : {T = gsl_odeiv2_step_msadams; break; }
133 case 10: {T = gsl_odeiv2_step_msbdf; break; }
134 default: ERROR(BAD_CODE);
135 }
136
137 Tode P;
138 P.f = f;
139 P.j = jac;
140 P.n = xin;
141
142 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
143
144 gsl_odeiv2_driver * d =
145 gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
146
147 double t = tsp[0];
148
149 double* y = (double*)calloc(xin,sizeof(double));
150 int i,j;
151 int status=0;
152 for(i=0; i< xin; i++) {
153 y[i] = xip[i];
154 solp[i] = xip[i];
155 }
156
157 for (i = 1; i < tsn ; i++)
158 {
159 double ti = tsp[i];
160
161 status = gsl_odeiv2_driver_apply (d, &t, ti, y);
162
163 if (status != GSL_SUCCESS) {
164 printf ("error in ode, return value=%d\n", status);
165 break;
166 }
167
168// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
169
170 for(j=0; j<xin; j++) {
171 solp[i*xin + j] = y[j];
172 }
173 }
174
175 free(y);
176 gsl_odeiv2_driver_free (d);
177
178 return status;
179}
180
181#endif
182