Pass complex numpy array to C++ in Cython

Question

I want to Cythonize portion of a pyx script which involves work with numpy arrays with complex numbers. The relevant portion of the python script looks like this:

M = np.dot(N , Q)

In my work, N, Q and M are numpy arrays with complex number entries.

Specifically, I want to transfer the matrices N and Q to a C++ code and do the matrix multiplication in C++.

While I know the method to transfer real valued numpy arrays using pointers to C++ script, followed by use of cython, I am a bit confused about how I should approach things for numpy arrays with complex values.

This is how I am trying to transfer the array from pyx to C++ presently.

import numpy as np
cimport numpy as np

cdef extern from "./matmult.h" nogil:
    void mult(double* M, double* N, double* Q)

def sim():    
    cdef:
        np.ndarray[np.complex128_t,ndim=2] N = np.zeros(( 2 , 2 ), dtype=np.float64)
        np.ndarray[np.complex128_t,ndim=2] Q = np.zeros(( 2 , 2 ), dtype=np.float64)  
        np.ndarray[np.complex128_t,ndim=2] M = np.zeros(( 2 , 2 ), dtype=np.float64)          

    N = np.array([[1.1 + 2j,2.2],[3.3,4.4]])
    Q = np.array([[3.3,4.4+5j],[5.5,6.6]])  

    mult(&M[0,0], &N[0,0], &Q[0,0])
    print M

This is my C++ code:

#include "matmult.h"
using namespace std;

int main(){}

void mult(double *M, double *N, double *Q)
{
  double P[2][2], A[2][2], B[2][2];

  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      A[i][j] = *( N + ((2*i) + j) );
      B[i][j] = *( Q + ((2*i) + j) );
      P[i][j] = 0;      
    }
  }

  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      for (int k=0; k<2; k++)
      {
         P[i][j] += A[i][k]*B[k][i];  
      }
    }
  }  

  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      *( M + ((2*i) + j) ) = P[i][j];      
    }
  }
}

When I compile this using cython, I get the following error

mat.pyx:17:27: Cannot assign type 'double complex *' to 'double *'

I will be grateful to have some help here.

evamicur · Accepted Answer · 2018-02-27 17:36:52Z

1

This error message is telling you what's wrong:

mat.pyx:17:27: Cannot assign type 'double complex *' to 'double *'

That is, you have a double complex pointer from numpy (pointer to complex128 numpy dtype) and you're trying to pass that into the C++ function using double pointers. C++ needs to be able to deal with the complex numbers, so if you change your double* -> std::complex this should fix your problem

void mult(double *M, double *N, double *Q)

becomes

#include <complex>
void mult(std::complex<double> *M, std::complex<double> *N, std::complex<double> *Q)

Does numpy matrix multiply not suffice for your use case? Cython might be overkill.

Edit: Ok I finally got something, there's something a bit weird dealing with C++ std::complex and C double _Complex types.

cppmul.pyx:

import numpy as np
cimport numpy as np

cdef extern from "./matmult.h" nogil:
    void mult(np.complex128_t* M, np.complex128_t* N, np.complex128_t* Q)

def sim():
    cdef:
        np.ndarray[np.complex128_t,ndim=2] N = np.zeros(( 2 , 2 ), dtype=np.complex128)
        np.ndarray[np.complex128_t,ndim=2] Q = np.zeros(( 2 , 2 ), dtype=np.complex128)
        np.ndarray[np.complex128_t,ndim=2] M = np.zeros(( 2 , 2 ), dtype=np.complex128)

    N = np.array([[1.1 + 2j,2.2],[3.3,4.4]])
    Q = np.array([[3.3,4.4+5j],[5.5,6.6]])

    mult(&M[0,0], &N[0,0], &Q[0,0])
    print M

matmul.c:

#include "matmult.h"

void mult(complex_t *M, complex_t *N, complex_t *Q)
{
  complex_t P[2][2], A[2][2], B[2][2];

  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      A[i][j] = *( N + ((2*i) + j) );
      B[i][j] = *( Q + ((2*i) + j) );
      P[i][j] = 0;
    }
  }

  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      for (int k=0; k<2; k++)
      {
         P[i][j] += A[i][k]*B[k][i];
      }
    }
  }

  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      *( M + ((2*i) + j) ) = P[i][j];
    }
  }
}

