Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit 274dc29

Browse files
hakonhaglandhakonhagland
committedOct 13, 2019
Fixes issue leto#156 for gsl_odeiv_system.
Implements some typemaps for Math::GSL::ODEIV to fix issue leto#156 for gsl_odeiv_system and gsl_odeiv_evolve_apply().
1 parent 88e195a commit 274dc29

File tree

3 files changed

+794
-6
lines changed

3 files changed

+794
-6
lines changed
 

‎pod/ODEIV.pod

+56-2
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,6 @@
11
%perlcode %{
2+
package Math::GSL::ODEIV;
3+
24
@EXPORT_OK = qw/
35
gsl_odeiv_step_alloc
46
gsl_odeiv_step_reset
@@ -47,7 +49,7 @@ Math::GSL::ODEIV - functions for solving ordinary differential equation (ODE) in
4749

4850
=head1 SYNOPSIS
4951

50-
use Math::GSL::ODEIV qw /:all/;
52+
use Math::GSL::ODEIV qw /:all/;
5153

5254
=head1 DESCRIPTION
5355

@@ -87,7 +89,7 @@ Here is a list of all the functions in this module :
8789

8890
=item * C<gsl_odeiv_evolve_alloc($dim)> - This function returns a pointer to a newly allocated instance of an evolution function for a system of $dim dimensions.
8991

90-
=item * C<gsl_odeiv_evolve_apply >
92+
=item * C<gsl_odeiv_evolve_apply($e, $c, $step, $dydt, \$t, $t1, \$h, $y)> - This function advances the system ($e, $dydt) from time $t and position $y using the stepping function $step. The new time and position are stored in $t and $y on output. The initial step-size is taken as $h, but this will be modified using the control function $c to achieve the appropriate error bound if necessary. The routine may make several calls to step in order to determine the optimum step-size. If the step-size has been changed the value of $h will be modified on output. The maximum time $t1 is guaranteed not to be exceeded by the time-step. On the final time-step the value of $t will be set to $t1 exactly.
9193

9294
=item * C<gsl_odeiv_evolve_reset($e)> - This function resets the evolution function $e. It should be used whenever the next use of $e will not be a continuation of a previous step.
9395

@@ -138,6 +140,58 @@ This module also includes the following constants :
138140
For more informations on the functions, we refer you to the GSL offcial
139141
documentation: L<http://www.gnu.org/software/gsl/manual/html_node/>
140142

143+
=head1 EXAMPLE
144+
145+
The example is taken from L<https://www.math.utah.edu/software/gsl/gsl-ref_367.html>.
146+
147+
use strict;
148+
use warnings;
149+
use Math::GSL::Errno qw($GSL_SUCCESS);
150+
use Math::GSL::ODEIV qw/ :all /;
151+
use Math::GSL::Matrix qw/:all/;
152+
use Math::GSL::IEEEUtils qw/ :all /;
153+
154+
sub func {
155+
my ($t, $y, $dydt, $params) = @_;
156+
my $mu = $params->{mu};
157+
$dydt->[0] = $y->[1];
158+
$dydt->[1] = -$y->[0] - $mu*$y->[1]*(($y->[0])**2 - 1);
159+
return $GSL_SUCCESS;
160+
}
161+
162+
sub jac {
163+
my ($t, $y, $dfdy, $dfdt, $params) = @_;
164+
165+
my $mu = $params->{mu};
166+
my $m = gsl_matrix_view_array($dfdy, 2, 2);
167+
gsl_matrix_set( $m, 0, 0, 0.0 );
168+
gsl_matrix_set( $m, 0, 1, 1.0 );
169+
gsl_matrix_set( $m, 1, 0, (-2.0 * $mu * $y->[0] * $y->[1]) - 1.0 );
170+
gsl_matrix_set( $m, 1, 1, -$mu * (($y->[0])**2 - 1.0) );
171+
$dfdt->[0] = 0.0;
172+
$dfdt->[1] = 0.0;
173+
return $GSL_SUCCESS;
174+
}
175+
176+
my $T = $gsl_odeiv_step_rk8pd;
177+
my $s = gsl_odeiv_step_alloc($T, 2);
178+
my $c = gsl_odeiv_control_y_new(1e-6, 0.0);
179+
my $e = gsl_odeiv_evolve_alloc(2);
180+
my $params = { mu => 10 };
181+
my $sys = Math::GSL::ODEIV::gsl_odeiv_system->new(\&func, \&jac, 2, $params );
182+
my $t = 0.0;
183+
my $t1 = 100.0;
184+
my $h = 1e-6;
185+
my $y = [ 1.0, 0.0 ];
186+
gsl_ieee_env_setup;
187+
while ($t < $t1) {
188+
my $status = gsl_odeiv_evolve_apply ($e, $c, $s, $sys, \$t, $t1, \$h, $y);
189+
last if $status != $GSL_SUCCESS;
190+
printf "%.5e %.5e %.5e\n", $t, $y->[0], $y->[1];
191+
}
192+
gsl_odeiv_evolve_free($e);
193+
gsl_odeiv_control_free($c);
194+
gsl_odeiv_step_free($s);
141195

142196

143197

‎swig/ODEIV.i

+508-4
Original file line numberDiff line numberDiff line change
@@ -3,15 +3,519 @@
33
%include "typemaps.i"
44
%include "renames.i"
55

6-
%newobject gsl_odeiv_step_alloc;
6+
//%newobject gsl_odeiv_step_alloc;
77
%newobject gsl_odeiv_driver_alloc_y_new;
8-
%newobject gsl_odeiv_evolve_alloc;
8+
//%newobject gsl_odeiv_evolve_alloc;
99

1010
%{
11-
#include "gsl/gsl_types.h"
12-
#include "gsl/gsl_odeiv.h"
11+
//1
12+
#include "gsl/gsl_types.h"
13+
#include "gsl/gsl_odeiv.h"
14+
#include <stdarg.h>
15+
#define SWIG_MATH_GSL_ODEIV_PACKAGE_NAME "Math::GSL::ODEIV"
16+
#define SWIG_MATH_GSL_ODEIV_GUTS \
17+
SWIG_MATH_GSL_ODEIV_PACKAGE_NAME "::_guts"
18+
19+
/* NOTE: We do not use C static variables to store the values since
20+
static variable are not thread-safe, and SWIG does not currently support the
21+
MY_CXT framework, see perldoc perlxs for more information.
22+
Perl package variables are on the other hand thread safe.
23+
*/
24+
void swig_math_gsl_odeiv_set_gut_pv(
25+
const char *varname, const char *value
26+
)
27+
{
28+
SV *pvname = newSVpvf( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);
29+
/*char *pvname = form( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);*/
30+
SV* sv = get_sv( SvPV_nolen(pvname), GV_ADD );
31+
SvREFCNT_dec(pvname);
32+
sv_setpv( sv, value );
33+
}
34+
35+
const char *swig_math_gsl_odeiv_get_gut_pv(const char *varname)
36+
{
37+
SV *pvname = newSVpvf( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);
38+
/*char *pvname = form( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);*/
39+
SV* sv = get_sv( SvPV_nolen(pvname), GV_ADD );
40+
SvREFCNT_dec(pvname);
41+
return SvPV_nolen(sv);
42+
}
43+
44+
void swig_math_gsl_odeiv_set_callback_error_param( const char *func )
45+
{
46+
swig_math_gsl_odeiv_set_gut_pv( "cbname", func );
47+
}
48+
49+
void swig_math_gsl_odeiv_set_error_param(
50+
const char *symname, const char *param
51+
)
52+
{
53+
swig_math_gsl_odeiv_set_gut_pv( "symname", symname );
54+
swig_math_gsl_odeiv_set_gut_pv( "param", param );
55+
}
56+
57+
void swig_math_gsl_odeiv_callback_error(
58+
const char *msg, ...
59+
)
60+
{
61+
const char *cbname = swig_math_gsl_odeiv_get_gut_pv("cbname");
62+
va_list args;
63+
va_start(args, msg);
64+
/*char *msg2 = form(
65+
"Math::GSL::ODEIV: callback function : %s() : %s", cbname, msg
66+
);*/
67+
SV *msg2 = newSVpvf(
68+
"Math::GSL::ODEIV: callback function : %s() : %s", cbname, msg
69+
);
70+
vcroak( SvPV_nolen(msg2), &args );
71+
/* NOTE: these two lines will never be reached */
72+
SvREFCNT_dec(msg2);
73+
va_end(args);
74+
}
75+
76+
void swig_math_gsl_odeiv_input_param_error(
77+
const char *msg, ...
78+
)
79+
{
80+
const char *subname = swig_math_gsl_odeiv_get_gut_pv("symname");
81+
const char *param = swig_math_gsl_odeiv_get_gut_pv("param");
82+
va_list args;
83+
va_start(args, msg);
84+
SV *msg2 = newSVpvf(
85+
"Math::GSL::ODEIV:%s() : parameter $%s : %s", subname, param, msg
86+
);
87+
vcroak( SvPV_nolen(msg2), &args );
88+
/* NOTE: these two lines will never be reached */
89+
SvREFCNT_dec(msg2);
90+
va_end(args);
91+
}
92+
93+
void swig_math_gsl_odeiv_input_error(
94+
const char *msg, ...
95+
)
96+
{
97+
const char *subname = swig_math_gsl_odeiv_get_gut_pv("symname");
98+
va_list args;
99+
va_start(args, msg);
100+
SV *msg2 = newSVpvf( "Math::GSL::ODEIV:%s() : %s", subname, msg);
101+
vcroak( SvPV_nolen(msg2), &args );
102+
/* NOTE: these two lines will never be reached */
103+
SvREFCNT_dec(msg2);
104+
va_end(args);
105+
}
106+
107+
SV *swig_math_gsl_odeiv_get_hash_sv(HV *hash, const char *key)
108+
{
109+
SV *key_sv = newSVpv(key, strlen (key));
110+
SV *value;
111+
if (hv_exists_ent(hash, key_sv, 0)) {
112+
HE *he = hv_fetch_ent(hash, key_sv, 0, 0);
113+
value = HeVAL(he);
114+
}
115+
else {
116+
swig_math_gsl_odeiv_input_param_error(
117+
"The hash key '%s' doesn't exist", key
118+
);
119+
}
120+
return value;
121+
}
122+
123+
IV swig_math_gsl_odeiv_get_hash_iv(HV *hash, const char *key) {
124+
SV *sv = swig_math_gsl_odeiv_get_hash_sv(hash, key);
125+
if (SvROK(sv)) {
126+
swig_math_gsl_odeiv_input_param_error(
127+
"Hash value for key '%s' is not a scalar value", key
128+
);
129+
}
130+
if (!SvIOK(sv)) {
131+
swig_math_gsl_odeiv_input_param_error(
132+
"Hash value for key '%s' is not an integer", key
133+
);
134+
}
135+
return SvIV(sv);
136+
}
137+
138+
SV *swig_math_gsl_odeiv_get_hash_hashref(HV *hash, const char *key) {
139+
SV *sv = swig_math_gsl_odeiv_get_hash_sv(hash, key);
140+
if (!SvROK(sv)) {
141+
swig_math_gsl_odeiv_input_param_error(
142+
"Hash value for key '%s' is not a reference", key
143+
);
144+
}
145+
if (SvTYPE(SvRV(sv)) != SVt_PVHV) {
146+
swig_math_gsl_odeiv_input_param_error(
147+
"Hash value for key '%s' is not a hashref", key
148+
);
149+
}
150+
return sv;
151+
}
152+
153+
SV *swig_math_gsl_odeiv_get_hash_coderef(HV *hash, const char *key) {
154+
SV *sv = swig_math_gsl_odeiv_get_hash_sv(hash, key);
155+
if (!SvROK(sv)) {
156+
swig_math_gsl_odeiv_input_param_error(
157+
"Hash value for key '%s' is not a reference", key
158+
);
159+
}
160+
if (SvTYPE(SvRV(sv)) != SVt_PVCV) {
161+
swig_math_gsl_odeiv_input_param_error(
162+
"Hash value for key '%s' is not a coderef", key
163+
);
164+
}
165+
return sv;
166+
}
167+
168+
void swig_math_gsl_odeiv_store_hash_ptr( HV *hash, const char *key, void *ptr)
169+
{
170+
SV *sv = newSViv(PTR2IV(ptr));
171+
/* Let the hash take ownership of the sv */
172+
if( !hv_store(hash, key, strlen(key), sv, 0) ) {
173+
SvREFCNT_dec(sv);
174+
swig_math_gsl_odeiv_input_param_error(
175+
"Cannot store internal information in the hash"
176+
);
177+
}
178+
}
179+
180+
void *swig_math_gsl_odeiv_get_hash_ptr(HV *hash, const char *key) {
181+
IV ptr = swig_math_gsl_odeiv_get_hash_iv(hash, key);
182+
return (void *) INT2PTR(SV*, ptr);
183+
}
184+
185+
void swig_math_gsl_odeiv_store_double_in_av( AV *array, SSize_t index, double val)
186+
{
187+
SV *sval = newSVnv(val);
188+
if( !av_store(array, index, sval) ) {
189+
SvREFCNT_dec(sval);
190+
swig_math_gsl_odeiv_callback_error(
191+
"Cannot store internal information in array"
192+
);
193+
}
194+
}
195+
196+
void swig_math_gsl_odeiv_copy_av_to_carray(AV *array, double *y, size_t dim)
197+
{
198+
SSize_t array_len = av_top_index(array) + 1;
199+
if (array_len != dim ) {
200+
swig_math_gsl_odeiv_callback_error(
201+
"Callback returned array of wrong dimension"
202+
);
203+
}
204+
for (int i = 0; i < dim; i++) {
205+
SV **sv_ptr = av_fetch( array, i, 0 );
206+
if (!sv_ptr) {
207+
swig_math_gsl_odeiv_callback_error(
208+
"Cannot extract values from returned array"
209+
);
210+
}
211+
SV *sv = *sv_ptr;
212+
if (SvROK(sv)) {
213+
swig_math_gsl_odeiv_callback_error(
214+
"Returned array value is not a scalar"
215+
);
216+
}
217+
/* TODO: Not sure if this check is necessary */
218+
if (SvTYPE(sv) >= SVt_PVAV) {
219+
swig_math_gsl_odeiv_callback_error(
220+
"Returned array value is not of scalar type"
221+
);
222+
}
223+
double val = (double ) SvNV(sv);
224+
y[i] = val;
225+
}
226+
}
227+
228+
void swig_math_gsl_odeiv_copy_doubles_to_av(AV *array, const double *y, size_t dim)
229+
{
230+
for (int i = 0; i < dim; i++) {
231+
swig_math_gsl_odeiv_store_double_in_av(array, i, y[i]);
232+
}
233+
}
234+
235+
typedef struct {
236+
SV *func;
237+
SV *jac;
238+
SV *params;
239+
size_t dim;
240+
gsl_odeiv_system *sys;
241+
} swig_math_gsl_odeiv_system;
242+
243+
int swig_math_gsl_odeiv_call_perl_jac (
244+
SV *callback, double t, const double y[], double *dfdy, double *dfdt,
245+
swig_math_gsl_odeiv_system *params
246+
)
247+
{
248+
AV *ay = (AV *)sv_2mortal((SV *)newAV());
249+
AV *a_dfdy = (AV *)sv_2mortal((SV *)newAV());
250+
AV *a_dfdt = (AV *)sv_2mortal((SV *)newAV());
251+
dSP; /* declares a local copy of stack pointer */
252+
ENTER;
253+
SAVETMPS;
254+
PUSHMARK(SP);
255+
EXTEND(SP, 5);
256+
mPUSHs(newSVnv(t));
257+
swig_math_gsl_odeiv_copy_doubles_to_av(ay, y, params->dim);
258+
mPUSHs((SV *)newRV_inc((SV *) ay));
259+
mPUSHs((SV *)newRV_inc((SV *) a_dfdy));
260+
mPUSHs((SV *)newRV_inc((SV *) a_dfdt));
261+
XPUSHs(params->params);
262+
PUTBACK;
263+
int count = call_sv(callback, G_SCALAR); /* call the Perl callback */
264+
SPAGAIN;
265+
if (count != 1) {
266+
swig_math_gsl_odeiv_callback_error(
267+
"Bad return value from callback: expected 1 value, got %d", count
268+
);
269+
}
270+
IV result = POPi; /* TODO: check ST(0) instead for valid value */
271+
swig_math_gsl_odeiv_copy_av_to_carray(a_dfdy, dfdy, (params->dim)*(params->dim));
272+
swig_math_gsl_odeiv_copy_av_to_carray(a_dfdt, dfdt, params->dim);
273+
PUTBACK;
274+
FREETMPS;
275+
LEAVE;
276+
return (int) result;
277+
}
278+
279+
int swig_math_gsl_odeiv_call_perl_func (
280+
SV *callback, double t, const double y[], double dydt[],
281+
swig_math_gsl_odeiv_system *params)
282+
{
283+
AV *ay = (AV *)sv_2mortal((SV *)newAV());
284+
AV *aj = (AV *)sv_2mortal((SV *)newAV());
285+
dSP; /* declares a local copy of stack pointer */
286+
ENTER;
287+
SAVETMPS;
288+
PUSHMARK(SP);
289+
EXTEND(SP, 4);
290+
mPUSHs(newSVnv(t));
291+
swig_math_gsl_odeiv_copy_doubles_to_av(ay, y, params->dim);
292+
mPUSHs((SV *)newRV_inc((SV *) ay));
293+
mPUSHs((SV *)newRV_inc((SV *) aj));
294+
XPUSHs(params->params);
295+
PUTBACK;
296+
int count = call_sv(callback, G_SCALAR); /* call the Perl callback */
297+
SPAGAIN;
298+
/* This should happen for G_SCALAR, see perldoc perlcall.
299+
* Even if the callback does not return anything, count will still be 1
300+
* since we are not the G_DISCARD flag
301+
*/
302+
if (count != 1) {
303+
swig_math_gsl_odeiv_callback_error(
304+
"Bad return value from callback: expected 1 value, got %d", count
305+
);
306+
}
307+
IV result = POPi; /* TODO: check ST(0) instead for valid value */
308+
swig_math_gsl_odeiv_copy_av_to_carray(aj, dydt, params->dim);
309+
PUTBACK;
310+
FREETMPS;
311+
LEAVE;
312+
return result;
313+
}
314+
315+
int swig_math_gsl_odeiv_callback_function(
316+
double t, const double y[], double dydt[], void *params)
317+
{
318+
swig_math_gsl_odeiv_set_callback_error_param( "func" );
319+
swig_math_gsl_odeiv_system *sparam = (swig_math_gsl_odeiv_system *) params;
320+
return swig_math_gsl_odeiv_call_perl_func(sparam->func, t, y, dydt, sparam);
321+
}
322+
323+
int swig_math_gsl_odeiv_callback_jac(
324+
double t, const double y[], double *dfdy, double *dfdt, void *params
325+
)
326+
{
327+
swig_math_gsl_odeiv_set_callback_error_param( "jac" );
328+
swig_math_gsl_odeiv_system *sparam = (swig_math_gsl_odeiv_system *) params;
329+
return swig_math_gsl_odeiv_call_perl_jac(sparam->jac, t, y, dfdy, dfdt, sparam);
330+
}
331+
332+
void swig_math_gsl_odeiv_fill_system_struct( HV *hash, gsl_odeiv_system *sys )
333+
{
334+
swig_math_gsl_odeiv_system *ssys;
335+
Newx(ssys, 1, swig_math_gsl_odeiv_system);
336+
ssys->func = swig_math_gsl_odeiv_get_hash_coderef( hash, "func" );
337+
ssys->jac = swig_math_gsl_odeiv_get_hash_coderef( hash, "jac" );
338+
ssys->params = swig_math_gsl_odeiv_get_hash_hashref( hash, "params" );
339+
ssys->dim = (size_t) swig_math_gsl_odeiv_get_hash_iv( hash, "dim" );
340+
sys->function = swig_math_gsl_odeiv_callback_function;
341+
sys->jacobian = swig_math_gsl_odeiv_callback_jac;
342+
sys->dimension = (size_t) swig_math_gsl_odeiv_get_hash_iv( hash, "dim" );
343+
sys->params = (void *) ssys;
344+
}
13345
%}
14346

347+
%typemap(in) const gsl_odeiv_system * {
348+
SV *sv;
349+
HV *hash;
350+
351+
swig_math_gsl_odeiv_set_error_param( "$symname", "$1_name" );
352+
353+
if (!sv_isobject($input)) {
354+
swig_math_gsl_odeiv_input_error(
355+
"Input parameter $$1_name is not a blessed reference!"
356+
);
357+
}
358+
if (!sv_isa($input, "Math::GSL::ODEIV::gsl_odeiv_system")) {
359+
swig_math_gsl_odeiv_input_error(
360+
"Input parameter $$1_name is not an object of type "
361+
"Math::GSL::ODEIV::gsl_odeiv_system!"
362+
);
363+
}
364+
sv = SvRV($input);
365+
if ((SvTYPE(sv) != SVt_PVHV)) {
366+
swig_math_gsl_odeiv_input_error(
367+
"Input parameter $$1_name is not a blessed hash reference!"
368+
);
369+
}
370+
gsl_odeiv_system *sys;
371+
Newx(sys, 1, gsl_odeiv_system);
372+
swig_math_gsl_odeiv_fill_system_struct( (HV *)sv, sys);
373+
$1 = sys;
374+
}
375+
376+
%typemap(freearg) const gsl_odeiv_system * {
377+
swig_math_gsl_odeiv_system *ssys = (swig_math_gsl_odeiv_system *)$1->params;
378+
Safefree(ssys);
379+
Safefree($1);
380+
}
381+
382+
%ignore gsl_odeiv_evolve_apply;
15383
%include "gsl/gsl_types.h"
16384
%include "gsl/gsl_odeiv.h"
385+
// unignore gsl_odeiv_evolve_apply()
386+
%rename("%s") gsl_odeiv_evolve_apply;
387+
/* We want handle the last parameter to gsl_odeiv_evolve_apply(...), see
388+
* include file gsl/gsl_odeiv.h for a definition. This parameter is of
389+
* type 'double []' and acts as an input-output array.
390+
*
391+
* NOTE: gsl_typemaps.i has typemaps for a float [] input-output array,
392+
* but note that that typemap also returns the array elements on the perl stack
393+
* (in addition to modifying the original array).
394+
* However, here we do not want to return the result on the stack, we only
395+
* want to modify the original array.
396+
*
397+
* TODO: These typemaps might warrant to be moved to gsl_typemaps.i at a later time,
398+
* where they could be reused by other interface files, however currently they are
399+
* regarded as specific to only gsl_odeiv_evolve_apply().
400+
*/
401+
402+
%typemap(in) double y[] {
403+
struct perl_array * p_array = 0; /* see gsl_typemaps.i for definition */
404+
I32 len;
405+
AV *array;
406+
int i;
407+
SV **tv;
408+
swig_math_gsl_odeiv_set_error_param( "$symname", "$1_name" );
409+
if (!SvROK($input)) {
410+
swig_math_gsl_odeiv_input_error(
411+
"Input parameter $$1_name is not a reference!"
412+
);
413+
}
414+
if (SvTYPE(SvRV($input)) != SVt_PVAV) {
415+
swig_math_gsl_odeiv_input_error(
416+
"Input parameter $$1_name is not an array reference!"
417+
);
418+
}
419+
array = (AV*)SvRV($input);
420+
len = av_len(array);
421+
p_array = (struct perl_array *)
422+
malloc((len+1)*sizeof(double)+sizeof(struct perl_array));
423+
p_array->len=len;
424+
p_array->array=array;
425+
$1 = (double *)&p_array[1];
426+
for (i = 0; i <= len; i++) {
427+
tv = av_fetch(array, i, 0);
428+
$1[i] = (double) SvNV(*tv);
429+
}
430+
}
431+
432+
%typemap(argout) double y[] {
433+
struct perl_array * p_array = 0; /* see gsl_typemaps.i for definition */
434+
int i;
435+
SV **tv;
436+
p_array=(struct perl_array *)(((char*)$1)-sizeof(struct perl_array));
437+
for (i = 0; i <= p_array->len; i++) {
438+
double val= $1[i];
439+
tv = av_fetch(p_array->array, i, 0);
440+
sv_setnv(*tv, val);
441+
}
442+
}
443+
444+
%typemap(freearg) double y[] {
445+
if ($1) free(((char*)$1)-sizeof(struct perl_array));
446+
}
447+
448+
%{
449+
typedef struct {
450+
double h;
451+
SV *sv;
452+
} swig_perl_double_in_out;
453+
454+
%}
455+
456+
%typemap(in) double *h {
457+
swig_math_gsl_odeiv_set_error_param( "$symname", "$1_name" );
458+
if (!SvROK($input)) {
459+
swig_math_gsl_odeiv_input_error(
460+
"Input parameter $$1_name is not a reference!"
461+
);
462+
}
463+
if (SvTYPE(SvRV($input)) >= SVt_PVAV) {
464+
swig_math_gsl_odeiv_input_error(
465+
"Input parameter $$1_name is not a scalar reference!"
466+
);
467+
}
468+
SV *sv = SvRV($input);
469+
swig_perl_double_in_out *h_wrap;
470+
Newx(h_wrap, 1, swig_perl_double_in_out);
471+
h_wrap->sv = sv;
472+
h_wrap->h = (double) SvNV(sv);
473+
$1 = (double *) &h_wrap->h;
474+
}
475+
476+
%typemap(argout) double *h {
477+
swig_perl_double_in_out *h_wrap = (swig_perl_double_in_out *) $1;
478+
SV *sv = h_wrap->sv;
479+
sv_setnv(sv, h_wrap->h);
480+
}
481+
482+
%typemap(freearg) double *h {
483+
Safefree( $1 );
484+
}
485+
486+
487+
%typemap(in) double *t = double *h;
488+
%typemap(argout) double *t = double *h;
489+
%typemap(freearg) double *t = double *h;
490+
491+
// define our own name for the last parameter to gsl_odeiv_evolve_apply()
492+
int gsl_odeiv_evolve_apply(gsl_odeiv_evolve *e, gsl_odeiv_control *con, gsl_odeiv_step *step, const gsl_odeiv_system *dydt, double *t, double t1, double *h, double y[]);
493+
494+
%perlcode %{
495+
package Math::GSL::ODEIV::gsl_odeiv_system;
496+
497+
sub new {
498+
my ($class, $func, $jac, $dim, $params ) = @_;
499+
500+
my $check_ref = sub {
501+
if ( (ref $_[0]) ne $_[1] ) {
502+
die sprintf 'Usage: %s:new( $func, $jac, $dim, $params ). '
503+
. 'Argument %s is not %s reference',
504+
__PACKAGE__, $_[2], $_[3];
505+
}
506+
};
507+
my $check_subref = sub {
508+
$check_ref->($_[0], "CODE", $_[1], "code");
509+
};
510+
my $check_hashref = sub {
511+
$check_ref->($_[0], "HASH", $_[1], "hash");
512+
};
513+
$check_subref->($func, '$func');
514+
$check_subref->($jac, '$jac');
515+
$check_hashref->($params, '$params');
516+
return bless { func => $func, jac => $jac, dim => $dim, params => $params },
517+
$class;
518+
}
519+
%}
520+
17521
%include "../pod/ODEIV.pod"

‎t/ODEIV_evolve.t

+230
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,230 @@
1+
package Math::GSL::ODEIV::Test;
2+
use strict;
3+
use base q{Test::Class};
4+
use Test::Exception;
5+
use Test::More;
6+
use Math::GSL qw/:all/;
7+
use Math::GSL::IEEEUtils qw/ :all /;
8+
use Math::GSL::ODEIV qw/:all/;
9+
use Math::GSL::Errno qw/:all/;
10+
use Math::GSL::Test qw/:all/;
11+
use Data::Dumper;
12+
13+
BEGIN { gsl_set_error_handler_off(); }
14+
15+
sub make_fixture : Test(setup) {
16+
my $self = shift;
17+
my %config;
18+
my $T = $config{T} = $gsl_odeiv_step_rk8pd;
19+
$config{s} = gsl_odeiv_step_alloc($T, 2);
20+
$config{c} = gsl_odeiv_control_y_new(1e-6, 0.0);
21+
$config{e} = gsl_odeiv_evolve_alloc(2);
22+
my $params = $config{params} = { mu => 10 };
23+
$config{sys} = Math::GSL::ODEIV::gsl_odeiv_system->new(\&func, \&jac, 2, $params );
24+
gsl_ieee_env_setup;
25+
$self->{config} = \%config;
26+
}
27+
28+
sub teardown : Test(teardown) {
29+
my $self = shift;
30+
31+
my $c = $self->{config};
32+
gsl_odeiv_evolve_free($c->{e});
33+
gsl_odeiv_control_free($c->{c});
34+
gsl_odeiv_step_free($c->{s});
35+
}
36+
37+
sub func {
38+
my ($t, $y, $dydt, $params) = @_;
39+
my $mu = $params->{mu};
40+
41+
$dydt->[0] = $y->[1];
42+
$dydt->[1] = -$y->[0] - $mu*$y->[1]*(($y->[0])**2 - 1);
43+
return $GSL_SUCCESS;
44+
}
45+
46+
sub jac {
47+
my ($t, $y, $dfdy, $dfdt, $params) = @_;
48+
49+
my $mu = $params->{mu};
50+
my $m = gsl_matrix_view_array($dfdy, 2, 2);
51+
gsl_matrix_set( $m, 0, 0, 0.0 );
52+
gsl_matrix_set( $m, 0, 1, 1.0 );
53+
gsl_matrix_set( $m, 1, 0, (-2.0 * $mu * $y->[0] * $y->[1]) - 1.0 );
54+
gsl_matrix_set( $m, 1, 1, -$mu * (($y->[0])**2 - 1.0) );
55+
$dfdt->[0] = 0.0;
56+
$dfdt->[1] = 0.0;
57+
return $GSL_SUCCESS;
58+
}
59+
60+
sub GSL_ODEIV_BAD_SYS_INPUT : Tests {
61+
my $self = shift;
62+
63+
my $c = $self->{config};
64+
my $t = 0.0;
65+
my $t1 = 100.0;
66+
my $h = 1e-6;
67+
my $y = [ 1.0, 0.0 ];
68+
my $foo = 2;
69+
throws_ok {
70+
gsl_odeiv_evolve_apply (
71+
$c->{e}, $c->{c}, $c->{s}, $foo, \$t, $t1, \$h, $y
72+
);
73+
} qr/Input parameter \$dydt is not a blessed reference/, 'wrong type 1';
74+
throws_ok {
75+
gsl_odeiv_evolve_apply (
76+
$c->{e}, $c->{c}, $c->{s}, $self, \$t, $t1, \$h, $y
77+
);
78+
} qr/not an object of type Math::GSL::ODEIV::gsl_odeiv_system/, 'wrong type 2';
79+
my $bar = bless [], 'Math::GSL::ODEIV::gsl_odeiv_system';
80+
throws_ok {
81+
gsl_odeiv_evolve_apply (
82+
$c->{e}, $c->{c}, $c->{s}, $bar, \$t, $t1, \$h, $y
83+
);
84+
} qr/not a blessed hash reference/, 'wrong type 3';
85+
}
86+
87+
sub GSL_ODEIV_BAD_SYS_INPUT2 : Tests {
88+
my $self = shift;
89+
90+
my $c = $self->{config};
91+
my $t = 0.0;
92+
my $t1 = 100.0;
93+
my $h = 1e-6;
94+
my $y = [ 1.0, 0.0 ];
95+
my $sys = $c->{sys};
96+
$sys->{func} = 2;
97+
throws_ok {
98+
gsl_odeiv_evolve_apply (
99+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
100+
);
101+
} qr/not a reference/, 'wrong func type 1';
102+
$sys->{func} = [];
103+
throws_ok {
104+
gsl_odeiv_evolve_apply (
105+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
106+
);
107+
} qr/not a coderef/, 'wrong func type 2';
108+
$sys->{func} = \&func;
109+
$sys->{params} = 5;
110+
throws_ok {
111+
gsl_odeiv_evolve_apply (
112+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
113+
);
114+
} qr/not a reference/, 'wrong param type 1';
115+
$sys->{params} = [];
116+
throws_ok {
117+
gsl_odeiv_evolve_apply (
118+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
119+
);
120+
} qr/not a hashref/, 'wrong param type 2';
121+
$sys->{params} = { mu => 10 };
122+
$sys->{dim} = [];
123+
throws_ok {
124+
gsl_odeiv_evolve_apply (
125+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
126+
);
127+
} qr/not a scalar value/, 'wrong dim type 1';
128+
$sys->{dim} = 2.1;
129+
throws_ok {
130+
gsl_odeiv_evolve_apply (
131+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
132+
);
133+
} qr/not an integer/, 'wrong dim type 2';
134+
}
135+
136+
sub GSL_ODEIV_BAD_Y_INPUT : Tests {
137+
my $self = shift;
138+
139+
my $c = $self->{config};
140+
my $t = 0.0;
141+
my $t1 = 100.0;
142+
my $h = 1e-6;
143+
my $y = 1.0;
144+
145+
throws_ok {
146+
gsl_odeiv_evolve_apply (
147+
$c->{e}, $c->{c}, $c->{s}, $c->{sys}, \$t, $t1, \$h, $y
148+
);
149+
} qr/not a reference/, 'wrong y type 1';
150+
$y = { a => 1 };
151+
throws_ok {
152+
gsl_odeiv_evolve_apply (
153+
$c->{e}, $c->{c}, $c->{s}, $c->{sys}, \$t, $t1, \$h, $y
154+
);
155+
} qr/not an array reference/, 'wrong y type 2';
156+
}
157+
158+
sub func2 {
159+
my ($t, $y, $dydt, $params) = @_;
160+
my $mu = $params->{mu};
161+
162+
$dydt->[0] = $y->[1];
163+
$dydt->[1] = -$y->[0] - $mu*$y->[1]*(($y->[0])**2 - 1);
164+
return $GSL_SUCCESS;
165+
}
166+
167+
sub GSL_ODEIV_BAD_FUNC_CALLBACK : Tests {
168+
my $self = shift;
169+
170+
my $c = $self->{config};
171+
my $t1 = 100.0;
172+
my $h = 1e-6;
173+
my $y = [ 1.0, 0.0 ];
174+
my $t = 0.0;
175+
my $sys = $c->{sys};
176+
$sys->{func} = sub { my ($t, $y, $dydt, $params) = @_; $dydt->[2] = 3; 0 };
177+
178+
throws_ok {
179+
gsl_odeiv_evolve_apply (
180+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
181+
);
182+
} qr/array of wrong dimension/, 'wrong callback return type 1';
183+
$sys->{func} = sub { my ($t, $y, $dydt, $params) = @_;
184+
$dydt->[0]=1; $dydt->[1]={}; 0 };
185+
throws_ok {
186+
gsl_odeiv_evolve_apply (
187+
$c->{e}, $c->{c}, $c->{s}, $sys, \$t, $t1, \$h, $y
188+
);
189+
} qr/array value is not a scalar/, 'wrong callback return type 2';
190+
}
191+
192+
sub GSL_ODEIV_BAD_T_INPUT : Tests {
193+
my $self = shift;
194+
195+
my $c = $self->{config};
196+
my $t1 = 100.0;
197+
my $h = 1e-6;
198+
my $y = [ 1.0, 0.0 ];
199+
my $t = 0.0;
200+
201+
throws_ok {
202+
gsl_odeiv_evolve_apply (
203+
$c->{e}, $c->{c}, $c->{s}, $c->{sys}, $t, $t1, \$h, $y
204+
);
205+
} qr/not a reference/, 'wrong t type 1';
206+
$t = { a => 1 };
207+
throws_ok {
208+
gsl_odeiv_evolve_apply (
209+
$c->{e}, $c->{c}, $c->{s}, $c->{sys}, $t, $t1, \$h, $y
210+
);
211+
} qr/not a scalar reference/, 'wrong t type 2';
212+
}
213+
214+
sub GSL_ODEIV_STEP1 : Tests {
215+
my $self = shift;
216+
217+
my $c = $self->{config};
218+
my $t = 0.0;
219+
my $t1 = 100.0;
220+
my $h = 1e-6;
221+
my $y = [ 1.0, 0.0 ];
222+
223+
my $status = gsl_odeiv_evolve_apply (
224+
$c->{e}, $c->{c}, $c->{s}, $c->{sys}, \$t, $t1, \$h, $y
225+
);
226+
ok_similar( $t, 1e-6, 'first_step_dt', 1e-8 );
227+
ok_similar( $h, 5e-6, 'first_step_h', 1e-8 );
228+
}
229+
230+
Test::Class->runtests;

0 commit comments

Comments
 (0)
Please sign in to comment.