Использование интерфейса Scipy cython_blas от Cython, не работающего на векторах Mx1 1xN

Question

Это должно иметь дело с подобной проблемой здесь: вызов BLAS / LAPACK напрямую с использованием интерфейса SciPy и Cython, но отличается, потому что я использую фактический код в примере SciPy здесь _test_dgemm: https://github.com/scipy/scipy/blob/master/scipy/linalg/cython_blas.pyx который очень быстрый (в 5 раз быстрее, чем numpy.dot когда выходная матрица вводится или в 20 раз быстрее, если нет). Он не дает результатов, если переданы векторы Mx1 1xN. Выдает те же значения, что и numpy.dot с пройденными матрицами. Я свернул код, так как ответы не были опубликованы для ясности. Вот dgemm.pyx.:

import numpy as np
cimport numpy as np
from scipy.linalg.cython_blas cimport dgemm
from cython cimport boundscheck

@boundscheck(False)
cpdef int fast_dgemm(double[:,::1] a, double[:,::1] b, double[:,::1] c, double alpha=1.0, double beta=0.0) nogil except -1:

    cdef:
        char *transa = 'n'
        char *transb = 'n'
        int m, n, k, lda, ldb, ldc
        double *a0=&a[0,0]
        double *b0=&b[0,0]
        double *c0=&c[0,0]

    ldb = (&a[1,0]) - a0 if a.shape[0] > 1 else 1
    lda = (&b[1,0]) - b0 if b.shape[0] > 1 else 1

    k = b.shape[0]
    if k != a.shape[1]:
        with gil:
            raise ValueError("Shape mismatch in input arrays.")
    m = b.shape[1]
    n = a.shape[0]
    if n != c.shape[0] or m != c.shape[1]:
        with gil:
            raise ValueError("Output array does not have the correct shape.")
    ldc = (&c[1,0]) - c0 if c.shape[0] > 1 else 1
    dgemm(transa, transb, &m, &n, &k, &alpha, b0, &lda, a0,
               &ldb, &beta, c0, &ldc)
    return 0

Вот пример тестового скрипта:

import numpy as np;

a=np.random.randn(1000);
b=np.random.randn(1000);
a.resize(len(a),1);
a=np.array(a, order='c');
b.resize(1,len(b)); 
b=np.array(b, order='c');
c = np.empty((a.shape[0],b.shape[1]), float, order='c'); 

from dgemm import _test_dgemm; 
_test_dgemm(a,b,c);

И если вы хотите поиграть с ним на Windows с Python 3.5 x64, вот setup.py построить его с помощью командной строки, набрав python setup.py build_ext --inplace --compiler=msvc

from Cython.Distutils import build_ext
import numpy as np
import os

try:
    from setuptools import setup
    from setuptools import Extension
except ImportError:
    from distutils.core import setup
    from distutils.extension import Extension

module = 'dgemm'

ext_modules = [Extension(module, sources=[module + '.pyx'],
              include_dirs=['C://Program Files (x86)//Windows Kits//10//Include//10.0.10240.0//ucrt','C://Program Files (x86)//Microsoft Visual Studio 14.0//VC//include','C://Program Files (x86)//Windows Kits//8.1//Include//shared'],
              library_dirs=['C://Program Files (x86)//Windows Kits//8.1//bin//x64', 'C://Windows//System32', 'C://Program Files (x86)//Microsoft Visual Studio 14.0//VC//lib//amd64', 'C://Program Files (x86)//Windows Kits//8.1//Lib//winv6.3//um//x64', 'C://Program Files (x86)//Windows Kits//10//Lib//10.0.10240.0//ucrt//x64'],
              extra_compile_args=['/Ot', '/favor:INTEL64', '/EHsc', '/GA'],
              language='c++')]

setup(
    name = module,
    ext_modules = ext_modules,
    cmdclass = {'build_ext': build_ext},
    include_dirs = [np.get_include(), os.path.join(np.get_include(), 'numpy')]
    )

Любая помощь высоко ценится!

Answer 1

Если я правильно понял, вы пытаетесь использовать fortran-рутины для массивов с c-memory-layout.

Даже если вам это очевидно известно, я хотел бы сначала уточнить порядок строк-мажоров (c-memory-layout) и порядок столбцов-мажоров (fortran-memory-layout), чтобы вывести мой ответ.,

Так что, если у нас есть 2x3 матрица (т.е. 2 строки и 3 столбца) A и сохранить его в некоторой непрерывной памяти, мы получим:

row-major-order(A) = A11, A12, A13, A21, A22, A23
col-major-order(A) = A11, A21, A12, A22, A13, A33

Это означает, что если мы получим непрерывную память, которая представляет матрицу в мажорном порядке строк, и интерпретируем ее как матрицу в мажорном столбце, мы получим совершенно другую матрицу!

Тем не менее, мы посмотрим на транспонированную матрицу A^t мы можем легко увидеть:

row-major-order(A) = col-major-order(A^t)
col-major-order(A) = row-major-order(A^t)

Это означает, что если мы хотим получить матрицу C в порядке мажорной строки как результат, подпрограмма blas должна записать транспонированную матрицу C в порядке-столбце-мажоре (после всего этого мы не можем изменить) в эту самую память. Тем не мение, C^t=(AB)^t=B^t*A^t а также B^t A^t являются исходными матрицами, переосмысленными в главном порядке столбцов.

Теперь пусть A быть n x k -матрица и B k x m -matrix вызов функции dgemm должен выглядеть следующим образом:

dgemm(transa, transb, &m, &n, &k, &alpha, b0, &m, a0, &k, &beta, c0, &m)

Как видите, вы переключили некоторые n а также m в вашем коде.