← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/integrate/__odepack.h

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/integrate/__odepack.h
1
 
/* This file should be included in the multipack module */
2
 
/* $Revision: 2055 $ */
3
 
/*  module_methods:
4
 
 {"odeint", (PyCFunction) odepack_odeint, METH_VARARGS|METH_KEYWORDS, doc_odeint},
5
 
 */
6
 
/* link libraries: (should be listed in separate lines)
7
 
   odepack
8
 
   linpack_lite
9
 
   blas
10
 
   mach
11
 
 */
12
 
/* python files: (to be appended to Multipack.py)
13
 
   odepack.py
14
 
 */
15
 
 
16
 
 
17
 
#if defined(NO_APPEND_FORTRAN)
18
 
#define LSODA  lsoda
19
 
#else
20
 
#define LSODA  lsoda_
21
 
#endif
22
 
 
23
 
void LSODA();
24
 
 
25
 
/*
26
 
void ode_function(int *n, double *t, double *y, double *ydot)
27
 
{
28
 
  ydot[0] = -0.04*y[0] + 1e4*y[1]*y[2];
29
 
  ydot[2] = 3e7*y[1]*y[1];
30
 
  ydot[1] = -ydot[0] - ydot[2];
31
 
  return;
32
 
}
33
 
*/
34
 
 
35
 
void ode_function(int *n, double *t, double *y, double *ydot)
36
 
{
37
 
 /* This is the function called from the Fortran code it should
38
 
        -- use call_python_function to get a multiarrayobject result
39
 
        -- check for errors and return -1 if any
40
 
        -- otherwise place result of calculation in ydot
41
 
  */
42
 
 
43
 
  PyArrayObject *result_array = NULL;
44
 
  PyObject *arg1, *arglist;
45
 
 
46
 
  /* Append t to argument list */
47
 
  if ((arg1 = PyTuple_New(1)) == NULL) {
48
 
    if (PyErr_Occurred())
49
 
      PyErr_Print();
50
 
    return;
51
 
  }
52
 
  PyTuple_SET_ITEM(arg1, 0, PyFloat_FromDouble(*t)); 
53
 
                /* arg1 now owns newly created reference */
54
 
  if ((arglist = PySequence_Concat( arg1, multipack_extra_arguments)) == NULL) {
55
 
    if (PyErr_Occurred())
56
 
      PyErr_Print();
57
 
    Py_DECREF(arg1);
58
 
    return;
59
 
  }
60
 
  Py_DECREF(arg1);    /* arglist has reference */
61
 
  
62
 
  result_array = (PyArrayObject *)call_python_function(multipack_python_function, *n, y, arglist, 1, odepack_error);
63
 
  if (result_array == NULL) {
64
 
    PyErr_Print();
65
 
    Py_DECREF(arglist);
66
 
    return;
67
 
  }
68
 
  memcpy(ydot, result_array->data, (*n)*sizeof(double));
69
 
  Py_DECREF(result_array);
70
 
  Py_DECREF(arglist);
71
 
  return;
72
 
}
73
 
 
74
 
 
75
 
int ode_jacobian_function(int *n, double *t, double *y, int *ml, int *mu, double *pd, int *nrowpd)
76
 
{
77
 
  /* This is the function called from the Fortran code it should
78
 
        -- use call_python_function to get a multiarrayobject result
79
 
        -- check for errors and return -1 if any (though this is ignored 
80
 
                                 by calling program).
81
 
        -- otherwise place result of calculation in pd
82
 
  */
83
 
 
84
 
  PyArrayObject *result_array;
85
 
  PyObject *arglist, *arg1;
86
 
 
87
 
  /* Append t to argument list */
88
 
  if ((arg1 = PyTuple_New(1)) == NULL) {
89
 
    if (PyErr_Occurred())
90
 
      PyErr_Print();
91
 
    return -1;
92
 
  }
93
 
  PyTuple_SET_ITEM(arg1, 0, PyFloat_FromDouble(*t)); 
94
 
                /* arg1 now owns newly created reference */
95
 
  if ((arglist = PySequence_Concat( arg1, multipack_extra_arguments)) == NULL) {
96
 
    if (PyErr_Occurred())
97
 
      PyErr_Print();
98
 
    Py_DECREF(arg1);
99
 
    return -1;
100
 
  }
101
 
  Py_DECREF(arg1);    /* arglist has reference */
102
 
 
103
 
  result_array = (PyArrayObject *)call_python_function(multipack_python_jacobian, *n, y, arglist, 2, odepack_error);
104
 
  if (result_array == NULL) {
105
 
    Py_DECREF(arglist);
106
 
    return -1;
107
 
  }
108
 
  if (multipack_jac_transpose == 1) 
109
 
    MATRIXC2F(pd, result_array->data, *n, *nrowpd)
110
 
  else
111
 
    memcpy(pd, result_array->data, (*n)*(*nrowpd)*sizeof(double));
112
 
 
113
 
  Py_DECREF(arglist);
114
 
  Py_DECREF(result_array);
115
 
  return 0;
116
 
}
117
 
 
118
 
 
119
 
