-
Notifications
You must be signed in to change notification settings - Fork 0
/
Lsoda.xs
169 lines (154 loc) · 4.56 KB
/
Lsoda.xs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#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;