1 /* MULTIPACK module by Travis Oliphant
2
3 Copyright (c) 2002 Travis Oliphant all rights reserved
4 Oliphant.Travis@altavista.net
5 Permission to use, modify, and distribute this software is given under the
6 terms of the SciPy (BSD style) license. See LICENSE.txt that came with
7 this distribution for specifics.
8
9 NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
10 */
11
12
13 /* This extension module is a collection of wrapper functions around
14 common FORTRAN code in the packages MINPACK, ODEPACK, and QUADPACK plus
15 some differential algebraic equation solvers.
16
17 The wrappers are meant to be nearly direct translations between the
18 FORTRAN code and Python. Some parameters like sizes do not need to be
19 passed since they are available from the objects.
20
21 It is anticipated that a pure Python module be written to call these lower
22 level routines and make a simpler user interface. All of the routines define
23 default values for little-used parameters so that even the raw routines are
24 quite useful without a separate wrapper.
25
26 FORTRAN Outputs that are not either an error indicator or the sought-after
27 results are placed in a dictionary and returned as an optional member of
28 the result tuple when the full_output argument is non-zero.
29 */
30
31 #include "Python.h"
32 #include "numpy/arrayobject.h"
33 #include "ccallback.h"
34
35 #define PYERR(errobj,message) {PyErr_SetString(errobj,message); goto fail;}
36 #define PYERR2(errobj,message) {PyErr_Print(); PyErr_SetString(errobj, message); goto fail;}
37
38 #define STORE_VARS() ccallback_t callback; int callback_inited = 0; jac_callback_info_t jac_callback_info;
39
40 #define INIT_FUNC(fun,arg,errobj) do { /* Get extra arguments or set to zero length tuple */ \
41 if (arg == NULL) { \
42 if ((arg = PyTuple_New(0)) == NULL) goto fail_free; \
43 } \
44 else \
45 Py_INCREF(arg); /* We decrement on exit. */ \
46 if (!PyTuple_Check(arg)) \
47 PYERR(errobj,"Extra Arguments must be in a tuple"); \
48 /* Set up callback functions */ \
49 if (!PyCallable_Check(fun)) \
50 PYERR(errobj,"First argument must be a callable function."); \
51 if (init_callback(&callback, fun, arg) != 0) \
52 PYERR(errobj,"Could not init callback");\
53 callback_inited = 1; \
54 } while(0)
55
56 #define INIT_JAC_FUNC(fun,Dfun,arg,col_deriv,errobj) do { \
57 if (arg == NULL) { \
58 if ((arg = PyTuple_New(0)) == NULL) goto fail_free; \
59 } \
60 else \
61 Py_INCREF(arg); /* We decrement on exit. */ \
62 if (!PyTuple_Check(arg)) \
63 PYERR(errobj,"Extra Arguments must be in a tuple"); \
64 /* Set up callback functions */ \
65 if (!PyCallable_Check(fun) || (Dfun != Py_None && !PyCallable_Check(Dfun))) \
66 PYERR(errobj,"The function and its Jacobian must be callable functions."); \
67 if (init_jac_callback(&callback, &jac_callback_info, fun, Dfun, arg, col_deriv) != 0) \
68 PYERR(errobj,"Could not init callback");\
69 callback_inited = 1; \
70 } while(0)
71
72 #define RESTORE_JAC_FUNC() do { \
73 if (callback_inited && release_callback(&callback) != 0) { \
74 goto fail_free; \
75 }} while(0)
76
77 #define RESTORE_FUNC() do { \
78 if (callback_inited && release_callback(&callback) != 0) { \
79 goto fail_free; \
80 }} while(0)
81
82 #define SET_DIAG(ap_diag,o_diag,mode) { /* Set the diag vector from input */ \
83 if (o_diag == NULL || o_diag == Py_None) { \
84 ap_diag = (PyArrayObject *)PyArray_SimpleNew(1,&n,NPY_DOUBLE); \
85 if (ap_diag == NULL) goto fail; \
86 diag = (double *)PyArray_DATA(ap_diag); \
87 mode = 1; \
88 } \
89 else { \
90 ap_diag = (PyArrayObject *)PyArray_ContiguousFromObject(o_diag, NPY_DOUBLE, 1, 1); \
91 if (ap_diag == NULL) goto fail; \
92 diag = (double *)PyArray_DATA(ap_diag); \
93 mode = 2; } }
94
95 #define MATRIXC2F(jac,data,m,n) {double *p1=(double *)(jac), *p2, *p3=(double *)(data);\
96 int i,j;\
97 for (j=0;j<(m);p3++,j++) \
98 for (p2=p3,i=0;i<(n);p2+=(m),i++,p1++) \
99 *p1 = *p2; }
100
101 typedef struct {
102 PyObject *Dfun;
103 PyObject *extra_args;
104 int jac_transpose;
105 } jac_callback_info_t;
106
call_python_function(PyObject * func,npy_intp n,double * x,PyObject * args,int dim,PyObject * error_obj,npy_intp out_size)107 static PyObject *call_python_function(PyObject *func, npy_intp n, double *x, PyObject *args, int dim, PyObject *error_obj, npy_intp out_size)
108 {
109 /*
110 This is a generic function to call a python function that takes a 1-D
111 sequence as a first argument and optional extra_arguments (should be a
112 zero-length tuple if none desired). The result of the function is
113 returned in a multiarray object.
114 -- build sequence object from values in x.
115 -- add extra arguments (if any) to an argument list.
116 -- call Python callable object
117 -- check if error occurred:
118 if so return NULL
119 -- if no error, place result of Python code into multiarray object.
120 */
121
122 PyArrayObject *sequence = NULL;
123 PyObject *arglist = NULL;
124 PyObject *arg1 = NULL;
125 PyObject *result = NULL;
126 PyArrayObject *result_array = NULL;
127 npy_intp fvec_sz = 0;
128
129 /* Build sequence argument from inputs */
130 sequence = (PyArrayObject *)PyArray_SimpleNewFromData(1, &n, NPY_DOUBLE, (char *)x);
131 if (sequence == NULL) PYERR2(error_obj,"Internal failure to make an array of doubles out of first\n argument to function call.");
132
133 /* Build argument list */
134 if ((arg1 = PyTuple_New(1)) == NULL) {
135 Py_DECREF(sequence);
136 return NULL;
137 }
138 PyTuple_SET_ITEM(arg1, 0, (PyObject *)sequence);
139 /* arg1 now owns sequence reference */
140 if ((arglist = PySequence_Concat( arg1, args)) == NULL)
141 PYERR2(error_obj,"Internal error constructing argument list.");
142
143 Py_DECREF(arg1); /* arglist has a reference to sequence, now. */
144 arg1 = NULL;
145
146 /* Call function object --- variable passed to routine. Extra
147 arguments are in another passed variable.
148 */
149 if ((result = PyEval_CallObject(func, arglist))==NULL) {
150 goto fail;
151 }
152
153 if ((result_array = (PyArrayObject *)PyArray_ContiguousFromObject(result, NPY_DOUBLE, dim-1, dim))==NULL)
154 PYERR2(error_obj,"Result from function call is not a proper array of floats.");
155
156 fvec_sz = PyArray_SIZE(result_array);
157 if(out_size != -1 && fvec_sz != out_size){
158 PyErr_SetString(PyExc_ValueError, "The array returned by a function changed size between calls");
159 Py_DECREF(result_array);
160 goto fail;
161 }
162
163 Py_DECREF(result);
164 Py_DECREF(arglist);
165 return (PyObject *)result_array;
166
167 fail:
168 Py_XDECREF(arglist);
169 Py_XDECREF(result);
170 Py_XDECREF(arg1);
171 return NULL;
172 }
173