(2012) Accessing LAPACK and BLAS from the numpy C API
Hello, In 2009 there was a thread in this mailing list concerning the access to BLAS from C extension modules. If I have properly understood the thread: http://mail.scipy.org/pipermail/numpy-discussion/2009-November/046567.html the answer by then was that those functions were not exposed (only f2py functions). I just wanted to know if the situation has changed since 2009 because it is not uncommon that to optimize some operations one has to sooner or later access BLAS functions that are already wrapped in numpy (either from ATLAS, from the Intel MKL, ...) Thanks for your time, Armando
On 05.03.2012 14:26, "V. Armando Solé" wrote:
In 2009 there was a thread in this mailing list concerning the access to BLAS from C extension modules.
If I have properly understood the thread:
http://mail.scipy.org/pipermail/numpy-discussion/2009-November/046567.html
the answer by then was that those functions were not exposed (only f2py functions).
I just wanted to know if the situation has changed since 2009 because it is not uncommon that to optimize some operations one has to sooner or later access BLAS functions that are already wrapped in numpy (either from ATLAS, from the Intel MKL, ...)
Why do you want to do this? It does not make your life easier to use NumPy or SciPy's Python wrappers from C. Just use BLAS directly from C instead. BLAS is a Fortran 77 library (although it might be implemented in C or assembly). Fortran 77 compilers do not use a predefined binary interface. f2py has knowledge about all the major compilers, and will generate different call statements depending on compiler. NumPy and SciPy does not use the standard C BLAS interface, nor is all of BLAS and LAPACK exposed. You can, however, get a C function pointer to the part of BLAS and LAPACK that SciPy does expose: import scipy as sp from scipy.linalg import get_blas_funcs DGEMM = get_blas_funcs('gemm', dtype=np.float64) Now DGEMM._cpointer is a PyCObject that wraps the C function pointer. You can extract it by calling PyCObject_AsVoidPtr in C and cast the return value to the correct function pointer type. But be aware that you must know the binary interface of the underlying Fortran version. Generally you can assume that all arguments to a Fortran 77 subroutine are pointers (character strings can be problematic). For MKL you can let f2py generate code for the Intel Fortran compiler and use the same call statement in C. But then comes the question of legality: If you don't have a developer's license for MKL you are probably not allowed to use it like that. And if you do, you can just use the header files and the C BLAS interface instead. Generally, I will recommend that you build GotoBLAS2 (or OpenBLAS) if you don't have a license for MKL -- or download ACML from AMD. You will in any case get a set of C headers you can use from C. In either case it is just extra work to hack into SciPy's BLAS functions using the _cpointer attribute, nor do you gain anything from doing it. Sturla
On Tue, Mar 6, 2012 at 2:57 PM, Sturla Molden <sturla@molden.no> wrote:
On 05.03.2012 14:26, "V. Armando Solé" wrote:
In 2009 there was a thread in this mailing list concerning the access to BLAS from C extension modules.
If I have properly understood the thread:
http://mail.scipy.org/pipermail/numpy-discussion/2009-November/046567.html
the answer by then was that those functions were not exposed (only f2py functions).
I just wanted to know if the situation has changed since 2009 because it is not uncommon that to optimize some operations one has to sooner or later access BLAS functions that are already wrapped in numpy (either from ATLAS, from the Intel MKL, ...)
Why do you want to do this? It does not make your life easier to use NumPy or SciPy's Python wrappers from C. Just use BLAS directly from C instead.
Of course it does make his life easier. This way he does not have to distribute his own BLAS/LAPACK/etc... Please stop presenting as truth things which are at best highly opiniated. You already made such statements many times, and it is not helpful at all. David
On 06.03.2012 21:33, David Cournapeau wrote:
Of course it does make his life easier. This way he does not have to distribute his own BLAS/LAPACK/etc...
Please stop presenting as truth things which are at best highly opiniated. You already made such statements many times, and it is not helpful at all.
I showed him how to grab those function pointers if he wants to. Sturla
On 06/03/2012 20:57, Sturla Molden wrote:
On 05.03.2012 14:26, "V. Armando Solé" wrote:
In 2009 there was a thread in this mailing list concerning the access to BLAS from C extension modules.
If I have properly understood the thread:
http://mail.scipy.org/pipermail/numpy-discussion/2009-November/046567.html
the answer by then was that those functions were not exposed (only f2py functions).
I just wanted to know if the situation has changed since 2009 because it is not uncommon that to optimize some operations one has to sooner or later access BLAS functions that are already wrapped in numpy (either from ATLAS, from the Intel MKL, ...) Why do you want to do this? It does not make your life easier to use NumPy or SciPy's Python wrappers from C. Just use BLAS directly from C instead.
Wow! It certainly makes my life much, much easier. I can compile and distribute my python extension *even without having ATLAS, BLAS or MKL installed*. Please note I am not using the python wrappers from C. That would make no sense. I am using the underlying libraries supplied with python from C. I had already used the information Robert Kern provided on the 2009 thread and obtained the PyCObject as: from scipy.linalg.blas import fblas dgemm = fblas.dgemm._cpointer sgemm = fblas.sgemm._cpointer but I did not find a way to obtain those pointers from numpy. That was the goal of my post. My extension needs SciPy installed just to fetch the pointer. It would be very nice to have a way to get similar information from numpy. I have made a test on a Debian machine with BLAS installed but no ATLAS-> Extension slow but working. Then the system maintainer has installed ATLAS -> The extension flies. So, one can distribute a python extension that works on its own but that can take profit of any advanced library the end user might have installed. Your point of view is valid if one is not going to distribute the extension module but I *have to* distribute the module for Linux and for windows. To have a proper fortran compiler for windows 64 bit compatible with python is already an issue. If I have to distribute my own ATLAS or MKL then it gets even worse. All those issues are solved just by using the pointer to the function. Concerning licenses, if the end user has the right to use MKL, then he has the right to use it via my extension. It is not me who is using MKL Armando PS. The only issue I see with the whole approach is safety because the extension might be used to call some nasty function.
Den 07.03.2012 21:02, skrev "V. Armando Solé":
I had already used the information Robert Kern provided on the 2009 thread and obtained the PyCObject as:
from scipy.linalg.blas import fblas dgemm = fblas.dgemm._cpointer sgemm = fblas.sgemm._cpointer
but I did not find a way to obtain those pointers from numpy. That was the goal of my post. My extension needs SciPy installed just to fetch the pointer. It would be very nice to have a way to get similar information from numpy.
By the way, here is a code I wrote to use DGELS instead of DGELSS for least-squares (i.e. QR instead of SVD). It shows how we can grab a LAPACK function from MKL when using EPD. Sturla import numpy as np import scipy as sp from scipy.linalg import LinAlgError try: import ctypes from numpy.ctypeslib import ndpointer _c_int_p = ctypes.POINTER(ctypes.c_int) _c_double_p = ctypes.POINTER(ctypes.c_double) _double_array_1d = ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS' ) _int_array_1d = ndpointer(dtype=np.int32, ndim=1, flags='C_CONTIGUOUS' ) _double_array_2d = ndpointer(dtype=np.float64, ndim=2, flags='F_CONTIGUOUS' ) intel_mkl = ctypes.CDLL('mk2_rt.dll') dgels = intel_mkl.DGELS dgels.restype = None dgels.argtypes = (ctypes.c_char_p, _c_int_p, _c_int_p, _c_int_p, # TRANS, M, N, NRHS _double_array_2d, _c_int_p, # A, LDA _double_array_2d, _c_int_p, # b, LDB _double_array_1d, _c_int_p, _c_int_p) # WORK, LWORK, INFO _one = ctypes.c_int(1) _minus_one = ctypes.c_int(-1) _no_transpose = ctypes.c_char_p("N") def copy_fortran(x, dtype=np.float64): return np.array(x, dtype=dtype, copy=True, order='F') def lstsq( X, y ): assert(X.ndim == 2) assert(y.ndim == 1) X = copy_fortran(X) y = copy_fortran(y[:,np.newaxis]) assert(X.shape[0] == y.shape[0]) m = ctypes.c_int(X.shape[0]) n = ctypes.c_int(X.shape[1]) ldx = ctypes.c_int(m.value) ldy = ctypes.c_int(m.value) info = ctypes.c_int(0) swork = np.zeros(1, dtype=np.float64) dgels(_no_transpose, ctypes.byref(m), ctypes.byref(n), ctypes.byref(_one), X, ctypes.byref(ldx), y, ctypes.byref(ldy), swork, ctypes.byref(_minus_one), ctypes.byref(info)) if info.value < 0: raise ValueError, 'illegal argument to lapack dgels: arg no. %d' % (-info.value,) if info.value > 0: raise LinAlgError, 'lapack error %d, X does not have full rank' % (info.value,) lwork = ctypes.c_int(int(swork[0])) work = np.zeros(lwork.value, dtype=np.float64) dgels(_no_transpose, ctypes.byref(m), ctypes.byref(n), ctypes.byref(_one), X, ctypes.byref(ldx), y, ctypes.byref(ldy), work, ctypes.byref(lwork), ctypes.byref(info)) if info.value < 0: raise ValueError, 'illegal argument to lapack dgels: arg no. %d' % (-info.value,) if info.value != 0: raise LinAlgError, 'lapack error %d, X does not have full rank' % (info.value,) return (y[:X.shape[1],0],) except: from scipy.linalg import lstsq def linreg(y, X): return lstsq(X,y)[0]
Den 07.03.2012 21:02, skrev "V. Armando Solé":
I had already used the information Robert Kern provided on the 2009 thread and obtained the PyCObject as:
from scipy.linalg.blas import fblas dgemm = fblas.dgemm._cpointer sgemm = fblas.sgemm._cpointer
but I did not find a way to obtain those pointers from numpy. That was the goal of my post. My extension needs SciPy installed just to fetch the pointer. It would be very nice to have a way to get similar information from numpy.
The problem here is that NumPy's lapack_lite is compiled to C with f2c, and there is a handwritten C extension module (lapack_litemodule.c) to use it. Unlike SciPy, it does not use f2py to call Fortran LAPACK directly. I am not sure why NumPy uses f2c'd routines instead of a dependency on BLAS and LAPACK like SciPy. And by Murphy's law, function pointers to the f2c'd LAPACK routines are not exported from lapack_lite.pyd. So it does not even help to load it with ctypes as an ordinary DLL. Sturla
Den 10.03.2012 22:56, skrev Sturla Molden:
I am not sure why NumPy uses f2c'd routines instead of a dependency on BLAS and LAPACK like SciPy.
Actually, np.dot does depend on the CBLAS interface to BLAS (_dotblas.c). But the lapack methods in lapack_lite seems to use f2c'd code. I am not sure if they will use an optimized BLAS or just link to f2c's BLAS in blas_lite.c. If the intention is to avoid Fortran dependency in NumPy, I am not sure why this is better than a dependency on CBLAS and LAPACKE. Sturla
10.03.2012 23:19, Sturla Molden kirjoitti: [clip]
If the intention is to avoid Fortran dependency in NumPy, I am not sure why this is better than a dependency on CBLAS and LAPACKE.
The CBLAS parts aren't compiled if the library is not available --- in that case Numpy just falls back to a totally naive implementation of the matrix multiplication. Pauli
Hi Sturla, Quoting Sturla Molden <sturla@molden.no>:
Den 10.03.2012 22:56, skrev Sturla Molden:
I am not sure why NumPy uses f2c'd routines instead of a dependency on BLAS and LAPACK like SciPy.
Actually, np.dot does depend on the CBLAS interface to BLAS (_dotblas.c).
But the lapack methods in lapack_lite seems to use f2c'd code. I am not sure if they will use an optimized BLAS or just link to f2c's BLAS in blas_lite.c.
If the intention is to avoid Fortran dependency in NumPy, I am not sure why this is better than a dependency on CBLAS and LAPACKE.
Thanks for the more complete information. Now I understand better why it is more difficult to access the underlying libraries when using numpy instead of scipy. My main objective was to avoid having to ship libraries with my python extension modules. Using those already available via numpy seemed to me the most natural option since my C extension depends on numpy but not on scipy. In addition, people using my extension module could benefit from better libraries than those I would be able to ship. Eliminating the fortran dependency appeared to me as an added bonus: I have a 64 bit intel fortran compiler license for the only reason of not being limited by a fortran dependency. Best regards, Armando
participants (5)
-
"V. Armando Solé" -
David Cournapeau -
Pauli Virtanen -
sole＠esrf.fr
-
Sturla Molden