BLAS dger filling matrix in fortran order?

Question

I am using cython to do a rank one update of a rectangular matrix A. I cannot get dger to do the update as I want, so I've isolated it in a function:

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)

which compiles just fine. However, doing:

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)

gives me

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

which has the desired n by m shape, but the values in the wrong order. I assum C order vs fortran order has something to do with this, but I have not been able to solve the problem by myself.

The result I'm expecting is what would be given by

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

Answer 1

This is C vs Fortran orders, indeed. Since matrix A dimensions in Fortran order is 4x3, x and y should be swapped. The solution is:

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

Now it works just fine:

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, x, y)
A
# array([[ 0.,  1.,  2.,  3.],
#       [ 0.,  4.,  8., 12.],
#       [ 0.,  0.,  0.,  0.]])