C extension using GSL

jesse jberwald at gmail.com
Fri Mar 27 07:10:06 CET 2009

I give up. I cannot find my memory leak! I'm hoping that someone out
there has come across something similar. Let me lay out the basic

I'm performing multiple simulations on a model. Each iteration
involves solving a system of differential equations. For this I use
the GNU Scientific Library (GSL) -- specifically the rk4imp ODE
solver. After the ODE is solved the array is returned to python and is
analyzed. As you may have guessed, my problem is that over the course
of the iterations the memory keeps climbing until python crashes.

Note: *The extension does not keep running.* It returns object (a
list) and is done until the next iteration. AFAIK, any memory
allocated during execution of the extension should be released.
Question: Since the extension was run from within python is memory
allocated within an extension part of python's heap?  Would this have
an adverse or unpredictable affect on any memory allocated during the
running of the extension? One hypothesis I have, since most other
possibilities have been eliminated, is that GSL's creation of it's own
data structures (eg., gsl_vector) is messing up Python's control of
the heap. Is this possible? If so, how would one be able to fix this

It may help some nice folks out there who are good enough to look at
this post if I layout the flow, broken up by where stuff happens
(i.e., in the Python code or C code):

1) Python: Set up the simulation and basic data structures (numpy
arrays with initial conditions, coupling matrices for the ODE's, dicts
with parameters, etc).
2) Python: Pass these to the C extension
3) C: Python objects passed in are converted to C arrays, floats,
4) C: A PyList object, L,  is created (new reference!). This will hold
the solution vector for the ODE
5) C: Initialize GSL ODE solver and stepper functions. Solve the ODE,
at each step use PyList_Append(L, current_state) to append the current
state to L.
6) C: After the ODE solver finishes, free GSL objects, free coupling
matrices, free last state vector.
7) C: Return L to Python with return Py_BuildValue("N", L).
8) Python: Convert returned list L to array A, delete L, work with A.
8.1) Python: Step 8) includes plotting. (I have a small suspicion that
matplotlib holds onto lots of data, but I've used clf() and close(fig)
on all plots, so I think I'm safe here. )
8.2) Python: save analysis results from A, save A. (At this point
there should be no more use of A. In fact, at point 8) in the next
iteration A is replaced by a new array.)
9) Python: Change any parameters or initial conditions and goto 1).

thanks for any help,

More information about the Python-list mailing list