matmult.h:

#include <complex.h>

typedef double _Complex complex_t;
void mult(complex_t *M, complex_t *N, complex_t *Q);

setup.py:

from distutils.core import setup
from Cython.Build import cythonize
from distutils.extension import Extension
import numpy as np

sourcefiles = ['cppmul.pyx', 'matmult.c']

extensions = [Extension("cppmul",
                        sourcefiles,
                        include_dirs=[np.get_include()],
                        extra_compile_args=['-O3']
)]

setup(
    ext_modules = cythonize(extensions)
)

after running python setup.py build_ext --inplace it imports and runs as expected

import cppmul
cppmul.sim()

result:

[[15.73 +6.6j 15.73 +6.6j]
 [43.56+16.5j 43.56+16.5j]]

edited Feb 27, 2018 at 17:36

answered Feb 27, 2018 at 8:26

evamicur

4033 silver badges8 bronze badges

Sign up to request clarification or add additional context in comments.

10 Comments

Commoner Over a year ago

evanmicur: Actually, I only gave the skeleton of my full code here. In my real work, I am working with matrices with hundreds of thousands of entries. Numpy matrix multiply is very slow for that, especially when the matrix entries are complex numbers. That's the reason why I wanted to use Cython here. It will be great for me if you can suggest a better way to approach the problem. Also, I tried to implement your suggestion (std::complex<double> *M). It may be close to solving the problem, but I am still getting this error in compilation: error: command 'gcc' failed with exit status 1

evamicur Over a year ago

Just to make sure: you're using the numpy.dot(A, b) or (A @ B in 3.4+) operation? Hand writing the loop will be slow. Numpy will be using a C/Fortran backend for the matrix multiply which I really don't think you have a good chance of beating by writing your own C code. But if you insist, that compile error is a bit hard to help with- error 1 just means it failed, so there's not enough detail to go here. Make sure you're consistent in types throughout the C++ code. You may want to consider a parallel matrix library like petsc4py or something if you really need a very fast matmul.

Commoner Over a year ago

evanmicur: This is the full error:

mat.pyx:6:18: Expected an identifier or literal building 'matmult_cy' extension gcc -fno-strict-aliasing -I/Users/anaconda2/include -arch x86_64 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/Users/anaconda2/lib/python2.7/site-packages/numpy/core/include -I/Users/anaconda2/include/python2.7 -c mat.cpp -o build/temp.macosx-10.6-x86_64-2.7/mat.o -O3 cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++   ^ error: command 'gcc' failed with exit status 1

Commoner Over a year ago

evanmicur: I have been using numpy(A,b). It is slow even after that. I have never used petsc4py. Do you know of a helpful resource from where I can learn how to that?

evamicur Over a year ago

Unfortunately the I've found the documentation less than helpful, but it may be on conda. It's a python wrapper for petsc, a c library. Gimme a minute I'll test some things and get back to you

|

user5928768 · Accepted Answer · 2018-02-27 07:49:17Z

0

try this

    #include "matmult.h"
 using namespace std;

int main(){}

void mult(double *M, double *N, double *Q)
{
 double P[2][2], A[2][2], B[2][2];

 for (int i=0; i<2; i++)
 {
for (int j=0; j<2; j++)
{
  A[i][j] = *( N + ((2*i) + j) );
  B[i][j] = *( Q + ((2*i) + j) );
  P[i][j] = 0;      
  }
  }

  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
    for (int k=0; k<2; k++)
    {
       P[i][j] += A[i][k]*B[k][i];  
     }
 }
}  

for (int i=0; i<2; i++)
 {
  for (int j=0; j<2; j++)
  {
  *(  ((2*i) + j) )+ M = P[i][j];      
  }
 }
}

answered Feb 27, 2018 at 7:49

user5928768

Collectives™ on Stack Overflow

Pass complex numpy array to C++ in Cython

2 Answers 2

10 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

10 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related