%module "Math::GSL::ODEIV"
%include "gsl_typemaps.i"
%include "typemaps.i"
%include "renames.i"
//%newobject gsl_odeiv_step_alloc;
%newobject gsl_odeiv_driver_alloc_y_new;
//%newobject gsl_odeiv_evolve_alloc;
%{
#include "gsl/gsl_types.h"
#include "gsl/gsl_odeiv.h"
#include <stdarg.h>
#define SWIG_MATH_GSL_ODEIV_PACKAGE_NAME "Math::GSL::ODEIV"
#define SWIG_MATH_GSL_ODEIV_GUTS \
SWIG_MATH_GSL_ODEIV_PACKAGE_NAME "::_guts"
/* NOTE: We do not use C static variables to store the values since
static variable are not thread-safe, and SWIG does not currently support the
MY_CXT framework, see perldoc perlxs for more information.
Perl package variables are on the other hand thread safe.
*/
void swig_math_gsl_odeiv_set_gut_pv(
const char *varname, const char *value
)
{
SV *pvname;
SV *sv;
pvname = newSVpvf( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);
/*char *pvname = form( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);*/
sv = get_sv( SvPV_nolen(pvname), GV_ADD );
SvREFCNT_dec(pvname);
sv_setpv( sv, value );
}
const char *swig_math_gsl_odeiv_get_gut_pv(const char *varname)
{
SV *pvname, *sv;
pvname = newSVpvf( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);
/*char *pvname = form( "%s::%s", SWIG_MATH_GSL_ODEIV_GUTS, varname);*/
sv = get_sv( SvPV_nolen(pvname), GV_ADD );
SvREFCNT_dec(pvname);
return SvPV_nolen(sv);
}
void swig_math_gsl_odeiv_set_callback_error_param( const char *func )
{
swig_math_gsl_odeiv_set_gut_pv( "cbname", func );
}
void swig_math_gsl_odeiv_set_error_param(
const char *symname, const char *param
)
{
swig_math_gsl_odeiv_set_gut_pv( "symname", symname );
swig_math_gsl_odeiv_set_gut_pv( "param", param );
}
void swig_math_gsl_odeiv_callback_error(
const char *msg, ...
)
{
const char *cbname;
va_list args;
SV *msg2;
cbname = swig_math_gsl_odeiv_get_gut_pv("cbname");
va_start(args, msg);
/*char *msg2 = form(
"Math::GSL::ODEIV: callback function : %s() : %s", cbname, msg
);*/
msg2 = newSVpvf(
"Math::GSL::ODEIV: callback function : %s() : %s", cbname, msg
);
vcroak( SvPV_nolen(msg2), &args );
/* NOTE: these two lines will never be reached */
SvREFCNT_dec(msg2);
va_end(args);
}
void swig_math_gsl_odeiv_input_param_error(
const char *msg, ...
)
{
const char *subname;
const char *param;
va_list args;
SV *msg2;
subname = swig_math_gsl_odeiv_get_gut_pv("symname");
param = swig_math_gsl_odeiv_get_gut_pv("param");
va_start(args, msg);
msg2 = newSVpvf(
"Math::GSL::ODEIV:%s() : parameter $%s : %s", subname, param, msg
);
vcroak( SvPV_nolen(msg2), &args );
/* NOTE: these two lines will never be reached */
SvREFCNT_dec(msg2);
va_end(args);
}
void swig_math_gsl_odeiv_input_error(
const char *msg, ...
)
{
const char *subname;
SV *msg2;
va_list args;
subname = swig_math_gsl_odeiv_get_gut_pv("symname");
va_start(args, msg);
msg2 = newSVpvf( "Math::GSL::ODEIV:%s() : %s", subname, msg);
vcroak( SvPV_nolen(msg2), &args );
/* NOTE: these two lines will never be reached */
SvREFCNT_dec(msg2);
va_end(args);
}
SV *swig_math_gsl_odeiv_get_hash_sv(HV *hash, const char *key)
{
SV *key_sv, *value;
HE *he;
key_sv = newSVpv(key, strlen (key));
if (hv_exists_ent(hash, key_sv, 0)) {
he = hv_fetch_ent(hash, key_sv, 0, 0);
value = HeVAL(he);
}
else {
swig_math_gsl_odeiv_input_param_error(
"The hash key '%s' doesn't exist", key
);
}
return value;
}
IV swig_math_gsl_odeiv_get_hash_iv(HV *hash, const char *key) {
SV *sv;
sv = swig_math_gsl_odeiv_get_hash_sv(hash, key);
if (SvROK(sv)) {
swig_math_gsl_odeiv_input_param_error(
"Hash value for key '%s' is not a scalar value", key
);
}
if (!SvIOK(sv)) {
swig_math_gsl_odeiv_input_param_error(
"Hash value for key '%s' is not an integer", key
);
}
return SvIV(sv);
}
SV *swig_math_gsl_odeiv_get_hash_hashref(HV *hash, const char *key) {
SV *sv;
sv = swig_math_gsl_odeiv_get_hash_sv(hash, key);
if (!SvROK(sv)) {
swig_math_gsl_odeiv_input_param_error(
"Hash value for key '%s' is not a reference", key
);
}
if (SvTYPE(SvRV(sv)) != SVt_PVHV) {
swig_math_gsl_odeiv_input_param_error(
"Hash value for key '%s' is not a hashref", key
);
}
return sv;
}
SV *swig_math_gsl_odeiv_get_hash_coderef(HV *hash, const char *key) {
SV *sv;
sv = swig_math_gsl_odeiv_get_hash_sv(hash, key);
if (!SvROK(sv)) {
swig_math_gsl_odeiv_input_param_error(
"Hash value for key '%s' is not a reference", key
);
}
if (SvTYPE(SvRV(sv)) != SVt_PVCV) {
swig_math_gsl_odeiv_input_param_error(
"Hash value for key '%s' is not a coderef", key
);
}
return sv;
}
void swig_math_gsl_odeiv_store_hash_ptr( HV *hash, const char *key, void *ptr)
{
SV *sv;
sv = newSViv(PTR2IV(ptr));
/* Let the hash take ownership of the sv */
if( !hv_store(hash, key, strlen(key), sv, 0) ) {
SvREFCNT_dec(sv);
swig_math_gsl_odeiv_input_param_error(
"Cannot store internal information in the hash"
);
}
}
void *swig_math_gsl_odeiv_get_hash_ptr(HV *hash, const char *key) {
IV ptr;
ptr = swig_math_gsl_odeiv_get_hash_iv(hash, key);
return (void *) INT2PTR(SV*, ptr);
}
void swig_math_gsl_odeiv_store_double_in_av( AV *array, SSize_t index, double val)
{
SV *sval;
sval = newSVnv(val);
if( !av_store(array, index, sval) ) {
SvREFCNT_dec(sval);
swig_math_gsl_odeiv_callback_error(
"Cannot store internal information in array"
);
}
}
void swig_math_gsl_odeiv_copy_av_to_carray(AV *array, double *y, size_t dim)
{
int i;
SSize_t array_len;
SV **sv_ptr;
SV *sv;
double val;
array_len = av_len(array) + 1;
if (array_len != dim ) {
swig_math_gsl_odeiv_callback_error(
"Callback returned array of wrong dimension"
);
}
for (i = 0; i < dim; i++) {
sv_ptr = av_fetch( array, i, 0 );
if (!sv_ptr) {
swig_math_gsl_odeiv_callback_error(
"Cannot extract values from returned array"
);
}
sv = *sv_ptr;
if (SvROK(sv)) {
swig_math_gsl_odeiv_callback_error(
"Returned array value is not a scalar"
);
}
/* TODO: Not sure if this check is necessary */
if (SvTYPE(sv) >= SVt_PVAV) {
swig_math_gsl_odeiv_callback_error(
"Returned array value is not of scalar type"
);
}
val = (double ) SvNV(sv);
y[i] = val;
}
}
void swig_math_gsl_odeiv_copy_doubles_to_av(AV *array, const double *y, size_t dim)
{
int i;
for (i = 0; i < dim; i++) {
swig_math_gsl_odeiv_store_double_in_av(array, i, y[i]);
}
}
typedef struct {
SV *func;
SV *jac;
SV *params;
size_t dim;
gsl_odeiv_system *sys;
} swig_math_gsl_odeiv_system;
int swig_math_gsl_odeiv_call_perl_jac (
SV *callback, double t, const double y[], double *dfdy, double *dfdt,
swig_math_gsl_odeiv_system *params
)
{
AV *ay, *a_dfdy, *a_dfdt;
int count;
IV result;
ay = (AV *)sv_2mortal((SV *)newAV());
a_dfdy = (AV *)sv_2mortal((SV *)newAV());
a_dfdt = (AV *)sv_2mortal((SV *)newAV());
dSP; /* declares a local copy of stack pointer */
ENTER;
SAVETMPS;
PUSHMARK(SP);
EXTEND(SP, 5);
PUSHs(sv_2mortal(newSVnv(t)));
swig_math_gsl_odeiv_copy_doubles_to_av(ay, y, params->dim);
PUSHs(sv_2mortal((SV *)newRV_inc((SV *) ay)));
PUSHs(sv_2mortal((SV *)newRV_inc((SV *) a_dfdy)));
PUSHs(sv_2mortal((SV *)newRV_inc((SV *) a_dfdt)));
PUSHs(params->params);
PUTBACK;
count = call_sv(callback, G_SCALAR); /* call the Perl callback */
SPAGAIN;
if (count != 1) {
swig_math_gsl_odeiv_callback_error(
"Bad return value from callback: expected 1 value, got %d", count
);
}
result = POPi; /* TODO: check ST(0) instead for valid value */
swig_math_gsl_odeiv_copy_av_to_carray(a_dfdy, dfdy, (params->dim)*(params->dim));
swig_math_gsl_odeiv_copy_av_to_carray(a_dfdt, dfdt, params->dim);
PUTBACK;
FREETMPS;
LEAVE;
return (int) result;
}
int swig_math_gsl_odeiv_call_perl_func (
SV *callback, double t, const double y[], double dydt[],
swig_math_gsl_odeiv_system *params)
{
AV *ay, *aj;
int count;
IV result;
ay = (AV *)sv_2mortal((SV *)newAV());
aj = (AV *)sv_2mortal((SV *)newAV());
dSP; /* declares a local copy of stack pointer */
ENTER;
SAVETMPS;
PUSHMARK(SP);
EXTEND(SP, 4);
PUSHs(sv_2mortal(newSVnv(t)));
swig_math_gsl_odeiv_copy_doubles_to_av(ay, y, params->dim);
PUSHs(sv_2mortal((SV *)newRV_inc((SV *) ay)));
PUSHs(sv_2mortal((SV *)newRV_inc((SV *) aj)));
PUSHs(params->params);
PUTBACK;
count = call_sv(callback, G_SCALAR); /* call the Perl callback */
SPAGAIN;
/* This should not happen for G_SCALAR, see perldoc perlcall.
* Even if the callback does not return anything, count will still be 1
* since we are not passing the G_DISCARD flag
*/
if (count != 1) {
swig_math_gsl_odeiv_callback_error(
"Bad return value from callback: expected 1 value, got %d", count
);
}
result = POPi; /* TODO: check ST(0) instead for valid value */
swig_math_gsl_odeiv_copy_av_to_carray(aj, dydt, params->dim);
PUTBACK;
FREETMPS;
LEAVE;
return result;
}
int swig_math_gsl_odeiv_callback_function(
double t, const double y[], double dydt[], void *params)
{
swig_math_gsl_odeiv_set_callback_error_param( "func" );
swig_math_gsl_odeiv_system *sparam = (swig_math_gsl_odeiv_system *) params;
return swig_math_gsl_odeiv_call_perl_func(sparam->func, t, y, dydt, sparam);
}
int swig_math_gsl_odeiv_callback_jac(
double t, const double y[], double *dfdy, double *dfdt, void *params
)
{
swig_math_gsl_odeiv_set_callback_error_param( "jac" );
swig_math_gsl_odeiv_system *sparam = (swig_math_gsl_odeiv_system *) params;
return swig_math_gsl_odeiv_call_perl_jac(sparam->jac, t, y, dfdy, dfdt, sparam);
}
void swig_math_gsl_odeiv_fill_system_struct( HV *hash, gsl_odeiv_system *sys )
{
swig_math_gsl_odeiv_system *ssys;
Newx(ssys, 1, swig_math_gsl_odeiv_system);
ssys->func = swig_math_gsl_odeiv_get_hash_coderef( hash, "func" );
ssys->jac = swig_math_gsl_odeiv_get_hash_coderef( hash, "jac" );
ssys->params = swig_math_gsl_odeiv_get_hash_hashref( hash, "params" );
ssys->dim = (size_t) swig_math_gsl_odeiv_get_hash_iv( hash, "dim" );
sys->function = swig_math_gsl_odeiv_callback_function;
sys->jacobian = swig_math_gsl_odeiv_callback_jac;
sys->dimension = (size_t) swig_math_gsl_odeiv_get_hash_iv( hash, "dim" );
sys->params = (void *) ssys;
}
%}
%typemap(in) const gsl_odeiv_system * {
SV *sv;
HV *hash;
swig_math_gsl_odeiv_set_error_param( "$symname", "$1_name" );
if (!sv_isobject($input)) {
swig_math_gsl_odeiv_input_error(
"Input parameter $$1_name is not a blessed reference!"
);
}
if (!sv_isa($input, "Math::GSL::ODEIV::gsl_odeiv_system")) {
swig_math_gsl_odeiv_input_error(
"Input parameter $$1_name is not an object of type "
"Math::GSL::ODEIV::gsl_odeiv_system!"
);
}
sv = SvRV($input);
if ((SvTYPE(sv) != SVt_PVHV)) {
swig_math_gsl_odeiv_input_error(
"Input parameter $$1_name is not a blessed hash reference!"
);
}
gsl_odeiv_system *sys;
Newx(sys, 1, gsl_odeiv_system);
swig_math_gsl_odeiv_fill_system_struct( (HV *)sv, sys);
$1 = sys;
}
%typemap(freearg) const gsl_odeiv_system * {
swig_math_gsl_odeiv_system *ssys = (swig_math_gsl_odeiv_system *)$1->params;
Safefree(ssys);
Safefree($1);
}
%ignore gsl_odeiv_evolve_apply;
%include "gsl/gsl_types.h"
%include "gsl/gsl_odeiv.h"
// unignore gsl_odeiv_evolve_apply()
%rename("%s") gsl_odeiv_evolve_apply;
/* We want handle the last parameter to gsl_odeiv_evolve_apply(...), see
* include file gsl/gsl_odeiv.h for a definition. This parameter is of
* type 'double []' and acts as an input-output array.
*
* NOTE: gsl_typemaps.i has typemaps for a float [] input-output array,
* but note that that typemap also returns the array elements on the perl stack
* (in addition to modifying the original array).
* However, here we do not want to return the result on the stack, we only
* want to modify the original array.
*
* TODO: These typemaps might warrant to be moved to gsl_typemaps.i at a later time,
* where they could be reused by other interface files, however currently they are
* regarded as specific to only gsl_odeiv_evolve_apply().
*/
%typemap(in) double y[] {
struct perl_array * p_array = 0; /* see gsl_typemaps.i for definition */
I32 len;
AV *array;
int i;
SV **tv;
swig_math_gsl_odeiv_set_error_param( "$symname", "$1_name" );
if (!SvROK($input)) {
swig_math_gsl_odeiv_input_error(
"Input parameter $$1_name is not a reference!"
);
}
if (SvTYPE(SvRV($input)) != SVt_PVAV) {
swig_math_gsl_odeiv_input_error(
"Input parameter $$1_name is not an array reference!"
);
}
array = (AV*)SvRV($input);
len = av_len(array);
p_array = (struct perl_array *)
malloc((len+1)*sizeof(double)+sizeof(struct perl_array));
p_array->len=len;
p_array->array=array;
$1 = (double *)&p_array[1];
for (i = 0; i <= len; i++) {
tv = av_fetch(array, i, 0);
$1[i] = (double) SvNV(*tv);
}
}
%typemap(argout) double y[] {
struct perl_array * p_array = 0; /* see gsl_typemaps.i for definition */
int i;
double val;
SV **tv;
p_array=(struct perl_array *)(((char*)$1)-sizeof(struct perl_array));
for (i = 0; i <= p_array->len; i++) {
val = $1[i];
tv = av_fetch(p_array->array, i, 0);
sv_setnv(*tv, val);
}
}
%typemap(freearg) double y[] {
if ($1) free(((char*)$1)-sizeof(struct perl_array));
}
%{
typedef struct {
double h;
SV *sv;
} swig_perl_double_in_out;
%}
%typemap(in) double *h {
SV *sv;
swig_perl_double_in_out *h_wrap;
swig_math_gsl_odeiv_set_error_param( "$symname", "$1_name" );
if (!SvROK($input)) {
swig_math_gsl_odeiv_input_error(
"Input parameter $$1_name is not a reference!"
);
}
if (SvTYPE(SvRV($input)) >= SVt_PVAV) {
swig_math_gsl_odeiv_input_error(
"Input parameter $$1_name is not a scalar reference!"
);
}
sv = SvRV($input);
Newx(h_wrap, 1, swig_perl_double_in_out);
h_wrap->sv = sv;
h_wrap->h = (double) SvNV(sv);
$1 = (double *) &h_wrap->h;
}
%typemap(argout) double *h {
swig_perl_double_in_out *h_wrap;
SV *sv;
h_wrap = (swig_perl_double_in_out *) $1;
sv = h_wrap->sv;
sv_setnv(sv, h_wrap->h);
}
%typemap(freearg) double *h {
Safefree( $1 );
}
%typemap(in) double *t = double *h;
%typemap(argout) double *t = double *h;
%typemap(freearg) double *t = double *h;
// define our own name for the last parameter to gsl_odeiv_evolve_apply()
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[]);
%perlcode %{
package Math::GSL::ODEIV::gsl_odeiv_system;
sub new {
my ($class, $func, $jac, $dim, $params ) = @_;
my $check_ref = sub {
if ( (ref $_[0]) ne $_[1] ) {
die sprintf 'Usage: %s:new( $func, $jac, $dim, $params ). '
. 'Argument %s is not %s reference',
__PACKAGE__, $_[2], $_[3];
}
};
my $check_subref = sub {
$check_ref->($_[0], "CODE", $_[1], "code");
};
my $check_hashref = sub {
$check_ref->($_[0], "HASH", $_[1], "hash");
};
$check_subref->($func, '$func');
$check_subref->($jac, '$jac');
$check_hashref->($params, '$params');
return bless { func => $func, jac => $jac, dim => $dim, params => $params },
$class;
}
%}
%include "../pod/ODEIV.pod"