int setup_extra_inputs(PyArrayObject **ap_rtol, PyObject *o_rtol, PyArrayObject **ap_atol, PyObject *o_atol, PyArrayObject **ap_tcrit, PyObject *o_tcrit, int *numcrit, int neq)
120
 
{
121
 
  int itol = 0;
122
 
  double tol=1.49012e-8;
123
 
  int one = 1;
124
 
 
125
 
  /* Setup tolerances */
126
 
  if (o_rtol == NULL) {
127
 
    *ap_rtol = (PyArrayObject *)PyArray_FromDims(1,&one,PyArray_DOUBLE);
128
 
    if (*ap_rtol == NULL) PYERR2(odepack_error,"Error constructing relative tolerance.");
129
 
    *(double *)(*ap_rtol)->data = tol;                /* Default */
130
 
  }
131
 
  else {
132
 
    *ap_rtol = (PyArrayObject *)PyArray_ContiguousFromObject(o_rtol,PyArray_DOUBLE,0,1);
133
 
    if (*ap_rtol == NULL) PYERR2(odepack_error,"Error converting relative tolerance.");
134
 
    if ((*ap_rtol)->nd == 0); /* rtol is scalar */
135
 
    else if ((*ap_rtol)->dimensions[0] == neq)
136
 
      itol |= 2;      /* Set rtol array flag */
137
 
    else
138
 
      PYERR(odepack_error,"Tolerances must be an array of the same length as the\n     number of equations or a scalar.");
139
 
  }
140
 
 
141
 
  if (o_atol == NULL) {
142
 
    *ap_atol = (PyArrayObject *)PyArray_FromDims(1,&one,PyArray_DOUBLE);
143
 
    if (*ap_atol == NULL) PYERR2(odepack_error,"Error constructing absolute tolerance");
144
 
    *(double *)(*ap_atol)->data = tol;
145
 
  }
146
 
  else {
147
 
    *ap_atol = (PyArrayObject *)PyArray_ContiguousFromObject(o_atol,PyArray_DOUBLE,0,1);
148
 
    if (*ap_atol == NULL) PYERR2(odepack_error,"Error converting absolute tolerance.");
149
 
    if ((*ap_atol)->nd == 0); /* atol is scalar */
150
 
    else if ((*ap_atol)->dimensions[0] == neq) 
151
 
      itol |= 1;        /* Set atol array flag */
152
 
    else
153
 
      PYERR(odepack_error,"Tolerances must be an array of the same length as the\n     number of equations or a scalar.");
154
 
  }
155
 
  itol++;             /* increment to get correct value */
156
 
 
157
 
 
158
 
  /* Setup t-critical */
159
 
  if (o_tcrit != NULL) {
160
 
    *ap_tcrit = (PyArrayObject *)PyArray_ContiguousFromObject(o_tcrit,PyArray_DOUBLE,0,1);
161
 
    if (*ap_tcrit == NULL) PYERR2(odepack_error,"Error constructing critical times.");
162
 
    *numcrit = PyArray_Size((PyObject *)(*ap_tcrit));
163
 
  }
164
 
  return itol;
165
 
 
166
 
 fail:       /* Needed for use of PYERR */
167
 
  return -1;
168
 
}
169
 
 
170
 
 
171
 
int compute_lrw_liw(int *lrw, int *liw, int neq, int jt, int ml, int mu, int mxordn, int mxords)
172
 
