1 | /* ode-initval2/driver.c |
2 | * |
3 | * Copyright (C) 2009, 2010 Tuomo Keskitalo |
4 | * |
5 | * This program is free software; you can redistribute it and/or modify |
6 | * it under the terms of the GNU General Public License as published by |
7 | * the Free Software Foundation; either version 3 of the License, or (at |
8 | * your option) any later version. |
9 | * |
10 | * This program is distributed in the hope that it will be useful, but |
11 | * WITHOUT ANY WARRANTY; without even the implied warranty of |
12 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
13 | * General Public License for more details. |
14 | * |
15 | * You should have received a copy of the GNU General Public License |
16 | * along with this program; if not, write to the Free Software |
17 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
18 | */ |
19 | |
20 | /* Driver routine for odeiv2. This is a wrapper for low level GSL |
21 | functions that allows a simple interface to step, control and |
22 | evolve layers. |
23 | */ |
24 | |
25 | #include <config.h> |
26 | #include <math.h> |
27 | #include <gsl/gsl_errno.h> |
28 | #include <gsl/gsl_odeiv2.h> |
29 | #include <gsl/gsl_machine.h> |
30 | |
31 | static gsl_odeiv2_driver * |
32 | driver_alloc (const gsl_odeiv2_system * sys, const double hstart, |
| 1Start Analysis. |
33 | const gsl_odeiv2_step_type * T) |
34 | { |
35 | /* Allocates and initializes an ODE driver system. Step and evolve |
36 | objects are allocated here, but control object is allocated in |
37 | another function. |
38 | */ |
39 | |
40 | gsl_odeiv2_driver *state = |
41 | (gsl_odeiv2_driver *) malloc (sizeof (gsl_odeiv2_driver)); |
| 2Call a function.    malloc(sizeof(gsl_odeiv2_driver)) |
| 33Memory allocated in heap space is not freed.    malloc(sizeof(gsl_odeiv2_driver)) |
42 | |
43 | if (state == NULL) |
| 3Take the false branch.    state == ((void *)0) |
44 | { |
45 | GSL_ERROR_NULL ("failed to allocate space for driver state", |
46 | GSL_ENOMEM); |
47 | } |
48 | |
49 | if (sys == NULL) |
| 4Take the false branch.    sys == ((void *)0) |
50 | { |
51 | GSL_ERROR_NULL ("gsl_odeiv2_system must be defined", GSL_EINVAL); |
52 | } |
53 | |
54 | { |
55 | const size_t dim = sys->dimension; |
56 | |
57 | if (dim == 0) |
| 5Take the false branch.    dim == 0 |
58 | { |
59 | GSL_ERROR_NULL |
60 | ("gsl_odeiv2_system dimension must be a positive integer", |
61 | GSL_EINVAL); |
62 | } |
63 | |
64 | state->sys = sys; |
65 | |
66 | state->s = gsl_odeiv2_step_alloc (T, dim); |
| 6Call a function.    gsl_odeiv2_step_alloc(T, dim) |
67 | |
68 | if (state->s == NULL) |
| 12Take the false branch.    state->s == ((void *)0) |
69 | { |
70 | free (state); |
71 | GSL_ERROR_NULL ("failed to allocate step object", GSL_ENOMEM); |
72 | } |
73 | |
74 | state->e = gsl_odeiv2_evolve_alloc (dim); |
| 13Call a function.    gsl_odeiv2_evolve_alloc(dim) |
75 | } |
76 | |
77 | if (state->e == NULL) |
| 26Take the false branch.    state->e == ((void *)0) |
78 | { |
79 | gsl_odeiv2_step_free (state->s); |
80 | free (state); |
81 | GSL_ERROR_NULL ("failed to allocate evolve object", GSL_ENOMEM); |
82 | } |
83 | |
84 | if (hstart > 0.0 || hstart < 0.0) |
| 27Take the false branch.    hstart > 0. |
| 28Take the false branch.    hstart > 0. || hstart < 0. |
85 | { |
86 | state->h = hstart; |
87 | } |
88 | else |
89 | { |
90 | GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL); |
| 29Call a function.    gsl_error("invalid hstart", "/home/xingming/workplace/experiment/fp/gsl-2.1/ode-initval2/driver.c", 90, GSL_EINVAL) |
91 | } |
92 | |
93 | state->h = hstart; |
94 | state->hmin = 0.0; |
95 | state->hmax = GSL_DBL_MAX; |
96 | state->nmax = 0; |
97 | state->n = 0; |
98 | state->c = NULL; |
99 | |
100 | return state; |
101 | } |
102 | |
103 | int |
104 | gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin) |
105 | { |
106 | /* Sets minimum allowed step size fabs(hmin) for driver. It is |
107 | required that hmin <= fabs(h) <= hmax. */ |
108 | |
109 | if ((fabs (hmin) > fabs (d->h)) || (fabs (hmin) > d->hmax)) |
110 | { |
111 | GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL); |
112 | } |
113 | |
114 | d->hmin = fabs (hmin); |
115 | |
116 | return GSL_SUCCESS; |
117 | } |
118 | |
119 | int |
120 | gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax) |
121 | { |
122 | /* Sets maximum allowed step size fabs(hmax) for driver. It is |
123 | required that hmin <= fabs(h) <= hmax. */ |
124 | |
125 | if ((fabs (hmax) < fabs (d->h)) || (fabs (hmax) < d->hmin)) |
126 | { |
127 | GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL); |
128 | } |
129 | |
130 | if (hmax > 0.0 || hmax < 0.0) |
131 | { |
132 | d->hmax = fabs (hmax); |
133 | } |
134 | else |
135 | { |
136 | GSL_ERROR_NULL ("invalid hmax", GSL_EINVAL); |
137 | } |
138 | |
139 | return GSL_SUCCESS; |
140 | } |
141 | |
142 | int |
143 | gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d, |
144 | const unsigned long int nmax) |
145 | { |
146 | /* Sets maximum number of allowed steps (nmax) for driver */ |
147 | |
148 | d->nmax = nmax; |
149 | |
150 | return GSL_SUCCESS; |
151 | } |
152 | |
153 | gsl_odeiv2_driver * |
154 | gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * sys, |
155 | const gsl_odeiv2_step_type * T, |
156 | const double hstart, |
157 | const double epsabs, const double epsrel) |
158 | { |
159 | /* Initializes an ODE driver system with control object of type y_new. */ |
160 | |
161 | gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T); |
162 | |
163 | if (state == NULL) |
164 | { |
165 | GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM); |
166 | } |
167 | |
168 | if (epsabs >= 0.0 && epsrel >= 0.0) |
169 | { |
170 | state->c = gsl_odeiv2_control_y_new (epsabs, epsrel); |
171 | |
172 | if (state->c == NULL) |
173 | { |
174 | gsl_odeiv2_driver_free (state); |
175 | GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM); |
176 | } |
177 | } |
178 | else |
179 | { |
180 | gsl_odeiv2_driver_free (state); |
181 | GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL); |
182 | } |
183 | |
184 | /* Distribute pointer to driver object */ |
185 | |
186 | gsl_odeiv2_step_set_driver (state->s, state); |
187 | gsl_odeiv2_evolve_set_driver (state->e, state); |
188 | gsl_odeiv2_control_set_driver (state->c, state); |
189 | |
190 | return state; |
191 | } |
192 | |
193 | gsl_odeiv2_driver * |
194 | gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * sys, |
195 | const gsl_odeiv2_step_type * T, |
196 | const double hstart, |
197 | const double epsabs, const double epsrel) |
198 | { |
199 | /* Initializes an ODE driver system with control object of type yp_new. */ |
200 | |
201 | gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T); |
202 | |
203 | if (state == NULL) |
204 | { |
205 | GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM); |
206 | } |
207 | |
208 | if (epsabs >= 0.0 && epsrel >= 0.0) |
209 | { |
210 | state->c = gsl_odeiv2_control_yp_new (epsabs, epsrel); |
211 | |
212 | if (state->c == NULL) |
213 | { |
214 | gsl_odeiv2_driver_free (state); |
215 | GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM); |
216 | } |
217 | } |
218 | else |
219 | { |
220 | gsl_odeiv2_driver_free (state); |
221 | GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL); |
222 | } |
223 | |
224 | /* Distribute pointer to driver object */ |
225 | |
226 | gsl_odeiv2_step_set_driver (state->s, state); |
227 | gsl_odeiv2_evolve_set_driver (state->e, state); |
228 | gsl_odeiv2_control_set_driver (state->c, state); |
229 | |
230 | return state; |
231 | } |
232 | |
233 | gsl_odeiv2_driver * |
234 | gsl_odeiv2_driver_alloc_standard_new (const gsl_odeiv2_system * sys, |
235 | const gsl_odeiv2_step_type * T, |
236 | const double hstart, |
237 | const double epsabs, |
238 | const double epsrel, const double a_y, |
239 | const double a_dydt) |
240 | { |
241 | /* Initializes an ODE driver system with control object of type |
242 | standard_new. |
243 | */ |
244 | |
245 | gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T); |
246 | |
247 | if (state == NULL) |
248 | { |
249 | GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM); |
250 | } |
251 | |
252 | if (epsabs >= 0.0 && epsrel >= 0.0) |
253 | { |
254 | state->c = |
255 | gsl_odeiv2_control_standard_new (epsabs, epsrel, a_y, a_dydt); |
256 | |
257 | if (state->c == NULL) |
258 | { |
259 | gsl_odeiv2_driver_free (state); |
260 | GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM); |
261 | } |
262 | } |
263 | else |
264 | { |
265 | gsl_odeiv2_driver_free (state); |
266 | GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL); |
267 | } |
268 | |
269 | /* Distribute pointer to driver object */ |
270 | |
271 | gsl_odeiv2_step_set_driver (state->s, state); |
272 | gsl_odeiv2_evolve_set_driver (state->e, state); |
273 | gsl_odeiv2_control_set_driver (state->c, state); |
274 | |
275 | return state; |
276 | } |
277 | |
278 | gsl_odeiv2_driver * |
279 | gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system * sys, |
280 | const gsl_odeiv2_step_type * T, |
281 | const double hstart, |
282 | const double epsabs, const double epsrel, |
283 | const double a_y, const double a_dydt, |
284 | const double scale_abs[]) |
285 | { |
286 | /* Initializes an ODE driver system with control object of type |
287 | scaled_new. |
288 | */ |
289 | |
290 | gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T); |
291 | |
292 | if (state == NULL) |
293 | { |
294 | GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM); |
295 | } |
296 | |
297 | if (epsabs >= 0.0 && epsrel >= 0.0) |
298 | { |
299 | state->c = gsl_odeiv2_control_scaled_new (epsabs, epsrel, a_y, a_dydt, |
300 | scale_abs, sys->dimension); |
301 | |
302 | if (state->c == NULL) |
303 | { |
304 | gsl_odeiv2_driver_free (state); |
305 | GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM); |
306 | } |
307 | } |
308 | else |
309 | { |
310 | gsl_odeiv2_driver_free (state); |
311 | GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL); |
312 | } |
313 | |
314 | /* Distribute pointer to driver object */ |
315 | |
316 | gsl_odeiv2_step_set_driver (state->s, state); |
317 | gsl_odeiv2_evolve_set_driver (state->e, state); |
318 | gsl_odeiv2_control_set_driver (state->c, state); |
319 | |
320 | return state; |
321 | } |
322 | |
323 | int |
324 | gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double *t, |
325 | const double t1, double y[]) |
326 | { |
327 | /* Main driver function that evolves the system from t to t1. In |
328 | beginning vector y contains the values of dependent variables at |
329 | t. This function returns values at t=t1 in y. In case of |
330 | unrecoverable error, y and t contains the values after the last |
331 | successful step. |
332 | */ |
333 | |
334 | int sign = 0; |
335 | d->n = 0; |
336 | |
337 | /* Determine integration direction sign */ |
338 | |
339 | if (d->h > 0.0) |
340 | { |
341 | sign = 1; |
342 | } |
343 | else |
344 | { |
345 | sign = -1; |
346 | } |
347 | |
348 | /* Check that t, t1 and step direction are sensible */ |
349 | |
350 | if (sign * (t1 - *t) < 0.0) |
351 | { |
352 | GSL_ERROR_NULL |
353 | ("integration limits and/or step direction not consistent", |
354 | GSL_EINVAL); |
355 | } |
356 | |
357 | /* Evolution loop */ |
358 | |
359 | while (sign * (t1 - *t) > 0.0) |
360 | { |
361 | int s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, d->sys, |
362 | t, t1, &(d->h), y); |
363 | |
364 | if (s != GSL_SUCCESS) |
365 | { |
366 | return s; |
367 | } |
368 | |
369 | /* Check for maximum allowed steps */ |
370 | |
371 | if ((d->nmax > 0) && (d->n > d->nmax)) |
372 | { |
373 | return GSL_EMAXITER; |
374 | } |
375 | |
376 | /* Set step size if maximum size is exceeded */ |
377 | |
378 | if (fabs (d->h) > d->hmax) |
379 | { |
380 | d->h = sign * d->hmax; |
381 | } |
382 | |
383 | /* Check for too small step size */ |
384 | |
385 | if (fabs (d->h) < d->hmin) |
386 | { |
387 | return GSL_ENOPROG; |
388 | } |
389 | |
390 | d->n++; |
391 | } |
392 | |
393 | return GSL_SUCCESS; |
394 | } |
395 | |
396 | int |
397 | gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double *t, |
398 | const double h, const unsigned long int n, |
399 | double y[]) |
400 | { |
401 | /* Alternative driver function that evolves the system from t using |
402 | * n steps of size h. In the beginning vector y contains the values |
403 | * of dependent variables at t. This function returns values at t = |
404 | * t + n * h in y. In case of an unrecoverable error, y and t |
405 | * contains the values after the last successful step. |
406 | */ |
407 | |
408 | unsigned long int i; |
409 | d->n = 0; |
410 | |
411 | /* Evolution loop */ |
412 | |
413 | for (i = 0; i < n; i++) |
414 | { |
415 | int s = gsl_odeiv2_evolve_apply_fixed_step (d->e, d->c, d->s, d->sys, |
416 | t, h, y); |
417 | |
418 | if (s != GSL_SUCCESS) |
419 | { |
420 | return s; |
421 | } |
422 | |
423 | d->n++; |
424 | } |
425 | |
426 | return GSL_SUCCESS; |
427 | } |
428 | |
429 | int |
430 | gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d) |
431 | { |
432 | /* Reset the driver. Resets evolve and step objects. */ |
433 | |
434 | { |
435 | int s = gsl_odeiv2_evolve_reset (d->e); |
436 | |
437 | if (s != GSL_SUCCESS) |
438 | { |
439 | return s; |
440 | } |
441 | } |
442 | |
443 | { |
444 | int s = gsl_odeiv2_step_reset (d->s); |
445 | |
446 | if (s != GSL_SUCCESS) |
447 | { |
448 | return s; |
449 | } |
450 | } |
451 | |
452 | return GSL_SUCCESS; |
453 | } |
454 | |
455 | int |
456 | gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart) |
457 | { |
458 | /* Resets current driver and sets initial step size to hstart */ |
459 | |
460 | gsl_odeiv2_driver_reset (d); |
461 | |
462 | if ((d->hmin > fabs (hstart)) || (fabs (hstart) > d->hmax)) |
463 | { |
464 | GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL); |
465 | } |
466 | |
467 | if (hstart > 0.0 || hstart < 0.0) |
468 | { |
469 | d->h = hstart; |
470 | } |
471 | else |
472 | { |
473 | GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL); |
474 | } |
475 | |
476 | return GSL_SUCCESS; |
477 | } |
478 | |
479 | void |
480 | gsl_odeiv2_driver_free (gsl_odeiv2_driver * state) |
481 | { |
482 | if (state->c != NULL) |
483 | { |
484 | gsl_odeiv2_control_free (state->c); |
485 | } |
486 | |
487 | gsl_odeiv2_evolve_free (state->e); |
488 | gsl_odeiv2_step_free (state->s); |
489 | free (state); |
490 | } |
491 | |