view release on metacpan or  search on metacpan

Lsoda.xs  view on Meta::CPAN

  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