BLAS заполнение матрицы джера в фортране?

Question

Я использую Cython для обновления первого уровня прямоугольной матрицы A. Я не могу заставить dger выполнить обновление так, как я хочу, поэтому я выделил его в функции:

from scipy.linalg.cython_blas cimport dger
cimport cython

def test_dger(double[:, :] A, double[:] x, double[:] y):
    cdef int inc = 1
    cdef double one = 1
    cdef int n_= A.shape[0]
    cdef int m = A.shape[1]
    dger(&n, &m, &one, &x[0], &inc, &y[0], &inc, &A[0, 0], &n)
    return np.array(A)

который компилируется просто отлично. Тем не менее, делая:

n = 3
m = 4
A = np.zeros([n, m])
y = np.arange(m, dtype=float)
x = np.array([1., 4, 0])
test_dger(A, x, y)

дает мне

array([[  0.,   0.,   0.,   1.],
       [  4.,   0.,   2.,   8.],
       [  0.,   3.,  12.,   0.]])

который имеет желаемую форму n на m, но значения в неправильном порядке. Я предполагаю, что C-Order против Fortran Order-то как-то связан с этим, но я не смог решить проблему сам.

Результат, который я ожидаю, - это то, что

np.dot(x[:, None], y[None, :])
array([[  0.,   1.,   2.,   3.],
       [  0.,   4.,   8.,  12.],
       [  0.,   0.,   0.,   0.]])

Answer 1

Это действительно заказы C против Fortran. Поскольку матрица A размеры в порядке Fortran - 4x3, x а также yнужно поменять местами. Решение:

      cimport cython
from scipy.linalg.cython_blas cimport dger

def test_dger(double[:, :] A, double[:] y, double[:] x):
    # x and y swapped!
    cdef int inc = 1
    cdef double one = 1.0
    cdef int m = A.shape[0] # Transposed shape!
    cdef int n = A.shape[1] # Transposed shape!
    dger(&n, &m, &one, &x[0], &inc, &y[0], &inc, &A[0, 0], &n)
    return A

Теперь все работает нормально:

      n, m = 3, 4
A = np.zeros([n, m])
x = np.array([1., 4, 0], dtype=float)
y = np.arange(m, dtype=float)

test_dger(A, y, x)
A
array([[ 0.,  1.,  2.,  3.],
       [ 0.,  4.,  8., 12.],
       [ 0.,  0.,  0.,  0.]])