{
173
 
  int lrn, lrs, nyh, lmat;
174
 
 
175
 
  if (jt == 1 || jt == 2)
176
 
    lmat = neq*neq + 2;
177
 
  else if (jt == 4 || jt == 5)
178
 
    lmat = (2*ml + mu + 1)*neq + 2;
179
 
  else PYERR(odepack_error,"Incorrect value for jt");
180
 
 
181
 
  if (mxordn < 0) PYERR(odepack_error,"Incorrect value for mxordn");
182
 
  if (mxords < 0) PYERR(odepack_error,"Incorrect value for mxords");
183
 
  nyh = neq;
184
 
 
185
 
  lrn = 20 + nyh*(mxordn+1) + 3*neq;
186
 
  lrs = 20 + nyh*(mxords+1) + 3*neq + lmat;
187
 
 
188
 
  *lrw = NPY_MAX(lrn,lrs);
189
 
  *liw = 20 + neq;
190
 
  return 0;
191
 
 
192
 
 fail:
193
 
  return -1;
194
 
}
195
 
 
196
 
static char doc_odeint[] = "[y,{infodict,}istate] = odeint(fun, y0, t, args=(), Dfun=None, col_deriv=0, ml=, mu=, full_output=0, rtol=, atol=, tcrit=, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0.0, mxstep=0.0, mxhnil=0, mxordn=0, mxords=0)\n  yprime = fun(y,t,...)";
197
 
 
198
 
static PyObject *odepack_odeint(PyObject *dummy, PyObject *args, PyObject *kwdict) {
199
 
  PyObject *fcn, *y0, *p_tout, *o_rtol=NULL, *o_atol=NULL;
200
 
  PyArrayObject *ap_y = NULL, *ap_yout= NULL;
201
 
  PyArrayObject *ap_rtol=NULL, *ap_atol=NULL;
202
 
  PyArrayObject *ap_tout = NULL;
203
 
  PyObject *extra_args = NULL;
204
 
  PyObject *Dfun = Py_None;
205
 
  int      neq, itol=1, itask=1, istate=1, iopt=0, lrw, *iwork, liw, jt=4;
206
 
  double   *y, t, *tout, *rtol, *atol, *rwork;
207
 
  double   h0=0.0, hmax=0.0, hmin=0.0;
208
 
  int      ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, ml=-1, mu=-1;
209
 
  PyObject *o_tcrit=NULL;
210
 
  PyArrayObject *ap_tcrit=NULL;
211
 
  PyArrayObject *ap_hu=NULL, *ap_tcur=NULL, *ap_tolsf=NULL, *ap_tsw=NULL;
212
 
  PyArrayObject *ap_nst=NULL, *ap_nfe=NULL, *ap_nje=NULL, *ap_nqu=NULL;
213
 
  PyArrayObject *ap_mused=NULL;
214
 
  int      imxer=0, lenrw=0, leniw=0, out_sz=0, col_deriv = 0;
215
 
  int      k, dims[2], ntimes, crit_ind=0;
216
 
  int      allocated = 0, full_output = 0, numcrit=0;
217
 
  double   *yout, *yout_ptr, *tout_ptr, *tcrit;
218
 
  double   *wa;
219
 
  static char *kwlist[] = {"fun","y0","t","args","Dfun","col_deriv","ml","mu","full_output","rtol","atol","tcrit","h0","hmax","hmin","ixpr","mxstep","mxhnil","mxordn","mxords",NULL};
220
 
 
221
 
  STORE_VARS();
222
 
 
223
 
  if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OOiiiiOOOdddiiiii", kwlist, &fcn, &y0, &p_tout, &extra_args, &Dfun, &col_deriv, &ml, &mu, &full_output, &o_rtol, &o_atol, &o_tcrit, &h0, &hmax, &hmin, &ixpr, &mxstep, &mxhnil, &mxordn, &mxords)) return NULL;
224
 
 
225
 
  if (o_tcrit == Py_None) {
226
 
    o_tcrit = NULL;
227
 
  }
228
 
  if (o_rtol == Py_None) {
229
 
    o_rtol = NULL;
230
 
  }
231
 
  if (o_atol == Py_None) {
232
 
    o_atol = NULL;
233
 
  }
234
 
 
235
 
 
236
 
  INIT_JAC_FUNC(fcn,Dfun,extra_args,col_deriv,odepack_error);
237
 
 
238
 
  /* Set up jt, ml, and mu */
239
 
  if (Dfun == Py_None) jt++;    /* set jt for internally generated */
240
 
  if (ml < 0 && mu < 0) jt -= 3;     /* neither ml nor mu given, 
241
 
                                        mark jt for full jacobian */
242
 
  if (ml < 0) ml = 0;    /* if one but not both are given */
243
 
  if (mu < 0) mu = 0;
244
 
 
245
 
  /* Initial input vector */
246
 
  ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y0, PyArray_DOUBLE, 0, 1);
