Вызов BLAS / LAPACK напрямую с использованием интерфейса SciPy и Cython

Question

Здесь было сообщение об этом: https://gist.github.com/JonathanRaiman/f2ce5331750da7b2d4e9 которое показывает значительное улучшение скорости, просто вызывая библиотеки Фортрана (BLAS / LAPACK / Intel MKL / OpenBLAS / все, что вы установили с NumPy). После многих часов работы над этим (из-за устаревших библиотек SciPy) я, наконец, получил его для компиляции без результатов. Это было в 2 раза быстрее, чем NumPy. К сожалению, как указал другой пользователь, подпрограмма Fortran всегда добавляет выходную матрицу к новым вычисленным результатам, поэтому она соответствует NumPy только при первом запуске. Т.е. A := alpha*x*y.T + A, Так что это еще предстоит решить с помощью быстрого решения.

[ОБНОВЛЕНИЕ: ДЛЯ ТОГО, ЧТО ВЫ ИЩЕТЕ ИСПОЛЬЗОВАТЬ ИНТЕРФЕЙС SCIPY, ПРОСТО ПРОЙДИТЕ ЗДЕСЬ https://github.com/scipy/scipy/blob/master/scipy/linalg/cython_blas.pyx КАК ОНИ УЖЕ ОПТИМИЗИРУЛИ ЗВОНОКИ ДЛЯ BLAS / LAP В ЗАЯВЛЕНИЯХ CPDEF ПРОСТО КОПИРУЙТЕ / ВСТАВЛЯЙТЕ В СВОЙ ЦИФОННЫЙ СКРИПТ # Python-accessible wrappers for testing: Также по ссылке выше cython_lapack.pyx доступен, но не имеет тестовых скриптов Cython]

Тестовый сценарий

import numpy as np;
from cyblas import outer_prod;
a=np.random.randint(0,100, 1000);
b=np.random.randint(0,100, 1000);
a=a.astype(np.float64)
b=b.astype(np.float64)
cy_outer=np.zeros((a.shape[0],b.shape[0]));
np_outer=np.zeros((a.shape[0],b.shape[0]));

%timeit outer_prod(a,b,cy_outer)
#%timeit outer_prod(a,b) #use with fixed version instead of above line, results will automatically update cy_outer
%timeit np.outer(a,b, np_outer)
100 loops, best of 3: 2.83 ms per loop
100 loops, best of 3: 6.58 ms per loop

# КОНЕЦ ИСПЫТАНИЯ

PYX-файл для компиляции cyblas.pyx (в основном версия np.ndarray)

import cython
import numpy as np
cimport numpy as np

from cpython cimport PyCapsule_GetPointer 
cimport scipy.linalg.cython_blas
cimport scipy.linalg.cython_lapack
import scipy.linalg as LA

REAL = np.float64
ctypedef np.float64_t REAL_t
ctypedef np.uint64_t  INT_t

cdef int ONE = 1
cdef REAL_t ONEF = <REAL_t>1.0

ctypedef void (*dger_ptr) (const int *M, const int *N, const double *alpha, const double *X, const int *incX, double *Y, const int *incY, double *A, const int * LDA) nogil
cdef dger_ptr dger=<dger_ptr>PyCapsule_GetPointer(LA.blas.dger._cpointer, NULL)  # A := alpha*x*y.T + A

cpdef outer_prod(_x, _y, _output):
#cpdef outer_prod(_x, _y): #comment above line & use this to use the reset output matrix to zeros
    cdef REAL_t *x = <REAL_t *>(np.PyArray_DATA(_x))
    cdef int M = _y.shape[0]
    cdef int N = _x.shape[0]
    #cdef np.ndarray[np.float64_t, ndim=2, order='c'] _output = np.zeros((M,N)) #slow fix to uncomment to reset output matrix to zeros
    cdef REAL_t *y = <REAL_t *>(np.PyArray_DATA(_y))
    cdef REAL_t *output = <REAL_t *>(np.PyArray_DATA(_output))
    with nogil:
        dger(&M, &N, &ONEF, y, &ONE, x, &ONE, output, &M)

Очень признателен. Надеюсь, это сэкономит время другим людям (это ПОЧТИ работает) - на самом деле, как я уже говорил, он работает 1x и соответствует NumPy, а затем каждый последующий вызов добавляет матрицу результатов снова. Если я сбрасываю выходную матрицу в 0 и перезапускаю результаты, то сопоставляю NumPy. Странно... хотя, если кто-то раскомментирует несколько строк выше, он будет работать, хотя только на скоростях NumPy. Альтернатива найдена memset и будет в другом посте... Я просто еще не выяснил, как именно это назвать.

Answer 1

По данным netlib dger(M, N, ALPHA, X INCX, Y, INCY, A, LDA) выполняет A := alpha*x*y**T + A, Так A должны быть все нули, чтобы получить внешнее произведение X а также Y,