[cython-users] numpy with multidimensional array

[Wonjun, Choi]

Hi

I have checked example which is wrapping c by using numpy.
there is two converting function( function for converting numpy to c,
function for converting c to numpy)
and they both use malloc function.
but if I do not use memory copy, how can I make the code?

----------------------------------------------------------------------------------
def matmul1(np.ndarray[DTYPE_t, ndim=2] a, np.ndarray[DTYPE_t, ndim=2]
b):
''' Matrix multiplication. Takes two square Float32 numpy arrays.
'''

cdef int N = a.shape[0]
cdef int i
cdef float **a_c
cdef float **b_c
cdef float **res

# check if square arrays:
if a.shape[1] != N or b.shape[0] != N or b.shape[1] != N:
raise ValueError, 'matmul1: need square arrays for
multiplication!'

# check if contiguous, if not force C contiguous arrays
if not (<object>a).flags["C_CONTIGUOUS"]:
a = a.copy('C')
if not (<object>b).flags["C_CONTIGUOUS"]:
b = b.copy('C')

# convert using the function
a_c = npy2c_float(a)
b_c = npy2c_float(b)

# allocate res
res = <float **> malloc(N*sizeof(float*))
for i in range(N):
res[i] = <float *> malloc(N * sizeof(float))

matmul(a_c,b_c,res,N)

free(a_c)
free(b_c)

# convert to numpy array and free res
result = c2npy_float(res,N,N)

return result
----------------------------------------------------------------------------------

[Robert Bradshaw]

On Thu, Sep 8, 2011 at 5:56 AM, Wonjun, Choi <wonjun...@gmail.com> wrote:
> Hi
>
> I have checked example which is wrapping c by using numpy.
> there is two converting function( function for converting numpy to c,
> function for converting c to numpy)
> and they both use malloc function.
> but if I do not use memory copy, how can I make the code?

If you view your float** as a pointer to rows, there's no need to copy
the rows themselves, just point to the original data in the numpy
array. E.g. (untested)

cdef float** npy2c_float2d(np.ndarray[float, ndim=2, mode=c] a):
cdef float** a_c = malloc(a.shape[0] * sizeof(float*))
for k in range(a.shape[0]):
a_c[k] = &a[k, 0]
return a_c

will do the conversion without copying any data. (The caller is
responsible for freeing the C array.) To return the result, do
something like

def mul(a, b):
# any checks you need to do
try:
res = np.ndarray((N, N), dtype=float)
cdef float** a_c = npy2c_float2d(a)
cdef float** b_c = npy2c_float2d(b)
cdef float** res_c = npy2c_float2d(res)
matmul(a_c, b_c, res_c, N)
return res
finally:
free(a_c)
free(b_c)
free(res_c)

- Robert

[최원준]

2011/9/9 Robert Bradshaw <robe...@math.washington.edu>

On Thu, Sep 8, 2011 at 5:56 AM, Wonjun, Choi <wonjun...@gmail.com> wrote:
> Hi
>
> I have checked example which is wrapping c by using numpy.
> there is two converting function( function for converting numpy to c,
> function for converting c to numpy)
> and they both use malloc function.
> but if I do not use memory copy, how can I make the code?

If you view your float** as a pointer to rows, there's no need to copy
the rows themselves, just point to the original data in the numpy
array. E.g. (untested)

===> it means if the souce is fortran, I need memcpy function and contiguous flag?

 
cdef float** npy2c_float2d(np.ndarray[float, ndim=2, mode=c] a):

   cdef float** a_c = malloc(a.shape[0] * sizeof(float*))
   for k in range(a.shape[0]):
       a_c[k] = &a[k, 0]
   return a_c

will do the conversion without copying any data. (The caller is
responsible for freeing the C array.) To return the result, do
something like

def mul(a, b):
   # any checks you need to do
   try:
       res = np.ndarray((N, N), dtype=float)
       cdef float** a_c = npy2c_float2d(a)
       cdef float** b_c = npy2c_float2d(b)
       cdef float** res_c = npy2c_float2d(res)
       matmul(a_c, b_c, res_c, N)
       return res
   finally:

       free(a_c)
       free(b_c)

       free(res_c)

- Robert

> ----------------------------------------------------------------------------------

===> I tested below code but I got an error.

---------------------------------------------------------------------------
import numpy as np
cimport numpy as np

cdef float **a_c
cdef float **b_c

cdef float **tmp_c
cdef float **res_c

cdef float** npy2c_float2d(np.ndarray[float, ndim=2, mode=c] a):

   tmp_c = malloc(a.shape[0] * sizeof(float*))

   for k in range(a.shape[0]):

       tmp_c[k] = &a[k, 0]
   return tmp_c

def mul(a, b):
   # any checks you need to do
   try:
       res = np.ndarray((N, N), dtype=float)

       a_c = npy2c_float2d(a)
       b_c = npy2c_float2d(b)

       res_c = npy2c_float2d(res)
       matmul(a_c, b_c, res_c, N)
       return res
   finally:
       free(a_c)
       free(b_c)
       free(res_c)

error :

Error compiling Cython file:
------------------------------------------------------------
...

       matmul(a_c, b_c, res_c, N)
       return res
   finally:
       free(a_c)
       free(b_c)
       free(res_c)

                ^
------------------------------------------------------------

matmul_cython.pyx:27:17: Cannot convert 'float **' to Python object
make: *** [matmul_cython.c] Error 1


---------------------------------------------------------------------------

[Robert Bradshaw]

On Fri, Sep 9, 2011 at 1:04 AM, 최원준 <wonjun...@gmail.com> wrote:
>
>
> 2011/9/9 Robert Bradshaw <robe...@math.washington.edu>
>>
>> On Thu, Sep 8, 2011 at 5:56 AM, Wonjun, Choi <wonjun...@gmail.com>
>> wrote:
>> > Hi
>> >
>> > I have checked example which is wrapping c by using numpy.
>> > there is two converting function( function for converting numpy to c,
>> > function for converting c to numpy)
>> > and they both use malloc function.
>> > but if I do not use memory copy, how can I make the code?
>>
>> If you view your float** as a pointer to rows, there's no need to copy
>> the rows themselves, just point to the original data in the numpy
>> array. E.g. (untested)
>>
> ===> it means if the souce is fortran, I need memcpy function and contiguous
> flag?

Maybe. Maybe your NumPy arrays are already in fortran contiguous
order. Or maybe your fortran libraries can handle C-style arrays. Or
maybe it doesn't matter because the product of the transpose is the
transpose of the reverse product.

You didn't cimport or declare "free."

- Robert

[최원준]

I see, thank you.