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.

scipy allready has blas-functions wrapped, you could just use it https://docs.scipy.org/doc/scipy-0.19.0/reference/linalg.cython_blas.html#module-scipy.linalg.cython_blas — ead, Oct 21 '18 at 05:17
But that's for using them in Cython (in the .pyx file). Can that somehow be passed to the C functions in the .c file? — anymous.asker, Oct 21 '18 at 05:20

ead · Accepted Answer · 2018-10-22T19:08:19.063

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:

Define an explicit function-pointer to ddot-functionality, called my_ddot in the snippet below.
Use my_ddot-pointer where you would use ddot-otherwise.
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.

anymous.asker · Answer 2 · 2021-01-07T13:28:23.810

I've finally figured out an ugly hack for this. I'm not sure if it will always work, but at least it works for cobminations of Windows (mingw and visual studio), Linux, MKL and OpenBlas. I'd still like to know if there are better alternatives, but if not, this will do it:

Edit: Corrected for visual studio now

Modify C files to account for names with underscores (do it for each BLAS function that is called) - need to declare each function twice and add an if for each one

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

daxpy_(int *N, double *DA, double *DX, int *INCX, double *DY, int *INCY); #define daxpy(N, DA, DX, INCX, DY, INCY) daxpy_(N, DA, DX, INCX, DY, INCY)

... etc
Extract library path from NumPy or SciPy and add it to the link arguments.
Detect if the compiler to be used is visual studio, in which case the linking arguments are quite different.

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
from sys import platform
import os

try:
    blas_path = numpy.distutils.system_info.get_info('blas')['library_dirs'][0]
except:
    if "library_dirs" in numpy.__config__.blas_mkl_info:
        blas_path = numpy.__config__.blas_mkl_info["library_dirs"][0]
    elif "library_dirs" in numpy.__config__.blas_opt_info:
        blas_path = numpy.__config__.blas_opt_info["library_dirs"][0]
    else:
        raise ValueError("Could not locate BLAS library.")
        

if platform[:3] == "win":
    if os.path.exists(os.path.join(blas_path, "mkl_rt.lib")):
        blas_file = "mkl_rt.lib"
    elif os.path.exists(os.path.join(blas_path, "mkl_rt.dll")):
        blas_file = "mkl_rt.dll"
    else:
        import re
        blas_file = [f for f in os.listdir(blas_path) if bool(re.search("blas", f))]
        if len(blas_file) == 0:
            raise ValueError("Could not locate BLAS library.")
        blas_file = blas_file[0]
        
elif platform[:3] == "dar":
    blas_file = "libblas.dylib"
else:
    blas_file = "libblas.so"

## https://stackoverflow.com/questions/724664/python-distutils-how-to-get-a-compiler-that-is-going-to-be-used
class build_ext_subclass( build_ext ):
    def build_extensions(self):
        compiler = self.compiler.compiler_type
        if compiler == 'msvc': # visual studio
            for e in self.extensions:
                e.extra_link_args += [os.path.join(blas_path, blas_file)]
        else: # gcc
            for e in self.extensions:
                e.extra_link_args += ["-L"+blas_path, "-l:"+blas_file]
        build_ext.build_extensions(self)


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

I think if you use `libraries`-options there is no need to construct the right name - distutils will do it for you: https://docs.python.org/2/distutils/setupscript.html#library-options — ead, Oct 22 '18 at 09:14

score 0 · Answer 3 · answered Sep 06 '21 at 17:11

As yet another alternative with more recent Cython versions, one can create a "public" Cython function (which will be made available to C code and auto-generate a public header) that would simply call the corresponding BLAS function:

from scipy.linalg.cython_blas cimport ddot
cdef public double ddot_(int *n, double *x, int *ldx, double *y, int *ldy):
    return ddot(n, x, ldx, y, ldy)

Then one simply declares it in the C code or includes the header, and the rest of the Cython extension builder will take care of linkage:

extern double ddot_(int *n, double *x, int *ldx, double *y, int *ldy);

Link Cython-wrapped C functions against BLAS from NumPy

3 Answers3

Linked