Passing multidimensional memoryviews to c-function

Question

The answers in the following thread did not help mee since I have multidimensional arrays Passing memoryview to C function.

I have the following test code:

Header file:

//cm.h
double per(double **marr,int rows,int cols);

File with C-function

//cm.c
#include <stdio.h>
double per(double **marr,int rows,int cols){
    int i,j;
    for (i=0;i<rows;i++){
        for (j=0;j<cols;j++){
            //Just update the array
            marr[i][j] = (double)(i+1)*(j+1);
        }
    }
}

Cython file:

#m.pyz
import numpy as np
cimport numpy as np
from cython.view cimport array as cvarray

cdef extern from 'cm.h':
    double per(double **marr,int rows,int cols);

def test(double[:,::1] x):
    cdef int rows = x.shape[0]
    cdef int cols = x.shape[1]
    per(x,rows,cols)
    for i in range(rows):
        for j in range(cols):
            print x[i,j]

And error message:

Error compiling Cython file:
------------------------------------------------------------
...
    double per(double **marr,int rows,int cols);

def test(double[:,::1] x):
    cdef int rows = x.shape[0]
    cdef int cols = x.shape[1]
    per(x,rows,cols)
        ^
------------------------------------------------------------

m.pyx:12:9: Cannot assign type 'double[:, ::1]' to 'double **'

I have read that typed memoryviews is the most modern way of handling Python arrays in Cython, but I cannot find out how to do this. I have some numerical recipes in C which operates on dynamicall made large multidimensional arrays. Is what I try to do a totally wrong way?

Answer 1

Internally the memoryview is actually stored as a 1D array together with some information about the sizes of the dimensions. See http://docs.cython.org/src/userguide/memoryviews.html#brief-recap-on-c-fortran-and-strided-memory-layouts

(As a slight side note you can have memoryviews with "indirect" dimensions which do store things as pointer to pointers. This only makes sense if they're a view of something that already has memory allocated like that - for example if you build a 2D array like that in C. You won't get those from (say) numpy objects so I'll ignore this detail).

You change the C to be

// pass a 1D array, and then calculate where we should be looking in it
// based on the stride
double per(double *marr,int rows,int cols, int row_stride, int col_stride){
    int i,j;
    for (i=0;i<rows;i++){
        for (j=0;j<cols;j++){
            //Just update the array
            marr[row_stride*i + col_stride*j] = (double)(i+1)(j+1);
        }
    }
}

The Cython code then needs changes to pass it the strides (which is stored in bytes, so divide by the itemsize to get the strides in "number of doubles" which C expects), together with the address of the first element

// also update the cdef defintion for per...
per(&x[0,0],x.shape[0],x.shape[1],
     x.strides[0]/x.itemsize,x.strides[1]/x.itemsize)