python - calling dot products and linear algebra operations in Cython? -

i'm trying use dot products, matrix inversion , other basic linear algebra operations available in numpy cython. functions numpy.linalg.inv (inversion), numpy.dot (dot product), x.t (transpose of matrix/array). there's large overhead calling numpy.* cython functions , rest of function written in cython, i'd avoid this.

if assume users have numpy installed, there way like:

#include "numpy/npy_math.h" 

as extern, , call these functions? or alternatively call blas directly (or whatever numpy calls these core operations)?

to give example, imagine have function in cython many things , in end needs make computation involving dot products , matrix inverses:

cdef myfunc(...):   # ... many things faster python   # ...   # compute 1 value using dot products , inv   # without using    #   import numpy np    #   np.*   val = gammaln(sum(v)) - sum(gammaln(v)) + dot((v - 1).t, log(x).t) 

how can done? if there's library implements these in cython already, can use that, have not found anything. if procedures less optimized blas directly, not having overhead of calling numpy python module cython still make things overall faster.

example functions i'd call:

• dot product (np.dot)
• matrix inversion (np.linalg.inv)
• matrix multiplication
• taking transpose (equivalent of x.t in numpy)
• gammaln function (like scipy.gammaln equivalent, should available in c)

i realize says on numpy mailing list (https://groups.google.com/forum/?fromgroups=#!topic/cython-users/xzjmvsiqnte) if call these functions on large matrices, there no point in doing cython, since calling numpy result in majority of time spent in optimized c code numpy calls. however, in case, have many calls these linear algebra operations on small matrices -- in case, overhead introduced repeatedly going cython numpy , cython far outweigh time spent computing operation blas. therefore, i'd keep @ c/cython level these simple operations , not go through python.

i'd prefer not go through gsl, since adds dependency , since it's unclear if gsl actively maintained. since i'm assuming users of code have scipy/numpy installed, can safely assume have associated c code goes along these libraries, want able tap code , call cython.

edit: found library wraps blas in cython (https://github.com/tokyo/tokyo) close not i'm looking for. i'd call numpy/scipy c functions directly (i'm assuming user has these installed.)

calling blas bundled scipy "fairly" straightforward, here's 1 example calling dgemm compute matrix multiplication: https://gist.github.com/pv/5437087 note blas , lapack expect arrays fortran-contiguous (modulo lda/b/c parameters), hence order="f" , double[::1,:] required correct functioning.

computing inverses can done applying lapack function dgesv on identity matrix. signature, see here. requires dropping down rather low-level coding, need allocate temporary work arrays etc etc. --- these can encapsulated own convenience functions, or reuse code tokyo replacing lib_* functions function pointers obtained scipy in above way.

if use cython's memoryview syntax (double[::1,:]) transpose same x.t usual. alternatively, can compute transpose writing function of own swaps elements of array across diagonal. numpy doesn't contain operation, x.t changes strides of array , doesn't move data around.

it possible rewrite tokyo module use blas/lapack exported scipy , bundle in scipy.linalg, from scipy.linalg.blas cimport dgemm. pull requests accepted if wants down it.

---

as can see, boils down passing function pointers around. alluded above, cython in fact provide own protocol exchanging function pointers. example, consider from scipy.spatial import qhull; print(qhull.__pyx_capi__) --- functions accessed via from scipy.spatial.qhull cimport xxxx in cython (they're private though, don't that).

however, @ present, scipy.special not offer c-api. in fact quite simple provide it, given interface module in scipy.special written in cython.

i don't think there @ moment sane , portable way access function doing heavy lifting gamln, (although snoop around ufunc object, that's not sane solution :), @ moment it's best grab relevant part of source code scipy.special , bundle project, or use e.g. gsl.