[SciPy-dev] Extensions linking to libraries used by numpy/scipy

[Robert Kern]

Tom Loredo wrote:
> Hi folks,
> 
> I asked about this some time ago and was told there was no 
> straightforward way to accomplish it, but that was over a
> year ago and before the recent numpy/scipy split, etc..  So
> I thought I'd ask again.
> 
> I'm writing a package with several C and Fortran extensions,
> some of which call functions or subroutines that do things
> already done by Scipy routines, e.g., evaluate a gamma function,
> perform an LU decomposition, etc..  There are repeated calls,
> in the course of a more complex calculation, and the most
> sensible way to do the calculations would be to call Scipy's library
> routines directly, and not go back through Python.  In some cases as 
> a placeholder I have a version derived from Numerical Recipes, and 
> this is a no-no as far as public distribution is concerned.  I'd like
> to be able to call & link to the C or Fortran libraries that Scipy
> is accessing, but I need to do it in a way that will build
> and operate reliably across the various platforms Scipy builds on.
> 
> As specific examples, I'd like to access functions in
> cephes/gamma.c, and lapack routines for LU decomposition
> and determinant calcuation, e.g., DGETRF, DGBTRF or DGTTRF.
> 
> I suspect this isn't hard for gamma.c, but since Scipy 
> is sometimes built with BLAS/LAPACK, and sometimes with ATLAS,
> I suspect it's harder for the LU decomposition.
> 
> So... is there a way to do this portably?  If so, is it
> documented somewhere (example code and an example setup.py)?  
> Are there examples I could look at?  Do some particular
> Scipy modules provide good examples of this (which ones)?

I don't think anything has changed since the last time you asked. Sorry.

To do this reliably across platforms, you have to do the equivalent of numpy's
import_array() trick. The "exporting" extension has to be written specifically
to do it, and nothing in scipy has been. Possibly, some should. It occurs to me
that it might be possible to modify f2py to do add the necessary function tables
to the extension module. That would be a neat trick.

Also, I remember asking Pearu to expose the raw function pointer in the fortran
object wrapper. The idea was to be able to allow f2py callback functions to be
f2py'ed subroutines themselves, and have everything go really, really fast.

In [6]: linalg.flapack.dgetrf._cpointer
Out[6]: <PyCObject object at 0xd650>

I believe this actually gives you the function pointer to DGETRF_. Or you could
just accept the fortran object itself, and access the function pointer from the
PyFortranObject structure itself.

Ah, now I've rambled myself to the right answer. Don't pay attention to anything
but the last sentence of the previous paragraph. AFAICT, this is a new approach,
so you can be our guinea pig.

For gamma, I recommend just yoinking the routine from specfun. It's tiny.

-- 
Robert Kern
robert.kern at gmail.com

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco