c malloc array pointer return in cython

Question

How does one to return a malloc array pointer (or numpy array pointer) in cython back to python3, efficiently.

The cython code works perfectly as long as I don't return the array pointer

I would like:

def double complex* randn_zig(int n):
  ...
  r = malloc(n*n*sizeof(double complex))
  ...
  return r

The c11 (gcc 11) equivalent is:

double complex* randn_zig(int n){

    r = malloc(n*n*sizeof(double complex))

    return r
}

I have tried <double complex*> randn_zig(int n):

and randn_zig(<double complex*> r, int n):

and other permutations without success so far. The c and cython code version is 5 times as fast as Numby/ pylab randn version if I can find a way to return a pointer to the large 10^6 to 10^10 double complex array.

Hawkbill

Can you provide more code for the part when you are trying to call this function? And which errors do you get? – Lundgren

There is no python equivalent to a pointer so it's not really the pointer you want to return, it's a some reference to (ie, some python object that has a pointer to) the memory block that's been allocated. Which kind of python object you choose will depend on what you intend to do with the pointer/memory block. – Postbellum

Answer 1

Numpy C API

Your question is similar to this post.

You can use the function below to pass a C pointer to Numpy array. The memory will be freed automatically when the Numpy array is recycled. If you want free the pointer mamully, you should not set NPY_OWNDATA flag.

import numpy as np
cimport numpy as np

cdef pointer_to_numpy_array_complex128(void * ptr, np.npy_intp size):
    '''Convert c pointer to numpy array.
    The memory will be freed as soon as the ndarray is deallocated.
    '''
    cdef extern from "numpy/arrayobject.h":
        void PyArray_ENABLEFLAGS(np.ndarray arr, int flags)
    cdef np.ndarray[np.complex128, ndim=1] arr = \
            np.PyArray_SimpleNewFromData(1, &size, np.NPY_COMPLEX128, ptr)
    PyArray_ENABLEFLAGS(arr, np.NPY_OWNDATA)
    return arr

For reference:

• PyArray_SimpleNewFromData
• Numpy Data Type API

Cython Typed Memoryviews

Of couse, you can also use cython memoryview.

import numpy as np
cimport numpy as np

cdef np.complex128_t[:,:] view = <np.complex128_t[:n,:n]> c_pointer
numpy_arr = np.asarray(view)

The code above will transfer C pointer to a numpy array. However this would not free memory automaticlly, you have to free the memory by yourself or it would lead to memory leak!

Ericaericaceous

1. Syrtis thanks will try shortly. I really need to resolve this now for a project. Please see comment below. – Hawkbill

Answer 2

A further option (in addition to the two options from the top answer: PyArray_SimpleNewFromData and just returning the typed memoryview without handling the memory) is to use the cython.view.array class.

This is a fairly low-level class that can be used to wrap existing memory. It has an attribute callback_free_data where you can set a function to be called on destruction so that it does free the memory (example code here is copied from the documentation):

cdef view.array my_array = view.array(..., mode="fortran", allocate_buffer=False)
my_array.data = <char *> my_data_pointer

# define a function that can deallocate the data (if needed)
my_array.callback_free_data = free

It exposes the buffer protocol so that you can index it, use it with typed memoryviews, or wrap it with a Numpy array (without copying) with np.asarray. The latter feature may be easier to use than PyArray_SimpleNewFromData.

Hawger

Answer 3

I think the best approach is to pass the pointer of an existing array created in Python via NumPy to Cython, otherwise it seems you have to copy the content of the array created by malloc to another array, like demonstrated in this toy example:

import numpy as np
cimport numpy as np

from libc.stdlib cimport malloc, free

def main():
  cdef int i, n=40
  cdef double complex *r
  cdef np.ndarray[np.complex128_t, ndim=1] a
  a = np.zeros(n*n, dtype=np.complex128)
  r = <double complex *>malloc(n*n*sizeof(double complex))
  for i in range(n*n):
      r[i] = 1.
  for i in range(n*n):
      a[i] = r[i]
  free(r)
  return a

Bozcaada

I have been busy on other projects, and am now back resolving this one. Thanks. In any event, python 3.4 and numpy appear to include the Meresenne prime generator which I assume is dsfmt(), default in Ubuntu 15.10. I haven't looked at the Python source for ppython3.4, but it was not there for 3.2. I am not sure if they have the Ziggurat or the Leva algorithm for randn(). I need to look. I will implement one if desired. 73 – Hawkbill

Answer 4

For gcc 5+ using the C-11 standard ( gcc -std=gnu11 . . . ), the syntax for multi-dimensional malloc and calloc arrays has changed significantly.

A main() procedure, for creating a 2-D, double, complex calloc array r[n][n] for n = 1024 is now:

long n = 1024;
complex double (*r)[n] = calloc(n, sizeof *r);

An example of a Gaussian random number generator randn_box_muller() using a pointer to this calloc array r[n][n] is:

inline static void randn_box_muller(long n, complex double r[][n])
{
    long i, j; 
    register double x, y;

    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){  
            x = 2.*M_PI*dsfmt_genrand_close_open(&dsfmt);
            y = sqrt(-2.*log(dsfmt_genrand_close_open(&dsfmt)));
            r[i][j] = (cos(x) + I*sin(x))*y;
        }
     }
     return;
}

This relatively new calloc allocation syntax is a bit strange. It works well for 1, 2 and even n dimensional calloc and malloc arrays. Hopefully this will also work in conjunction with Python3. I hope to be testing this shortly.

Hawkbill