Link Cython-wrapped C functions against BLAS from NumPy

Question

I want to use inside a Cython extension some C functions defined in .c files that uses BLAS subroutines, e.g.

cfile.c

double ddot(int *N, double *DX, int *INCX, double *DY, int *INCY);

double call_ddot(double* a, double* b, int n){
    int one = 1;
    return ddot(&n, a, &one, b, &one);
}

(Let’s say the functions do more than just call one BLAS subroutine)

pyfile.pyx

cimport numpy as np
import numpy as np

cdef extern from "cfile.c":
    double call_ddot(double* a, double* b, int n)

def pyfun(np.ndarray[double, ndim=1] a):
    return call_ddot(&a[0], &a[0], <int> a.shape[0])

setup.py:

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

setup(
    name  = "wrapped_cfun",
    packages = ["wrapped_cfun"],
    cmdclass = {'build_ext': build_ext},
    ext_modules = [Extension("wrapped_cfun.cython_part", sources=["pyfile.pyx"], include_dirs=[numpy.get_include()])]
)

I want this package to link against the same BLAS library that the installed NumPy or SciPy are using, and would like it to be installable from PIP under different operating systems using numpy or scipy as dependencies, without any additional BLAS-related dependency.

Is there any hack for setup.py that would allow me to accomplish this, in a way that it could work with any BLAS implementation?

Update: With MKL, I can make it work by modifying the Extension object to point to libmkl_rt, which can be extracted from numpy if MKL is installed, e.g.: Extension("wrapped_cfun.cython_part", sources=["pyfile.pyx"], include_dirs=[numpy.get_include()], extra_link_args=["-L{path to python's lib dir}", "-l:libmkl_rt.{so, dll, dylib}"]) However, the same trick does not work for OpenBLAS (e.g. -l:libopenblasp-r0.2.20.so). Pointing to libblas.{so,dll,dylib} will not work if that file is a link to libopenblas, but works fine it it's a link to libmkl_rt.

Update 2: It seems OpenBLAS names their C functions with an underscore at the end, e.g. not ddot but ddot_. The code above with l:libopenblas will work if I change ddot to ddot_ in the .c file. I'm still wondering if there is some (ideally run-time) mechanism to detect which name should be used in the c file.

Answer 1

An alternative to depending on linker/loader to provide the right blas-functionality, would be to emulate resolution of the necessary blas-symbols (e.g. ddot) and to use the wrapped blas-function provided by scipy during the runtime.

Not sure, this approach is superior to the "normal way" of building, but wanted to bring it to your attention, even if only because I find this approach interesting.

The idea in a nutshell:

1. Define an explicit function-pointer to ddot-functionality, called my_ddot in the snippet below.
2. Use my_ddot-pointer where you would use ddot-otherwise.
3. Initialize my_ddot-pointer when the cython-module is loaded with the functionality provided by scipy.

Here is a working prototype (I use C-code-verbatim to make the snippet standalone and easily testable in a jupiter-notebook, trust you to transform it to format you need/like):

%%cython
# h-file:
cdef extern from *:
    """
    // blas-functionality,
    // will be initialized by cython when module is loaded:
    typedef double (*ddot_t)(int *N, double *DX, int *INCX, double *DY, int *INCY);
    extern ddot_t my_ddot;

    double call_ddot(double* a, double* b, int n);
    """
    ctypedef double (*ddot_t)(int *N, double *DX, int *INCX, double *DY, int *INCY)
    ddot_t my_ddot
    double call_ddot(double* a, double* b, int n)    

# init the functions of the c-library
# with blas-function provided by scipy
from scipy.linalg.cython_blas cimport ddot
my_ddot=ddot

# a simple function to demonstrate, that it works
def ddot_mult(double[:]a, double[:]b):
    cdef int n=len(a)
    return call_ddot(&a[0], &b[0], n)

#-------------------------------------------------
# c-file, added so the example is complete    
cdef extern from *:
    """  
    ddot_t my_ddot;
    double call_ddot(double* a, double* b, int n){
        int one = 1;
        return my_ddot(&n, a, &one, b, &one);
    }
    """
    pass

And now ddot_mult can be used:

import numpy as np
a=np.arange(4, dtype=float)

ddot_mult(a,a)  # 14.0 as expected!

An advantage of this approach is, that there is no hustle with distutils and you have a guarantee, to use the same blas-functionality as scipy.

Another perk: One could switch the used engine (mkl, open_blas or even an own implementation) during the runtime without the need to recompile/relink.

On there other hand, there is some additional amount of boilerplate-code and also the danger, that initialization of some symbols will be forgotten.