[Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API
Sturla Molden
sturla at molden.no
Sat Mar 10 16:27:42 EST 2012
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]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20120310/8c1d4d61/attachment.html>
More information about the NumPy-Discussion
mailing list