[Numpy-discussion] (2012) Accessing LAPACK and BLAS from the numpy C API

[Sturla Molden]

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>