您可能会对使用scipy.linalg.cython_lapack感兴趣。
它提供了访问LAPACK函数
dgetri
等的方法。而且,
好消息是:
这使得第三方Cython模块可以从SciPy的BLAS和LAPACK中获得访问权限,而无需显式地链接库。这意味着像scikit-learn和statsmodels这样的项目不需要维护BLAS和LAPACK的单独构建依赖关系。
一个使用dger
的示例可在Calling BLAS / LAPACK directly using the SciPy interface and Cython上找到。还可以查看Improving Cython Lapack performance with internal array definitions?。我在MPI python-Open-MPI的答案中详细介绍了如何使用cython_blas,因此以下是如何将其调整为dgetri:
代码的关键部分是用Cython
编写的,位于专用文件myinverse.pyx
中。
通过Cython将此文件转换为myinverse.c
文件。
此c文件由您喜欢的C编译器gcc
编译以构建共享库myinverse.so
。
优化的函数可以在导入myinverse之后在您的程序中使用。
这是一个Cython模块,应放置在.pyx文件中:
import numpy
cimport numpy
cimport scipy.linalg.cython_lapack
ctypedef numpy.float64_t DTYPE_t
cimport cython
from libc.stdlib cimport malloc, free
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def invert(numpy.ndarray[DTYPE_t, ndim=2] array):
cdef int rows = array.shape[0]
cdef int cols = array.shape[1]
cdef int info = 0
if cols !=rows:
return array,1,"not a square matrix"
cdef int* ipiv = <int *> malloc(rows * sizeof(int))
if not ipiv:
raise MemoryError()
scipy.linalg.cython_lapack.dgetrf(&cols,&rows,&array[0,0],&rows,ipiv,&info)
if info !=0:
free(ipiv)
return array,info,"dgetrf failed, INFO="+str(info)
cdef double workl
cdef int lwork=-1
scipy.linalg.cython_lapack.dgetri(&cols,&array[0,0],&rows,ipiv,&workl,&lwork,&info)
if info !=0:
free(ipiv)
return array,info,"dgetri failed, workspace query, INFO="+str(info)
lwork= int(workl)
cdef double* work = <double *> malloc(lwork * sizeof(double))
if not work:
raise MemoryError()
scipy.linalg.cython_lapack.dgetri(&cols,&array[0,0],&rows,ipiv,work,&lwork,&info)
if info !=0:
free(ipiv)
free(work)
return array,info,"dgetri failed, INFO="+str(info)
free(ipiv)
free(work)
return array,info,""
为了将.pyx文件转换为Cython并编译,可以使用以下makefile(希望您正在使用Linux...)
all: myinverse myinverseb
myinverse: myinverse.pyx
cython -a myinverse.pyx
myinverseb: myinverse.c
gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.7 -o myinverse.so myinverse.c
在主Python文件中调用新的Python myinverse函数,该函数使用LAPACK的dgetrf()和dgetri()进行链式操作:
import numpy as np
import myinverse
n=42
A=np.random.rand(n,n)
Am,info,string=myinverse.invert(A.copy())
if info==0:
print np.linalg.norm(A.dot(Am)-np.identity(n), np.inf)
else :
print "inversion failed, info=",info, string