diff --git a/pyradau13.c b/pyradau13.c
index 5b24e21768f1909227ae5e7676239f0948029e9a..de80c002bc22ee30bb14721e6287b20de12f6699 100644
--- a/pyradau13.c
+++ b/pyradau13.c
@@ -8,7 +8,7 @@
typedef void (*rhs_t)(int *n, double *x, double *y, double *f, float *rpar, int *ipar);
typedef void (*jac_t)(int *n, double *x, double *y, double *dfy, int *ldfy, float *rpar, int *ipar);
typedef void (*mas_t)(int *n, double *am, int *lmas, float *rpar, int *ipar);
-typedef void (*sol_t)(int *nr, double *xold, double *x, double *y, int *cont, int *lrc, int *n,
+typedef void (*sol_t)(int *nr, double *xold, double *x, double *y, double *cont, int *lrc, int *n,
float *rpar, int *ipar, int *irtrn);
#ifdef __cplusplus
@@ -49,6 +49,13 @@ extern "C" {
// IDID=-3 STEP SIZE BECOMES TOO SMALL,
// IDID=-4 MATRIX IS REPEATEDLY SINGULAR.
);
+
+ double contra_(
+ int *I, // Dimension
+ double *X, // Time
+ double *CONT, // Continous data
+ int *LRC // Dimensionality of CONT
+ );
#ifdef __cplusplus
}
#endif
@@ -59,6 +66,85 @@ struct radau_options {
PyArrayObject *y_out;
};
+typedef struct {
+ PyObject_HEAD
+
+ PyArrayObject *y;
+ double t0;
+ double t1;
+ double *cont;
+ int lrc;
+} RadauDenseEvaluator;
+
+static PyObject *RadauDenseEvaluator_call(RadauDenseEvaluator *self, PyObject *args, PyObject *kwargs) {
+ double time;
+ if(self->cont == NULL) {
+ PyErr_SetString(PyExc_RuntimeError, "This object is only valid in the dense callback function");
+ return NULL;
+ }
+ if(!PyArg_ParseTuple(args, "d", &time)) return NULL;
+ if(time < self->t0 || time > self->t1) {
+ char err[255];
+ sprintf(err, "Time must fulfil %e <= t <= %e", self->t0, self->t1);
+ PyErr_SetString(PyExc_ValueError, err);
+ return NULL;
+ }
+
+ PyArrayObject *retval = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(self->y), PyArray_DIMS(self->y), NPY_DOUBLE);
+ double *data = PyArray_DATA(retval);
+
+ int n;
+ for(n = 1; n<PyArray_SIZE(self->y)+1; n++) {
+ data[n-1] = contra_(&n, &time, self->cont, &(self->lrc));
+ }
+
+ return (PyObject *)retval;
+}
+
+static PyTypeObject RadauDenseEvaluatorType = {
+ PyObject_HEAD_INIT(NULL)
+ 0, /* ob_size*/
+ "RadauDenseEvaluator", /* tp_name*/
+ sizeof(RadauDenseEvaluator), /* tp_basicsize*/
+ 0, /* tp_itemsize*/
+ 0, /* tp_dealloc*/
+ 0, /* tp_print*/
+ 0, /* tp_getattr*/
+ 0, /* tp_setattr*/
+ 0, /* tp_compare*/
+ 0, /* tp_repr*/
+ 0, /* tp_as_number*/
+ 0, /* tp_as_sequence*/
+ 0, /* tp_as_mapping*/
+ 0, /* tp_hash */
+ (ternaryfunc)RadauDenseEvaluator_call, /* tp_call*/
+ 0, /* tp_str*/
+ 0, /* tp_getattro*/
+ 0, /* tp_setattro*/
+ 0, /* tp_as_buffer*/
+ Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /* tp_flags*/
+ "Dense evaluator", /* tp_doc */
+ 0, /* tp_traverse */
+ 0, /* tp_clear */
+ 0, /* tp_richcompare */
+ 0, /* tp_weaklistoffset */
+ 0, /* tp_iter */
+ 0, /* tp_iternext */
+ 0, /* tp_methods */
+ 0, /* tp_members */
+ 0, /* tp_getset */
+ 0, /* tp_base */
+ 0, /* tp_dict */
+ 0, /* tp_descr_get */
+ 0, /* tp_descr_set */
+ 0, /* tp_dictoffset */
+ 0, /* tp_init */
+ 0, /* tp_alloc */
+ 0, /* tp_new */
+};
+
+
+
static void radau_rhs(int *n, double *x, double *y, double *f, float *rpar, int *ipar) {
struct radau_options *options = (struct radau_options *)ipar;
@@ -76,24 +162,34 @@ static void radau_rhs(int *n, double *x, double *y, double *f, float *rpar, int
Py_DECREF(rhs_retval);
}
-static void radau_dense_feedback(int *nr, double *xold, double *x, double *y, int *cont, int *lrc, int *n, float *rpar, int *ipar, int *irtrn) {
+static void radau_dense_feedback(int *nr, double *xold, double *x, double *y, double *cont, int *lrc, int *n, float *rpar, int *ipar, int *irtrn) {
struct radau_options *options = (struct radau_options *)ipar;
PyObject *y_current = PyArray_SimpleNewFromData(PyArray_NDIM(options->y_out), PyArray_DIMS(options->y_out), NPY_DOUBLE, y);
if(y_current == NULL) return;
PyObject *y_copy = (PyObject *)PyArray_Copy((PyArrayObject *)y_current);
if(y_copy == NULL) return;
- PyObject *rhs_retval = PyObject_CallFunction(options->dense_callback, "dO", *x, y_copy);
+
+ RadauDenseEvaluator *evaluator = PyObject_NEW(RadauDenseEvaluator, &RadauDenseEvaluatorType);
+ evaluator->t1 = *x;
+ evaluator->t0 = *xold;
+ evaluator->lrc = *lrc;
+ evaluator->cont = cont;
+ evaluator->y = (PyArrayObject *)y_copy;
+
+ PyObject *rhs_retval = PyObject_CallFunction(options->dense_callback, "ddOO", *xold, *x, y_copy, evaluator);
Py_DECREF(y_copy);
Py_DECREF(y_current);
+ evaluator->cont = NULL;
+ Py_DECREF(evaluator);
- if(PyObject_IsTrue(rhs_retval)) {
+ if(rhs_retval == NULL || PyObject_IsTrue(rhs_retval)) {
*irtrn = -1;
}
else {
*irtrn = 1;
}
- Py_DECREF(rhs_retval);
+ if(rhs_retval) Py_DECREF(rhs_retval);
}
static PyObject *radau13(PyObject *self, PyObject *args, PyObject *kwargs) {
@@ -211,9 +307,11 @@ static PyObject *radau13(PyObject *self, PyObject *args, PyObject *kwargs) {
Py_DECREF(y0_array);
if(idid != 1 || y_out == NULL) {
- char err[30];
- sprintf(err, "radau failed with idid = %d", idid);
- PyErr_SetString(PyExc_RuntimeError, err);
+ if(PyErr_Occurred() == NULL) {
+ char err[30];
+ sprintf(err, "radau failed with idid = %d", idid);
+ PyErr_SetString(PyExc_RuntimeError, err);
+ }
if(time_levels > 1) {
Py_DECREF(list_retval);
}
@@ -250,8 +348,9 @@ static PyMethodDef methods[] = {
" - order must be one of 13, 9 or 5.\n"
" - max_steps denotes the maximal step count to take.\n"
" - dense_callback must be either not set, or a callable of the\n"
- " type dense_callback(t, y). It is called after each internal\n"
- " integration step.\n"
+ " type dense_callback(told, t, y, cont). It is called after each
+ " internal integration step. cont is a callable that takes an argument\n"
+ " from [told, t] and returns an approximation to the value at that time.\n"
" If dense_callback returns a value that evaluates to True, integration\n"
" is halted.\n"
" - t0 denotes the initial time.\n"
@@ -274,6 +373,10 @@ static PyMethodDef methods[] = {
};
PyMODINIT_FUNC initpyradau13(void) {
- (void)Py_InitModule("pyradau13", methods);
- import_array();
+ RadauDenseEvaluatorType.tp_new = PyType_GenericNew;
+ if(PyType_Ready(&RadauDenseEvaluatorType) < 0)
+ return;
+
+ (void)Py_InitModule("pyradau13", methods);
+ import_array();
}
diff --git a/setup.py b/setup.py
index 48689122cdc462f3e6834d806808fcd889823d93..cb47bf08faa868a60bff8057d7a74c5ec8f99487 100755
--- a/setup.py
+++ b/setup.py
@@ -16,7 +16,7 @@ Hairer, Norsett and Wanner (1993): Solving Ordinary Differential Equations.
Nonstiff Problems. 2nd edition. Springer Series in Comput. Math., vol. 8.
This package does not use f2py, but instead calls the fortran function from a C
-extension. Dense output is available. Supplying a Jacobian function is
+extension. Dense/continous output is available. Supplying a Jacobian function is
not (yet).
"""
diff --git a/test.py b/test.py
index aab3e27a6c7b504fad6353e2f642bc90e5f48ed2..5380b5155797eed325cf41d0b53722e773409454 100755
--- a/test.py
+++ b/test.py
@@ -22,7 +22,7 @@ class TestIntegration(unittest.TestCase):
def test_dense_feedback(self):
_x = []
_t = []
- def _dense_cb(t, x):
+ def _dense_cb(told, t, x, cont):
_x.append(x[0])
_t.append(t)
radau13(self._pendulum_rhs, [1, 0], 10, dense_callback=_dense_cb)
@@ -36,6 +36,12 @@ class TestIntegration(unittest.TestCase):
y = radau13(lambda t, x: x, 1, X)
self.assertAlmostEqual(max(abs(np.exp(X) - y)), 0)
+ def test_continous_output(self):
+ def _dense_cb(told, t, x, cont):
+ tavg = (told + t) / 2
+ self.assertAlmostEqual(cont(tavg), np.exp(tavg))
+ radau13(lambda t, x: x, 1, 1, dense_callback=_dense_cb)
+
if __name__ == '__main__':
unittest.main()