我试图使用点产品，矩阵求逆，并且是从用Cython numpy的其他可用的基本线性代数运算。 类似功能numpy.linalg.inv （反转）， numpy.dot （点积）， Xt （转置矩阵/阵列）。 有一个大的开销调用numpy.*从用Cython函数和函数的其余部分是写在用Cython，所以我想避免这种情况。
如果我假设用户有numpy安装，有没有办法做这样的事情：
#include "numpy/npy_math.h"
作为extern ，并调用这些功能呢？ 或可替代直接调用BLAS（或者不管它是什么，对这些核心业务numpy的电话）？
举个例子，假设你在用Cython的功能，做很多事情，到底需要作出计算，涉及的点积和矩阵求逆：
cdef myfunc(...):
# ... do many things faster than Python could
# ...
# compute one value using dot products and inv
# without using
# import numpy as np
# np.*
val = gammaln(sum(v)) - sum(gammaln(v)) + dot((v - 1).T, log(x).T)
如何才能做到这一点？ 如果有，在用Cython已经实现了这些库，我也可以使用，但没有发现任何东西。 即使这些程序比BLAS优化程度不够直接，没有调用的开销numpy距离用Cython Python模块仍将使事情更快的整体。
例如功能我想打电话：
- 点积（
np.dot ） - 矩阵求逆（
np.linalg.inv ） - 矩阵乘法
- 服用转置（相当于
xT在numpy的） - GAMMALN函数（像
scipy.gammaln当量，这应该在C提供）
我知道，因为它说numpy的邮件列表（ https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE ），如果你调用大型矩阵的这些功能，有在任何点从用Cython这样做，因为从numpy的调用它只会导致多数在优化的C代码，与NumPy调用花费的时间。 然而，在我的情况，我有很多呼吁对小矩阵这些线性代数运算 -在这种情况下，开销由用Cython反复回到numpy的和回用Cython将远远超过所花的时间实际上是从计算操作介绍BLAS。 因此，我希望把一切都在C /用Cython水平，这些简单的操作，而不是经过蟒蛇。
我不想去通过GSL，因为这增加了另一个依赖，因为目前还不清楚，如果GSL积极维护。 由于我假设的代码的用户已经SciPy的/ numpy的安装，我可以有把握地认为他们拥有所有关联的C代码，使用这些库一起去，所以我只是希望能够挖掘到的代码，并调用它从用Cython。
编辑 ：我发现，包装BLAS在用Cython库（ https://github.com/tokyo/tokyo ），这是接近，但不是我要找的。 我想直接调用numpy的/ SciPy的C函数（我假设用户已经安装了这些。）
调用带有SciPy的捆绑BLAS是“相当”直白，这里的调用DGEMM计算矩阵乘法一个例子： https://gist.github.com/pv/5437087注意，BLAS和LAPACK期望所有的数组为Fortran的连续（模将LDA / b / C参数），因此， order="F"和double[::1,:]其所需的正确功能。
计算逆可以通过应用LAPACK功能可以类似地完成dgesv上单位矩阵。 对于签名，看这里 。 所有这些都需要下降到相当低的水平编码，您需要分配临时工作阵列自己等等等等---但是，这些可以被封装到自己方便的功能，或者只是重用的代码tokyo通过更换lib_*功能与以上述方式从SciPy的获得函数指针。
如果你用用Cython的memoryview语法（ double[::1,:]你调换相同xT如常。 或者，您可以通过编写自己的函数，交换的对角线的数组元素计算转置。 NumPy的实际上不包含此操作， xT仅改变阵列的步伐，不会四处移动数据。
或许，这将有可能重写tokyo模块使用BLAS / LAPACK通过SciPy的出口，并捆绑其scipy.linalg ，这样你可以只是做from scipy.linalg.blas cimport dgemm 。 引入请求如果有人想踏踏实实地它被接受。
正如你所看到的，这一切都归结到各地传递函数指针。 如上面提到的，用Cython事实上确实提供它自己的协议用于交换函数指针。 举一个例子，可以考虑from scipy.spatial import qhull; print(qhull.__pyx_capi__) from scipy.spatial import qhull; print(qhull.__pyx_capi__) ---这些功能可以通过以下方式访问from scipy.spatial.qhull cimport XXXX在用Cython（他们是私人的，所以不这样做）。
然而，目前， scipy.special不提供此C-API。 它将但是实际上是非常简单的，提供它，因为在scipy.special接口模块是写在用Cython。
我不认为这是目前进入功能做了繁重的任务，任何理智的和可移植的方法gamln ，（虽然你可以围绕UFunc对象窥探，但是这不是一个明智的解决方案:)，所以目前它的可能最好还是抓住从scipy.special源代码的相关部分，并与项目捆绑，或使用如GSL。
也许最简单的方式，如果你使用GSL是使用这个GSL-接受>用Cython接口https://github.com/twiecki/CythonGSL ，并从那里调用BLAS（见例如https://github.com/twiecki /CythonGSL/blob/master/examples/blas2.pyx ）。 还应该照顾Fortran语言对C的顺序。 不是有很多新的GSL的功能，但你可以安全地假设它是积极的维护。 相较于东京CythonGSL较为齐全; 例如，它的特点是在numpy的缺席对称矩阵产品。
正如我刚才遇到同样的问题，并写了一些额外的功能，我会在这里包括他们在情况下别人发现它们非常有用。 我编写了一些矩阵乘法，并且也调用LAPACK函数矩阵求逆，行列式和Cholesky分解。 但是，你应该考虑尝试做线性代数的东西任何循环外，如果您有任何和我一样在这里 。 顺便说一句，这里的决定因素功能不够，如果您有任何建议工作。 另外，请注意，我没有做任何检查，看是否输入是顺应性的。
from scipy.linalg.cython_lapack cimport dgetri, dgetrf, dpotrf
cpdef void double[:, ::1] inv_c(double[:, ::1] A, double[:, ::1] B,
double[:, ::1] work, double[::1] ipiv):
'''invert float type square matrix A
Parameters
----------
A : memoryview (numpy array)
n x n array to invert
B : memoryview (numpy array)
n x n array to use within the function, function
will modify this matrix in place to become the inverse of A
work : memoryview (numpy array)
n x n array to use within the function
ipiv : memoryview (numpy array)
length n array to use within function
'''
cdef int n = A.shape[0], info, lwork
B[...] = A
dgetrf(&n, &n, &B[0, 0], &n, &ipiv[0], &info)
dgetri(&n, &B[0,0], &n, &ipiv[0], &work[0,0], &lwork, &info)
cpdef double det_c(double[:, ::1] A, double[:, ::1] work, double[::1] ipiv):
'''obtain determinant of float type square matrix A
Notes
-----
As is, this function is not yet computing the sign of the determinant
correctly, help!
Parameters
----------
A : memoryview (numpy array)
n x n array to compute determinant of
work : memoryview (numpy array)
n x n array to use within function
ipiv : memoryview (numpy array)
length n vector use within function
Returns
-------
detval : float
determinant of matrix A
'''
cdef int n = A.shape[0], info
work[...] = A
dgetrf(&n, &n, &work[0,0], &n, &ipiv[0], &info)
cdef double detval = 1.
cdef int j
for j in range(n):
if j != ipiv[j]:
detval = -detval*work[j, j]
else:
detval = detval*work[j, j]
return detval
cdef void chol_c(double[:, ::1] A, double[:, ::1] B):
'''cholesky factorization of real symmetric positive definite float matrix A
Parameters
----------
A : memoryview (numpy array)
n x n matrix to compute cholesky decomposition
B : memoryview (numpy array)
n x n matrix to use within function, will be modified
in place to become cholesky decomposition of A. works
similar to np.linalg.cholesky
'''
cdef int n = A.shape[0], info
cdef char uplo = 'U'
B[...] = A
dpotrf(&uplo, &n, &B[0,0], &n, &info)
cdef int i, j
for i in range(n):
for j in range(n):
if j > i:
B[i, j] = 0
cpdef void dotmm_c(double[:, :] A, double[:, :] B, double[:, :] out):
'''matrix multiply matrices A (n x m) and B (m x l)
Parameters
----------
A : memoryview (numpy array)
n x m left matrix
B : memoryview (numpy array)
m x r right matrix
out : memoryview (numpy array)
n x r output matrix
'''
cdef Py_ssize_t i, j, k
cdef double s
cdef Py_ssize_t n = A.shape[0], m = A.shape[1]
cdef Py_ssize_t l = B.shape[0], r = B.shape[1]
for i in range(n):
for j in range(r):
s = 0
for k in range(m):
s += A[i, k]*B[k, j]
out[i, j] = s
