Distributing Cython based extensions using LAPACK

Question

I am writing a Python module that includes Cython extensions and uses LAPACK (and BLAS). I am open to using either clapack or lapacke, or some kind of f2c or f2py solution if necessary. What is important is that I am able to call lapack and blas routines from Cython in tight loops without Python call overhead.

I've found one example here. However, that example depends on SAGE. I want my module to be installable without installing SAGE, since my users are not likely to want or need SAGE for anything else. My users are likely to have packages like numpy, scipy, pandas, and scikit learn installed, so those would be reasonable dependencies. What is the best combination of interfaces to use, and what would the minimal setup.py file look like that could fetch the necessary information (from numpy, scipy, etc.) for compilation?

EDIT: Here is what I ended up doing. It works on my macbook, but I have no idea how portable it is. Surely there's a better way.

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy
from Cython.Build import cythonize
from numpy.distutils.system_info import get_info

# TODO: This cannot be the right way
blas_include = get_info('blas_opt')['extra_compile_args'][1][2:]
includes = [blas_include,numpy.get_include()]

setup(
    cmdclass = {'build_ext': build_ext},
    ext_modules = cythonize([Extension("cylapack", ["cylapack.pyx"],
                                       include_dirs = includes,
                                       libraries=['blas','lapack'])
                   ])
)

This works because, on my macbook, the clapack.h header file is in the same directory as cblas.h. I can then do this in my pyx file:

ctypedef np.int32_t integer

cdef extern from "cblas.h":
    double cblas_dnrm2(int N,double *X, int incX)
cdef extern from "clapack.h":
    integer dgelsy_(integer *m, integer *n, integer *nrhs, 
    double *a, integer *lda, double *b, integer *ldb, integer *
    jpvt, double *rcond, integer *rank, double *work, integer *
    lwork, integer *info)

Answer 1

If I have understood the question correctly, you could make use of SciPy's Cython wrappers for BLAS and LAPACK routines. These wrappers are documented here:

• BLAS
• LAPACK

As the documentation states, you are responsible for checking that any arrays that you pass to these functions are aligned correctly for the Fortran routines. You can simply import and use these functions as needed in your .pyx file. For instance:

from scipy.linalg.cython_blas cimport dnrm2 
from scipy.linalg.cython_lapack cimport dgelsy 

Given that this is well-tested, widely-used code that runs on different platforms, I'd argue that it is a good candidate for reliably distributing Cython extensions that directly call BLAS and LAPACK routines.

---

If you do not want your code to have a dependency on the entirety of SciPy, you can find many of the relevant files for these wrapper functions in SciPy's linalg directory here. A useful reference is these lines of setup.py which list the source and header files. Note that a Fortran compiler is required!

In theory it should be possible to isolate only the source files here that are needed to compile the BLAS and LAPACK Cython wrappers and then bundle them as an independent extension with your module.

In practice this is very fiddly to do. The build process for the linalg submodule requires some Python functions to aid the compilation on different platforms (e.g. from here). Building also relies upon other C and Fortran source files (here), the paths of which are hard-coded into these Python functions.

Clearly a lot of work has gone into making sure that SciPy compiles properly on different operating systems and architectures.

I'm sure it is possible to do, but after shuffling files about and tweaking paths, I have not yet found the right way to build this part of the linalg submodule independently from the rest of SciPy. Should I find the correct way, I'll be sure to update this answer.