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
|
# File 'ext/ode/ode.c', line 35
VALUE excute_lsoda(VALUE self, VALUE t_out, VALUE func, VALUE jac, VALUE init_t, VALUE init_y, VALUE fargs, VALUE opts){
int itol; double rtol, atol;
int *iwork; double *rwork;
int lrw, liw;
int neq = RARRAY_LEN(init_y);
double *y = (double *)ALLOC_N(double, neq);
double t = NUM2DBL(init_t);
double tout = NUM2DBL(t_out);
int itask = 1; // run solver untill t == tout
int istate = 1; // this is the first call of the problem
int iopt = 0; // no optional inputs
int jt = 2; // jacobian type indicator. 1: user provides full jacobian. 2: interanally generated jacobian.
VALUE ret_arr = rb_ary_new2(neq);
double h0 = 0, hmax =0, hmin = 0;
int ixpr=0, max_step=0, max_hnil=0, max_ordn=12, max_ords=5;
int i;
RESTORE_NEQ(neq)
RESTORE_FUNC(func)
RESTORE_FARGS(fargs)
for(i=0; i<neq; i++){
y[i] = NUM2DBL(rb_ary_entry(init_y, i));
}
// parse options (TODO: accept array)
itol = NUM2INT(rb_hash_lookup(opts, ID2SYM(rb_intern("itol"))));
rtol = NUM2DBL(rb_hash_aref(opts, ID2SYM(rb_intern("rtol"))));
atol = NUM2DBL(rb_hash_aref(opts, ID2SYM(rb_intern("atol"))));
// decide lrw, liw
liw = 20 + neq;
iwork = (int *)ALLOC_N(int, liw);
if(jac != Qnil){
jt = 1;
RESTORE_JAC(jac);
lrw = max(20+16*neq, 22+9*neq+neq*neq);
}
else{
int ml=0, mu=0; VALUE val;
if(val = rb_hash_aref(opts, ID2SYM(rb_intern("ml"))) != Qnil)ml = NUM2INT(val);
if(val = rb_hash_aref(opts, ID2SYM(rb_intern("mu"))) != Qnil)mu = NUM2INT(val);
iwork[0] = mu;
iwork[1] = ml;
lrw = max(20+16*neq, 22+10*neq+(2*ml+mu)*neq);
}
rwork = (double *)ALLOC_N(double, lrw);
rwork[4] = h0;
rwork[5] = hmax;
rwork[6] = hmin;
iwork[4] = ixpr;
iwork[5] = max_step;
iwork[6] = max_hnil;
iwork[7] = max_ordn;
iwork[8] = max_ords;
// run lsoda
lsoda_(ode_function, &neq, y, &t, &tout, &itol, &rtol, &atol, &itask, &istate, &iopt, rwork, &lrw, iwork, &liw, ode_jacobian_function, &jt);
// check error
if(istate<0){
rb_raise(rb_eTypeError, "Something went wrong.");
}
for(i=0; i<neq; i++){
VALUE num = DBL2NUM(y[i]);
rb_ary_store(ret_arr, i, num);
}
xfree(y);
xfree(iwork);
xfree(rwork);
return ret_arr;
}
|