247
 
  if (ap_y == NULL) goto fail;
248
 
  y = (double *) ap_y->data;
249
 
  neq = PyArray_Size((PyObject *)ap_y);
250
 
  dims[1] = neq;
251
 
 
252
 
  /* Set of output times for integration */
253
 
  ap_tout = (PyArrayObject *)PyArray_ContiguousFromObject(p_tout, PyArray_DOUBLE, 0, 1);
254
 
  if (ap_tout == NULL) goto fail;
255
 
  tout = (double *)ap_tout->data;
256
 
  ntimes = PyArray_Size((PyObject *)ap_tout);
257
 
  dims[0] = ntimes;
258
 
  t = tout[0];
259
 
 
260
 
  /* Setup array to hold the output evaluations*/
261
 
  ap_yout= (PyArrayObject *)PyArray_FromDims(2,dims,PyArray_DOUBLE);
262
 
  if (ap_yout== NULL) goto fail;
263
 
  yout = (double *) ap_yout->data;
264
 
  /* Copy initial vector into first row of output */
265
 
  memcpy(yout, y, neq*sizeof(double));  /* copy intial value to output */
266
 
  yout_ptr = yout + neq;    /* set output pointer to next position */
267
 
 
268
 
  itol = setup_extra_inputs(&ap_rtol, o_rtol, &ap_atol, o_atol, &ap_tcrit, o_tcrit, &numcrit, neq);
269
 
  if (itol < 0 ) goto fail;  /* Something didn't work */
270
 
  rtol = (double *) ap_rtol->data;
271
 
  atol = (double *) ap_atol->data;
272
 
  if (o_tcrit != NULL) tcrit = (double *)(ap_tcrit->data);
273
 
 
274
 
  /* Find size of working arrays*/
275
 
  if (compute_lrw_liw(&lrw, &liw, neq, jt, ml, mu, mxordn, mxords) < 0) goto fail;
276
 
 
277
 
  if ((wa = (double *)malloc(lrw*sizeof(double) + liw*sizeof(int)))==NULL) {
278
 
    PyErr_NoMemory();
279
 
    goto fail;
280
 
  }
281
 
  allocated = 1;
282
 
  rwork = wa;
283
 
  iwork = (int *)(wa + lrw);
284
 
 
285
 
  iwork[0] = ml; iwork[1] = mu;      /* ignored if not needed */
286
 
 
287
 
  if (h0 != 0.0 || hmax != 0.0 || hmin != 0.0 || ixpr != 0 || mxstep != 0 || mxhnil != 0 || mxordn != 0 || mxords != 0) {
288
 
  rwork[4] = h0; rwork[5] = hmax; rwork[6] = hmin;
289
 
  iwork[4] = ixpr; iwork[5] = mxstep; iwork[6] = mxhnil;
290
 
  iwork[7] = mxordn; iwork[8] = mxords;
291
 
  iopt = 1;
292
 
  }
293
 
  istate = 1;
294
 
  k = 1;
295
 
 
296
 
  /* If full output make some useful output arrays */
297
 
  if (full_output) {
298
 
    out_sz = ntimes-1;
299
 
    ap_hu = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_DOUBLE);
300
 
    ap_tcur = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_DOUBLE);
301
 
    ap_tolsf = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_DOUBLE);
302
 
    ap_tsw = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_DOUBLE);
303
 
    ap_nst = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_INT);
304
 
    ap_nfe = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_INT);
305
 
    ap_nje = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_INT);
306
 
    ap_nqu = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_INT);
307
 
    ap_mused = (PyArrayObject *)PyArray_FromDims(1,&out_sz,PyArray_INT);
308
 
    if (ap_hu == NULL || ap_tcur == NULL || ap_tolsf == NULL || ap_tsw == NULL || ap_nst == NULL || ap_nfe == NULL || ap_nje == NULL || ap_nqu == NULL || ap_mused == NULL) goto fail;
309
 
  }
