0

我正在使用 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 x m 形状，但值的顺序错误。我认为 C 顺序与 fortran 顺序与此有关，但我自己无法解决问题。

我期待的结果是

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

1 回答 1

0

这确实是 C 与 Fortran 的命令。由于 Fortran 顺序中的矩阵A尺寸是 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.]])
于 2021-05-05T12:04:00.733 回答