calling dot products and linear algebra operations in Cython?
Calling BLAS bundled with Scipy is "fairly" straightforward, here's one example for calling DGEMM to compute matrix multiplication: https://gist.github.com/pv/5437087 Note that BLAS and LAPACK expect all arrays to be Fortran-contiguous (modulo the lda/b/c parameters), hence order="F" and double[::1,:] which are required for correct functioning.
Computing inverses can be similarly done by applying the LAPACK function dgesv on the identity matrix. For the signature, see here. All this requires dropping down to rather low-level coding, you need to allocate temporary work arrays yourself etc etc. --- however these can be encapsulated into your own convenience functions, or just reuse the code from tokyo by replacing the lib_* functions with function pointers obtained from Scipy in the above way.
If you use Cython's memoryview syntax (double[::1,:]) you transpose is the same x.T as usual. Alternatively, you can compute the transpose by writing a function of your own that swaps elements of the array across the diagonal. Numpy doesn't actually contain this operation, x.T only changes the strides of the array and doesn't move the data around.
It would probably be possible to rewrite the tokyo module to use the BLAS/LAPACK exported by Scipy and bundle it in scipy.linalg, so that you could just do from scipy.linalg.blas cimport dgemm. Pull requests are accepted if someone wants to get down to it.
As you can see, it all boils down to passing function pointers around. As alluded to above, Cython does in fact provide its own protocol for exchanging function pointers. For an example, consider from scipy.spatial import qhull; print(qhull.__pyx_capi__) --- those functions could be accessed via from scipy.spatial.qhull cimport XXXX in Cython (they're private though, so don't do that).
However, at the present, scipy.special does not offer this C-API. It would however in fact be quite simple to provide it, given that the interface module in scipy.special is written in Cython.
I don't think there is at the moment any sane and portable way to access the function doing the heavy lifting for gamln, (although you could snoop around the UFunc object, but that's not a sane solution :), so at the moment it's probably best to just grab the relevant part of source code from scipy.special and bundle it with your project, or use e.g. GSL.
Perhaps the easiest way if you do accept using the GSL would be to use this GSL->cython interface https://github.com/twiecki/CythonGSL and call BLAS from there (see the example https://github.com/twiecki/CythonGSL/blob/master/examples/blas2.pyx). It should also take care of the Fortran vs C ordering.There aren't many new GSL features, but you can safely assume it is actively maintained. The CythonGSL is more complete compared to tokyo; e.g., it features symmetric-matrix products that are absent in numpy.
As I've just encountered the same problem, and wrote some additional functions, I'll include them here in case someone else finds them useful. I code up some matrix multiplication, and also call LAPACK functions for matrix inversion, determinant and cholesky decomposition. But you should consider trying to do linear algebra stuff outside any loops, if you have any, like I do here. And by the way, the determinant function here isn't quite working if you have suggestions. Also, please note that I don't do any checking to see if inputs are conformable.
from scipy.linalg.cython_lapack cimport dgetri, dgetrf, dpotrfcpdef 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 detvalcdef 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