The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
#ifdef __cplusplus
extern "C" {
#endif
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "ppport.h"
#include <math.h>
  typedef int (*deriv_func) (const int* dim, const double* t, double* x, double* y);
  typedef int (*jacobi_func) (const int* dim, const double* t, double* x, const int* ml, const int* mu, double* pd, const int* nrowpd);
  void dlsoda_(deriv_func f, const int* dim, double* x,
               double* t, const double* tout, const int* itol, const double* rtol,
               const double* atol, const int* itask, int* istate,
               const int* iopt, double* rwork, const int* lrw, int* iwork,
               const int* liw, jacobi_func jac, const int* jt);
#ifdef __cplusplus
}
#endif

typedef PerlIO * OutputStream;

#define XS_STATE(type, x)                           \
  INT2PTR(type, SvROK(x) ? SvIV(SvRV(x)) : SvIV(x))

#define XS_STRUCT2OBJ(sv, class, obj)           \
  if (obj == NULL) {                            \
    sv_setsv(sv, &PL_sv_undef);                 \
  } else {                                      \
    sv_setref_pv(sv, class, (void *) obj);      \
  }

static SV* func = (SV*)NULL;

void call_func(SV* func_name, const double* t, SV* x, SV* y)
{
  dSP;
  ENTER;
  SAVETMPS;
  PUSHMARK(SP);
  XPUSHs(sv_2mortal(newSVnv(*t)));
  XPUSHs(sv_2mortal(newSVsv(x)));
  XPUSHs(sv_2mortal(newSVsv(y)));
  PUTBACK;
  call_sv(func_name, G_ARRAY);
  SPAGAIN;
  PUTBACK;
  FREETMPS;
  LEAVE;
}

int equation(const int* dim, const double* t, double* x, double* y)
{
  int size = *dim, i;
  SV* sv_x[size];
  SV* sv_y[size];
  for(i=0;i<size;i++){
    sv_x[i] = sv_2mortal(newSVnv(x[i]));
    sv_y[i] = sv_2mortal(newSVnv(y[i]));
  }
  SV* sv_x1 = newRV_noinc((SV*)av_make(size, sv_x));
  SV* sv_y1 = newRV_noinc((SV*)av_make(size, sv_y));

  call_func(func, t, sv_x1, sv_y1);
  AV* av_y = (AV*)SvRV(sv_y1);
  for(i=0;i<size;i++){
    SV** pv = av_fetch(av_y,i,0);
    y[i] = SvNV(*pv);
  }
}

int required_size(double start, double end, double dt)
{
  if(dt!=0){
    int size = (int)ceil( ((end - start)/dt) ) + 1;
    if(size > 0) return size;
  }
  return -1;
}

int
solve(SV* func_name, AV* x, double t, double tout, double dt, AV* rtol, AV* atol, OutputStream stream)
{
  func = func_name;
  int dim = av_len(x) + 1;
  int i, step;
  int maxvalue = 16 > (dim+9) ? 16 : (dim+9);
  int lrw = 22 + dim * maxvalue;
  int liw = 20 + dim;
  int itol = 2;
  int itask = 1;
  int istate = 1;
  int iopt = 0;
  int jt = 2;
  double a_x[dim], a_rtol[dim], a_atol[dim];
  double rwork[lrw];
  int iwork[liw];
  FILE *fp = PerlIO_findFILE(stream);
  int maxStep = required_size(t, tout, dt);
  double t1 = t + dt;
  for(i=0;i<dim;i++){
    SV** pv = av_fetch(x,i,0);
    a_x[i] = SvNV(*pv);
    pv = av_fetch(rtol,i,0);
    a_rtol[i] = SvNV(*pv);
    pv= av_fetch(atol,i,0);
    a_atol[i] = SvNV(*pv);
  }

  fprintf(fp, "%g", t);
  for(i=0;i<dim;i++){
    fprintf(fp, "\t%g", a_x[i]);
  }
  fprintf(fp, "\n");

  for(step=1;step<maxStep;step++){
    dlsoda_(equation, &dim, a_x, &t, &t1, &itol, a_rtol, a_atol, &itask, &istate, &iopt, rwork, &lrw, iwork, &liw, NULL, &jt);

    fprintf(fp, "%g", t);
    for(i=0;i<dim;i++){
    fprintf(fp, "\t%g", a_x[i]);
    }
    fprintf(fp, "\n");

    t1 = t1 + dt;
    switch (istate) {
    case 2: /* successful exit */
      break;
    case -1:
      fprintf(stderr, "[WARNING]: excess work done at t =%14.6e (perhaps wrong JT).\n", t);
      /* istate = 1; */
      return istate;
    case -2:
      fprintf(stderr, "[ERROR]: excess accuracy requested at t =%14.6e (tolerances too small).\n", t);
      return istate;
    case -3:
      fprintf(stderr, "[ERROR]: illegal input detected at t =%14.6e (see printed message).\n", t);
      return istate;
    case -4:
      fprintf(stderr, "[ERROR]: repeated error test failures at t =%14.6e (check all inputs).\n", t);
      return istate;
    case -5:
      fprintf(stderr, "[ERROR]: repeated convergence failures at t =%14.6e (perhaps bad Jacobian supplied or wrong choice of JT or tolerances).\n", t);
      return istate;
    case -6:
      fprintf(stderr, "[ERROR]: error weight became zero at t =%14.6e (Solution component i vanished, and ATOL or ATOL(i) = 0.)\n", t);
      return istate;
    case -7:
      fprintf(stderr, "[ERROR]: work space insufficient at t =%14.6e\n", t);
      return istate;
    default:
      fprintf(stderr, "[ERROR]: unknown error at t =%14.6e\n", t);
      return istate;
    }
  }
  return istate;
}

MODULE = Math::Lsoda  PACKAGE = Math::Lsoda

int
solve(func_name, x, t, tout, dt, rtol, atol, stream)
SV* func_name;
AV* x;
double t;
double tout;
double dt;
AV* rtol;
AV* atol;
OutputStream stream;