Внедрение более быстрого внутреннего продукта python с BLAS

Я нашел это полезное руководство по использованию низкоуровневых BLAS-функций (реализованных в Cython), чтобы получить большие улучшения скорости по сравнению с стандартными процедурами линейной алгебры numpy в python. Теперь я успешно получил работу векторного продукта. Сначала я сохраняю следующее как linalg.pyx:

import cython
import numpy as np
cimport numpy as np

from libc.math cimport exp
from libc.string cimport memset

from scipy.linalg.blas import fblas

REAL = np.float64
ctypedef np.float64_t REAL_t

cdef extern from "/home/jlorince/flda/voidptr.h":
    void* PyCObject_AsVoidPtr(object obj)

ctypedef double (*ddot_ptr) (const int *N, const double *X, const int *incX, const double *Y, const int *incY) nogil
cdef ddot_ptr ddot=<ddot_ptr>PyCObject_AsVoidPtr(fblas.ddot._cpointer)  # vector-vector multiplication 

cdef int ONE = 1
def vec_vec(syn0, syn1, size):
    cdef int lSize = size
    f = <REAL_t>ddot(&lSize, <REAL_t *>(np.PyArray_DATA(syn0)), &ONE, <REAL_t *>(np.PyArray_DATA(syn1)), &ONE)
    return f

(исходный код для voidptr.h доступен здесь)

Как только я его скомпилирую, он отлично работает и определенно быстрее, чем np.inner:

In [1]: import linalg
In [2]: import numpy as np
In [3]: x = np.random.random(100)
In [4]: %timeit np.inner(x,x)
1000000 loops, best of 3: 1.61 µs per loop
In [5]: %timeit linalg.vec_vec(x,x,100)
1000000 loops, best of 3: 483 ns per loop
In [8]: np.all(np.inner(x,x)==linalg.vec_vec(x,x,100))
Out[8]: True

Теперь это замечательно, но работает только для случая вычисления точечного/внутреннего произведения (эквивалентного в этом случае) двух векторов. Что мне нужно сделать сейчас, реализовать аналогичные функции (которые, я надеюсь, предложит аналогичные ускорения) для выполнения векторных матричных внутренних продуктов. То есть, я хочу реплицировать функции np.inner при передаче 1D-массива и 2D-матрицы:

In [4]: x = np.random.random(5)
In [5]: y = np.random.random((5,5))
In [6]: np.inner(x,y)
Out[6]: array([ 1.42116225,  1.13242989,  1.95690196,  1.87691992,  0.93967486])

Это эквивалентно вычислению внутренних/точечных продуктов (опять же, эквивалентных для 1D массивов) 1D-массива и каждой строки матрицы:

In [32]: [np.inner(x,row) for row in y]
Out[32]:
[1.4211622497461549, 1.1324298918119025, 1.9569019618096966,1.8769199192990056, 0.93967485730285505]

Из того, что я видел в документации BLAS, я думаю, мне нужно начать с чего-то вроде этого (используя dgemv):

ctypedef double (*dgemv_ptr) (const str *TRANS, const int *M, const int *N, const double *ALPHA, const double *A, const int *LDA, const double *X, const int *incX, const double *BETA, const double *Y, const int *incY)
cdef dgemv_ptr dgemv=<dgemv>PyCObject_AsVoidPtr(fblas.dgemv._cpointer)  # matrix vector multiplication

Но мне нужна помощь (а), определяющая фактическую функцию, которую я могу использовать в Python (т.е. функцию vec-matrix, аналогичную выше vec_vec), и (b) знание того, как использовать ее для правильной репликации поведения np.inner, что мне нужно для модели, которую я реализую.

Изменить: Здесь - ссылка на соответствующую документацию BLAS для dgemv, которую мне нужно использовать, что подтверждается здесь:

In [13]: np.allclose(scipy.linalg.blas.fblas.dgemv(1.0,y,x), np.inner(x,y))
Out[13]: True

Но использование этого из коробки, как это, на самом деле медленнее, чем чистый np.inner, поэтому мне все еще нужна помощь в реализации Cython.

Edit2 Здесь моя последняя попытка, которая компилируется отлично, но сбой при разбиении на python с ошибкой сегментации, когда я пытаюсь запустить его:

cdef int ONE = 1
cdef char tr = 'n'
cdef REAL_t ZEROF = <REAL_t>0.0
cdef REAL_t ONEF = <REAL_t>1.0
def mat_vec(mat,vec,mat_rows,mat_cols):
    cdef int m = mat_rows
    cdef int n = mat_cols
    out = <REAL_t>dgemv(&tr, &m, &n, &ONEF, <REAL_t *>(np.PyArray_DATA(mat)), &m, <REAL_t *>(np.PyArray_DATA(vec)), &ONE, &ZEROF, NULL, &ONE)
    return out

После компиляции я пытаюсь запустить linalg.mat_vec(y,x,5,5) (используя те же x и y, что и выше), но это просто сбой. Я думаю, что я рядом, но не знаю, что еще изменить...

Ответ 1