您可能有兴趣使用scipy.linalg.cython_lapack。
它提供对 LAPACK 函数dgetri 等的访问。而the good news 是:
这使得使用来自任何第 3 方的 SciPy 的 BLAS 和 LAPACK 成为可能
Cython 模块没有明确地与库链接。这表示
scikit-learn 和 statsmodels 等项目不需要维护
对 BLAS 和 LAPACK 的单独构建依赖。
Calling BLAS / LAPACK directly using the SciPy interface and Cython 提供了使用 dger 的示例。另见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
优化后的函数可以在import 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)
#workspace query
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)
#allocation workspace
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文件进行cythonize和编译,可以使用以下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 myinverse 函数chainng LAPACK 的dgetrf() 和dgetri() 在python 主文件中被调用:
import numpy as np
import myinverse
n=42
#A=np.zeros((n,n))
#for i in range(n):
# A[i,i]=10
A=np.random.rand(n,n)
#A=np.zeros((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