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.]])

I change the line `A = np.zeros([n, m])` to `A = np.zeros([n, m], order='f')`, and get the result you are expecting. — oz1, Feb 24 '17 at 03:15
As as side note, you can change your `test_dger` function parameter `double[:, :] A` to `double[::1, :] A`, for explicitly expecting an fortran contiguous array. If not, `ValueError ` will be raised. — oz1, Feb 24 '17 at 03:25
@oz1 Thank you, I tried that at some point (noticing that raveling the first result gives the same thing as transposing the second one them raveling it). However, the array I'm manipulating is C ordered, so I cannot do that. — P. Camilleri, Feb 24 '17 at 07:37
Well, the `dger` is actually a fortran subroutine, which expects fortran contiguous arrays, if you only want to work with c-contiguous array, you can try `cblas` in cython. — oz1, Feb 24 '17 at 08:10
http://stackoverflow.com/questions/7688304/how-to-force-numpy-array-order-to-fortran-style — DavidW, Feb 24 '17 at 08:26

score 1 · Answer 1 · edited May 05 '22 at 17:44

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.]])

BLAS dger filling matrix in fortran order?

1 Answers1