310
 
 
311
 
  if (o_tcrit != NULL) {itask = 4; rwork[0] = *tcrit;}  /* There are critical points */
312
 
  while (k < ntimes && istate > 0) {    /* loop over desired times */
313
 
 
314
 
    tout_ptr = tout + k;
315
 
    /* Use tcrit if relevant */
316
 
    if (itask == 4 && *tout_ptr > *(tcrit + crit_ind)) {crit_ind++; rwork[0] = *(tcrit+crit_ind);}
317
 
    if (crit_ind >= numcrit) itask = 1;  /* No more critical values */
318
 
 
319
 
    LSODA(ode_function, &neq, y, &t, tout_ptr, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &lrw, iwork, &liw, ode_jacobian_function, &jt);
320
 
    if (full_output) {
321
 
      *((double *)ap_hu->data + (k-1)) = rwork[10];
322
 
      *((double *)ap_tcur->data + (k-1)) = rwork[12];
323
 
      *((double *)ap_tolsf->data + (k-1)) = rwork[13];
324
 
      *((double *)ap_tsw->data + (k-1)) = rwork[14];
325
 
      *((int *)ap_nst->data + (k-1)) = iwork[10];
326
 
      *((int *)ap_nfe->data + (k-1)) = iwork[11];
327
 
      *((int *)ap_nje->data + (k-1)) = iwork[12];
328
 
      *((int *)ap_nqu->data + (k-1)) = iwork[13];
329
 
      imxer = iwork[15];
330
 
      lenrw = iwork[16];
331
 
      leniw = iwork[17];
332
 
      *((int *)ap_mused->data + (k-1)) = iwork[18];
333
 
    }
334
 
    if (PyErr_Occurred()) goto fail;
335
 
    memcpy(yout_ptr, y, neq*sizeof(double));  /* copy integration result to output*/
336
 
    yout_ptr += neq;  k++;
337
 
  }
338
 
 
339
 
  RESTORE_JAC_FUNC();
340
 
 
341
 
  Py_DECREF(extra_args);
342
 
  Py_DECREF(ap_atol);
343
 
  Py_DECREF(ap_rtol);
344
 
  Py_XDECREF(ap_tcrit);
345
 
  Py_DECREF(ap_y);
346
 
  Py_DECREF(ap_tout);
347
 
  free(wa);
348
 
 
349
 
  /* Do Full output */
350
 
    if (full_output) {
351
 
      return Py_BuildValue("N{s:N,s:N,s:N,s:N,s:N,s:N,s:N,s:N,s:i,s:i,s:i,s:N}i",PyArray_Return(ap_yout),"hu",PyArray_Return(ap_hu),"tcur",PyArray_Return(ap_tcur),"tolsf",PyArray_Return(ap_tolsf),"tsw",PyArray_Return(ap_tsw),"nst",PyArray_Return(ap_nst),"nfe",PyArray_Return(ap_nfe),"nje",PyArray_Return(ap_nje),"nqu",PyArray_Return(ap_nqu),"imxer",imxer,"lenrw",lenrw,"leniw",leniw,"mused",PyArray_Return(ap_mused),istate);
352
 
    }
353
 
    else {
354
 
      return Py_BuildValue("Ni",PyArray_Return(ap_yout),istate);
355
 
    }
356
 
    
357
 
 fail:
358
 
  RESTORE_JAC_FUNC();
359
 
  Py_XDECREF(extra_args);
360
 
  Py_XDECREF(ap_y);
361
 
  Py_XDECREF(ap_rtol);
362
 
  Py_XDECREF(ap_atol);
363
 
  Py_XDECREF(ap_tcrit);
364
 
  Py_XDECREF(ap_tout);
365
 
  Py_XDECREF(ap_yout);
366
 
  if (allocated) free(wa);
367
 
  if (full_output) {
368
 
    Py_XDECREF(ap_hu);
369
 
    Py_XDECREF(ap_tcur);
370
 
    Py_XDECREF(ap_tolsf);
371
 
    Py_XDECREF(ap_tsw);
372
 
    Py_XDECREF(ap_nst);
373
 
    Py_XDECREF(ap_nfe);
374
 
    Py_XDECREF(ap_nje);
375
 
    Py_XDECREF(ap_nqu);
376
 
    Py_XDECREF(ap_mused);
377
 
  }
378
 
  return NULL;
379
 
  
380
 
}