I am trying to create a function that can either beat numexpr or perform comparably for the vectorized mathematical operation k * exp(ee), where k and ee are 1D arrays of complex numbers. I had tried using OpenMP in Cython directly, but I have resorted to creating a DLL in C++ and then wrapping with Cython. Here are the components (of current interest is the function multi_mlt_exp_c4).I am using Windows 10, 64 bit, Intel processor, and Python 3.x.C++ DLL, Header:#define EXTERN_DLL_EXPORT extern C __declspec(dllexport)#include <complex>#include <stdlib.h>EXTERN_DLL_EXPORT void cexp_scale_c4(std::complex<float> k, std::complex<float> ee,                                     std::complex<float>*res);EXTERN_DLL_EXPORT void multi_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee,                                        int sz, std::complex<float>* res, int threads);EXTERN_DLL_EXPORT void multi_mlt_exp_c8(std::complex<double>* k, std::complex<double>* ee,                                        int sz, std::complex<double>* res, int threads);EXTERN_DLL_EXPORT void serial_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee,                                         int sz, std::complex<float>* res);C++ DLL, Source#include stdafx.h#include <stdio.h>#include <omp.h>#include C_Fun_Lib.hvoid cexp_scale_c4(std::complex<float> k, std::complex<float> ee, std::complex<float>*res){    *res = k*exp(ee);}void multi_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee, int sz, std::complex<float>* res, int threads){    #pragma omp parallel num_threads(threads)    {        #pragma omp for        for (int i = 0; i < sz; i++) res[i] = k[i] * exp(ee[i]);    }}void multi_mlt_exp_c8(std::complex<double>* k, std::complex<double>* ee, int sz, std::complex<double>* res, int threads){    #pragma omp parallel num_threads(threads)    {        #pragma omp for        for (int i = 0; i < sz; i++) res[i] = k[i] * exp(ee[i]);    }}void serial_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee, int sz, std::complex<float>* res){    for (int i = 0; i < sz; i++) res[i] = k[i] * exp(ee[i]);}Here are the relevant Python / Cython components:PXDcdef extern from C_Fun_Lib.h:    void multi_mlt_exp_c4(float complex* k, float complex* ee, int sz,                          float complex* res, int threads);    void serial_mlt_exp_c4(float complex* k, float complex* ee, int sz,                           float complex* res);    void multi_mlt_exp_c8(double complex* k, double complex* ee, int sz,                          double complex* res, int threads);PYXcimport csampleimport numpy as npcimport numpy as np# Import some functionality from Python and the C stdlibfrom cpython.pycapsule cimport *from libc.stdlib cimport malloc, freedef mlt_exp_c4(np.ndarray[np.complex64_t, ndim=1] k,               np.ndarray[np.complex64_t, ndim=1] ee, int sz,               np.ndarray[np.complex64_t, ndim=1] res, int threads):    csample.multi_mlt_exp_c4(&k[0], &ee[0], sz, &res[0], threads)def mlt_exp_c4_serial(np.ndarray[np.complex64_t, ndim=1] k,               np.ndarray[np.complex64_t, ndim=1] ee, int sz,               np.ndarray[np.complex64_t, ndim=1] res):    csample.serial_mlt_exp_c4(&k[0], &ee[0], sz, &res[0])def mlt_exp_c8(np.ndarray[np.complex128_t, ndim=1] k,               np.ndarray[np.complex128_t, ndim=1] ee, int sz,               np.ndarray[np.complex128_t, ndim=1] res, int threads):    csample.multi_mlt_exp_c8(&k[0], &ee[0], sz, &res[0], threads)Setup.pyfrom distutils.core import setupfrom distutils.extension import Extensionfrom Cython.Distutils import build_extimport numpy as npext_modules = [    Extension('sample',              ['sample.pyx'],              language=c++,              libraries=['C_Fun_Lib'],              library_dirs=['.'])    ]setup(    name = 'sample',    cmdclass = {'build_ext': build_ext},    ext_modules = ext_modules,    include_dirs=[np.get_include()])And then finally here is my timing code and results. I am specifically trying to beat Numexpr because, currently, that application up-casts complex64 arrays to complex128 arrays, and I am trying to save memory. Writing my own code gives me control over the memory, but I can't quite match the speed of Numexpr.import sampleimport numpy as npimport numexpr as neimport timek = np.ones(int(2.5e8), dtype='complex64')*1.1234 + 1.1234jee = ksz = k.szres = np.zeros(int(2.5e8), dtype='complex64')  # Resultstic = time.time()sample.mlt_exp_c4(k, ee, sz, res, 8)toc = time.time()print('Elasped MP: %f' % (toc - tic))# Reset resultsres = np.zeros(int(2.5e8), dtype='complex64')tic = time.time()res = ne.evaluate('k*exp(ee)')toc = time.time()print('Elasped NE: %f' % (toc - tic))>>>>>>>Elasped MP: 2.292317Elasped NE: 1.184981Things I have triedFirst, out of all the things I've tried, this is my current best. I have experimented to see if using complex128s makes a difference (since I'm on a 64 bit machine). It does not. I have tried doing the complex exponential and multiplication directly and using Euler's formula instead of using the C++ built in complex type and functions... that is slower. I have tried using Cython code directly with OpenMP support... that is slower.It may very well be that trying to match the highly optimized Numexpr code is a fool's errand. I just wanted you all to review this code to see if there were any obvious speed ups that I am missing.
Calculating complex exponential numbers using OpenMP in Cython
python;c++;numpy;openmp;native code