Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
Hi Sturla, Quoting Sturla Molden : > 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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
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]
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
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. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
On Tue, Mar 6, 2012 at 2:57 PM, 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. 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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
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
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
