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
