BLAS dger以fortran顺序填充矩阵？答案

【问题标题】：BLAS dger filling matrix in fortran order?BLAS dger以fortran顺序填充矩阵？
【发布时间】：2017-02-23 18:52:22
【问题描述】：

我正在使用 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.]])

【问题讨论】：

我把A = np.zeros([n, m])换成A = np.zeros([n, m], order='f')，得到了你想要的结果。
作为旁注，您可以将test_dger 函数参数double[:, :] A 更改为double[::1, :] A，以明确期望一个fortran 连续数组。如果没有，ValueError 将被提升。
@oz1 谢谢，我在某个时候尝试过（注意到拆解第一个结果与转置第二个结果相同）。但是，我正在操作的数组是 C 排序的，所以我不能这样做。
嗯，dger实际上是一个fortran子程序，它需要fortran连续数组，如果你只想使用c连续数组，你可以在cython中尝试cblas。
stackoverflow.com/questions/7688304/…

标签： python cython blas cblas

【解决方案1】：

这确实是 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.]])

【讨论】：