Hi Mathieu, Mathieu Blondel wrote:

Hello,

About one year ago, a high-level, objected-oriented SIMD API was added to Mono. For example, there is a class Vector4f for vectors of 4 floats and this class implements methods such as basic operators, bitwise operators, comparison operators, min, max, sqrt, shuffle directly using SIMD operations. You can have a look at the following pages for further details:

http://tirania.org/blog/archive/2008/Nov-03.html (blog post)

I am not sure how this could be applied to numpy case ? From what I can understand, this cannot be directly applied to python: the described changes are vm changes, and we cannot do anything at python vm level (I would guess the python vm to be too primitive to implement this kind of things anyway). I don't see how the high level API at the assembly level (Mono.Simd) would work either: the overhead of python and numpy to deal with 4 or 8 items in python would make this API useless from a speed POV. Implementing some numpy internal code in SIMD, and having a 'object oriented' C API for SIMD would indeed be nice - gcc provides SSE intrinsics, as well as visual studio (although the later seems quite buggy if I believe this link: http://www.virtualdub.org/blog/pivot/entry.php?id=162), which would make this in principle relatively easy. This is only my opinion (read other numpy dev may disagree), but I think that the numpy C code should be cleaned up before adding this kind of features: there is still too much coupling between the pure C core and the python machinery. Also, any use of SIMD code should be done at runtime IMHO (so that one binary can be used on multiple architectures), which has some issues on its own from a cross platform POV. David

The licenses look all hodge-podge:

[...]

However, if the good stuff is in the class libraries, that looks OK. But that still leaves it in C#, no?

I was mentioning Mono just to show that "this has been done" and also their API reference can serve as inspiration to design Numpy's own API. Mathieu

Mathieu Blondel wrote:

I saw the video of Peter Norvig at the last Scipy conference who was suggesting to merge Numpy into Cython. The SIMD API would be an argument in favor of this too because of the possible interactions between such a SIMD API and an array API.

Hm, I don't remember this - I guess I would have to look at the video. Do you know at which point of the presentation he discussed about SIMD ?

My original idea was to write the code in C with Intel/Alvitec/Neon intrinsics and have this code binded to be able to call it from Python. So the SIMD code would be compiled already, ready to be called from Python. Like you said, there's a risk that the overhead of calling Python is bigger than the benefit of using SIMD instructions. If it's worth trying out, an experiment can be made with Vector4f to see if it's even worth continuing with other types.

I am quite confident that the overhead will be way too significant for this approach to be useful. If you have two python objects, using + on it will induce at least one function call, and most likely several function calls at the python level. Python function calls are painfully slow (several thousand cycles per call in the most optimistic case). Python overhead is several order of magnitude bigger than what you can earn between SIMD and straightforward C. The only way I can see to make this work is to generate SIMD code from python (which would be a poor man's replacement for a JIT in a way), there was a presentation following this direction at scipy 09 conference.

I recently used SIMD instructions for a project and I realized that they cannot be activated in a standard Debian package, because the package has to remain general-purpose. So people who want to benefit the speed up have to compile my project from source...

Yes - that's unacceptable IMHO. The real solution is to include all the code at build time, detect at *runtime* which ISA is supported, and select the functions accordingly. The problem is that loading shared code at runtime in a cross platform way is complicated - python already does it, but unfortunately does not provide a C API for it AFAIK, so we would have to re-implement it in python. cheers, David

Wed, 21 Oct 2009 16:48:22 +0900, Mathieu Blondel wrote: [clip]

My original idea was to write the code in C with Intel/Alvitec/Neon intrinsics and have this code binded to be able to call it from Python. So the SIMD code would be compiled already, ready to be called from Python. Like you said, there's a risk that the overhead of calling Python is bigger than the benefit of using SIMD instructions. If it's worth trying out, an experiment can be made with Vector4f to see if it's even worth continuing with other types.

The overhead is quickly checked for multiplication with numpy arrays of varying size, without SSE: Overhead per iteration (ms): 1.6264549101 Time per array element (ms): 0.000936947636565 Cross-over point: 1735.90801303 #---------------------------------------------- import numpy as np from scipy import optimize import time import matplotlib.pyplot as plt def main(): data = [] for n in np.unique(np.logspace(0, 5, 20).astype(int)): print n m = 100 reps = 5 times = [] for rep in xrange(reps): x = np.zeros((n,), dtype=np.float_) start = time.time() #------------------ for k in xrange(m): x *= 1.1 #------------------ end = time.time() times.append(end - start) t = min(times) data.append((n, t)) data = np.array(data) def model(z): n, t = data.T overhead, per_elem = z return np.log10(t) - np.log10(overhead + per_elem * n) z, ier = optimize.leastsq(model, [1., 1.]) overhead, per_elem = z print "" print "Overhead per iteration (ms):", overhead*1e3 print "Time per array element (ms):", per_elem*1e3 print "Cross-over point: ", overhead/per_elem n = np.logspace(0, 5, 500) plt.loglog(data[:,0], data[:,0]/data[:,1], 'x', label=r'measured') plt.loglog(n, n/(overhead + per_elem*n), 'k-', label=r'fit to $t = a + b n$') plt.xlabel(r'$n$') plt.ylabel(r'ops/second') plt.grid(1) plt.legend() plt.show() if __name__ == "__main__": main()

On Wed, Oct 21, 2009 at 5:23 PM, David Cournapeau wrote:

Ah, I think you are mistaken, then - he referred to merging numpy and scipy into python during his talk, not cython.

Oh, I meant to say CPython (the default implementation of Python), not Cython. I didn't realize that they were different projects. So the method dispatch seems to be a great obstacle to an object-oriented SIMD API. That would seem more feasible in C++ with non-virtual methods. Java has final methods, which can be useful information to the JIT. C# seems to have "sealed" methods. Interestingly, the Mono.SIMD API uses static methods, which I guess is to avoid the dispatch problem. But it makes the code look uglier. For example, instead of a + b, you have to do Vector4f.Addition(a, b). Mathieu

On Wed, Oct 21, 2009 at 6:12 PM, Francesc Alted wrote:

A Wednesday 21 October 2009 07:44:39 Mathieu Blondel escrigué:

...

Hello,

About one year ago, a high-level, objected-oriented SIMD API was added to Mono. For example, there is a class Vector4f for vectors of 4 floats and this class implements methods such as basic operators, bitwise operators, comparison operators, min, max, sqrt, shuffle directly using SIMD operations. [clip]

It is important to stress out that all the above operations, except probably sqrt, are all memory-bound operations, and that implementing them for numpy would not represent a significant improvement at all.

This is because numpy is a package that works mainly with arrays in an element-wise way, and in this scenario, the time to transmit data to CPU dominates, by and large, over the time to perform operations.

Is it general, or just for simple operations in numpy and ufunc ? I remember that for music softwares, SIMD used to matter a lot, even for simple bus mixing (which is basically a ax+by with a, b scalars and x y the input arrays). Do you have any interest in adding SIMD to some core numpy (transcendental functions). If so, I would try to go back to the problem of runtime SSE detection and loading of optimized shared library in a cross-platform way - that's something which should be done at some point in numpy, and people requiring it would be a good incentive. David

A Wednesday 21 October 2009 14:27:46 David Cournapeau escrigué:

...

This is because numpy is a package that works mainly with arrays in an element-wise way, and in this scenario, the time to transmit data to CPU dominates, by and large, over the time to perform operations.

Is it general, or just for simple operations in numpy and ufunc ? I remember that for music softwares, SIMD used to matter a lot, even for simple bus mixing (which is basically a ax+by with a, b scalars and x y the input arrays).

This is general, as long as the dataset has to be brought from memory to CPU, and operations to be done are element-wise and simple (i.e. not transcendental). SIMD does matter in general when the dataset: 1) is already in cache 2) you have to perform costly operations (mainly transcendental) 3) a combination of the above I don't know the case for music software, but if you say that ax+by are accelerated by SIMD, I'd say that case 1) is happening.

Do you have any interest in adding SIMD to some core numpy (transcendental functions). If so, I would try to go back to the problem of runtime SSE detection and loading of optimized shared library in a cross-platform way - that's something which should be done at some point in numpy, and people requiring it would be a good incentive.

I don't personally have a lot of interest implementing this for numpy. But in case anyone does, I find the next library: http://gruntthepeon.free.fr/ssemath/ very interesting. Perhaps there could be other (free) implementations... -- Francesc Alted

On Wed, Oct 21, 2009 at 2:14 PM, Pauli Virtanen

wrote:

Wed, 21 Oct 2009 14:47:02 +0200, Francesc Alted wrote: [clip]

...

...

Do you have any interest in adding SIMD to some core numpy (transcendental functions). If so, I would try to go back to the problem of runtime SSE detection and loading of optimized shared library in a cross-platform way - that's something which should be done at some point in numpy, and people requiring it would be a good incentive.

I don't personally have a lot of interest implementing this for numpy. But in case anyone does, I find the next library:

http://gruntthepeon.free.fr/ssemath/

very interesting. Perhaps there could be other (free) implementations...

Optimized transcendental functions could be interesting. For example for tanh, call overhead is overcome already for ~30-element arrays.

Since these are ufuncs, I suppose the SSE implementations could just be put in a separate module, which is always compiled. Before importing the module, we could simply check from Python side that the CPU supports the necessary instructions. If everything is OK, the accelerated implementations would then just replace the Numpy routines.

This type of project could probably also be started outside Numpy, and just monkey-patch the Numpy routines on import.

-- Pauli Virtanen

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

Anyone seen the corepy numpy gsoc project? http://numcorepy.blogspot.com/ It implements a number of functions with the corepy runtime assembler. The project showed nice simd speedups for numpy. I've been following the liborc project... which is a runtime assembler that uses a generic assembly language and supports many different simd assembly languages (eg SSE, MMX, ARM, Altivec). It's the replacement for the liboil library (used in gstreamer etc). http://code.entropywave.com/projects/orc/ cu!

On Wed, Oct 21, 2009 at 11:38 AM, Gregor Thalhammer wrote:

I once wrote a module that replaces the built in transcendental functions of numpy by optimized versions from Intels vector math library. If someone is interested, I can publish it. In my experience it was of little use since real world problems are limited by memory bandwidth. Therefore extending numexpr with optimized transcendental functions was the better solution. Afterwards I discovered that I could have saved the effort of the first approach since gcc is able to use optimized functions from Intels vector math library or AMD's math core library, see the doc's of -mveclibabi. You just need to recompile numpy with proper compiler arguments.

Do you have a link to the documentation for -mveclibabi? I can't find this anywhere and I'm *very* interested. Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma

...

I once wrote a module that replaces the built in transcendental functions of numpy by optimized versions from Intels vector math library. If someone is interested, I can publish it. In my experience it was of little use since real world problems are limited by memory bandwidth. Therefore extending numexpr with optimized transcendental functions was the better solution. Afterwards I discovered that I could have saved the effort of the first approach since gcc is able to use optimized functions from Intels vector math library or AMD's math core library, see the doc's of -mveclibabi. You just need to recompile numpy with proper compiler arguments.

I'm interested. I'd like to try AMD rather than intel, because AMD is easier to obtain. I'm running on intel machine, I hope that doesn't matter too much. What exactly do I need to do? I see that numpy/site.cfg has an MKL section. I'm assuming I should not touch that, but just mess with gcc flags?

sigh; yet another email dropped by the list. David Warde-Farley wrote:

On 21-Oct-09, at 9:14 AM, Pauli Virtanen wrote:

...

Since these are ufuncs, I suppose the SSE implementations could just be put in a separate module, which is always compiled. Before importing the module, we could simply check from Python side that the CPU supports the necessary instructions. If everything is OK, the accelerated implementations would then just replace the Numpy routines.

Am I mistaken or wasn't that sort of the goal of Andrew Friedley's CorePy work this summer?

Looking at his slides again, the speedups are rather impressive. I wonder if these could be usefully integrated into numpy itself?

Yes, my GSoC project is closely related, though I didn't do the CPU detection part, that'd be easy to do. Also I wrote my code specifically for 64-bit x86. I didn't focus so much on the transcendental functions, though they wouldn't be too hard to implement. There's also the possibility to provide implementations with differing tradeoffs between accuracy and performance. I think the blog link got posted already, but here's relevant info: http://numcorepy.blogspot.com http://www.corepy.org/wiki/index.php?title=CoreFunc I talked about this in my SciPy talk and up-coming paper, as well. Also people have just been talking about x86 in this thread -- other architectures could be supported too; eg PPC/Altivec or even Cell SPU and other accelerators. I actually wrote a quick/dirty implementation of addition and vector normalization ufuncs for Cell SPU recently. Basic result is that overall performance is very roughly comparable to a similar speed x86 chip, but this is a huge win over just running on the extremely slow Cell PPC cores. Andrew

Mathieu Blondel skrev:

Hello,

About one year ago, a high-level, objected-oriented SIMD API was added to Mono. For example, there is a class Vector4f for vectors of 4 floats and this class implements methods such as basic operators, bitwise operators, comparison operators, min, max, sqrt, shuffle directly using SIMD operations. I think you are confusing SIMD with Intel's MMX/SSE instruction set.

SIMD means "single instruction - multiple data". NumPy is interherently an object-oriented SIMD API: array1[:] = array2 + array3 is a SIMD instruction by definition. SIMD instructions in hardware for length-4 vectors are mostly useful for 3D graphics. But they are not used a lot for that purpose, because GPUs are getting common. SSE is mostly for rendering 3D graphics without a GPU. There is nothing that prevents NumPy from having a Vector4f dtype, that internally stores four float32 and is aligned at 16 byte boundaries. But it would not be faster than the current float32 dtype. Do you know why? The reason is that memory access is slow, and computation is fast. Modern CPUs are starved. The speed of NumPy is not limited by not using MMX/SSE whenever possible. It is limited from having to create and delete temporary arrays all the time. You are suggesting to optimize in the wrong place. There is a lot that can be done to speed up computation: There are optimized BLAS libraries like ATLAS and MKL. NumPy uses BLAS for things like matrix multiplication. There are OpenMP for better performance on multicores. There are OpenCL and CUDA for moving computation from CPUs to GPU. But the main boost you get from going from NumPy to hand-written C or Fortran comes from reduced memory use.

existing discussion here. Memory-alignment is an import related issue since non-aligned movs can tank the performance.

You can align an ndarray on 16-byte boundary like this: def aligned_array(N, dtype): d = dtype() tmp = numpy.zeros(N * d.nbytes + 16, dtype=numpy.uint8) address = tmp.__array_interface__['data'][0] offset = (16 - address % 16) % 16 return tmp[offset:offset+N].view(dtype=dtype) Sturla Molden

On Wed, Oct 21, 2009 at 22:32, Mathieu Blondel wrote:

On Thu, Oct 22, 2009 at 11:31 AM, Sturla Molden wrote:

...

Mathieu Blondel skrev:

...

Hello,

About one year ago, a high-level, objected-oriented SIMD API was added to Mono. For example, there is a class Vector4f for vectors of 4 floats and this class implements methods such as basic operators, bitwise operators, comparison operators, min, max, sqrt, shuffle directly using SIMD operations. I think you are confusing SIMD with Intel's MMX/SSE instruction set.

OK, I should have said "Object-oriented SIMD API that is implemented using hardware SIMD instructions".

No, I think you're right. Using "SIMD" to refer to numpy-like operations is an abuse of the term not supported by any outside community that I am aware of. Everyone else uses "SIMD" to describe hardware instructions, not the application of a single syntactical element of a high level language to a non-trivial data structure containing lots of atomic data elements. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco

...

OK, I should have said "Object-oriented SIMD API that is implemented using hardware SIMD instructions".

No, I think you're right. Using "SIMD" to refer to numpy-like operations is an abuse of the term not supported by any outside community that I am aware of. Everyone else uses "SIMD" to describe hardware instructions, not the application of a single syntactical element of a high level language to a non-trivial data structure containing lots of atomic data elements.

I agree with Sturla, for instance nVidia GPUs do SIMD computations with blocs of 16 values at a time, but the hardware behind can't compute on so much data at a time. It's SIMD from our point of view, just like Numpy does ;) Matthieu -- Information System Engineer, Ph.D. Website: http://matthieu-brucher.developpez.com/ Blogs: http://matt.eifelle.com and http://blog.developpez.com/?blog=92 LinkedIn: http://www.linkedin.com/in/matthieubrucher

Mathieu Blondel skrev:

Peter Norvig suggested to merge Numpy into Cython but he didn't mention SIMD as the reason (this one is from me).

I don't know what Norvig said or meant. However: There is NumPy support in Cython. Cython has a general syntax applicable to any PEP 3118 buffer. (As NumPy is not yet PEP 3118 compliant, NumPy arrays are converted to Py_buffer structs behind the scenes.) Support for optimized vector expressions might be added later. Currently, slicing works as with NumPy in Python, producing slice objects and invoking NumPy's own code, instead of being converted to fast inlined C. The PEP 3118 buffer syntax in Cython can be used to port NumPy to Py3k, replacing the current C source. That might be what Norvig meant if he suggested merging NumPy into Cython. S.M.

Mathieu Blondel skrev:

As I wrote earlier in this thread, I confused Cython and CPython. PN was suggesting to include Numpy in the CPython distribution (not Cython). The reason why was also given earlier.

First, that would currently not be possible, as NumPy does not support Py3k. Second, the easiest way to port NumPy to Py3k is Cython, which would prevent adoption in the Python standard library. At least they have to change their current policy. Also with NumPy in the standard library, any modification to NumPy would require a PEP. But Python should have a PEP 3118 compliant buffer object in the standard library, which NumPy could subclass. S.M.

Robert Kern wrote:

On Wed, Oct 21, 2009 at 22:32, Mathieu Blondel wrote:

...

On Thu, Oct 22, 2009 at 11:31 AM, Sturla Molden wrote:

...

Mathieu Blondel skrev:

...

Hello,

About one year ago, a high-level, objected-oriented SIMD API was added to Mono. For example, there is a class Vector4f for vectors of 4 floats and this class implements methods such as basic operators, bitwise operators, comparison operators, min, max, sqrt, shuffle directly using SIMD operations.

I think you are confusing SIMD with Intel's MMX/SSE instruction set.

OK, I should have said "Object-oriented SIMD API that is implemented using hardware SIMD instructions".

No, I think you're right. Using "SIMD" to refer to numpy-like operations is an abuse of the term not supported by any outside community that I am aware of. Everyone else uses "SIMD" to describe hardware instructions, not the application of a single syntactical element of a high level language to a non-trivial data structure containing lots of atomic data elements.

BTW, is there any term for this latter concept that's not SIMD or "vector operation"? It would be good to have a word to distinguish this concept from both CPU instructions and linear algebra. (Personally I think describing NumPy as SIMD and use "SSE/MMX" for CPU instructions makes best sense, but I'm happy to yield to conventions...) Dag Sverre

On Thu, Oct 22, 2009 at 02:35, Sturla Molden wrote:

Robert Kern skrev:

...

No, I think you're right. Using "SIMD" to refer to numpy-like operations is an abuse of the term not supported by any outside community that I am aware of. Everyone else uses "SIMD" to describe hardware instructions, not the application of a single syntactical element of a high level language to a non-trivial data structure containing lots of atomic data elements.

Then you should pick up a book on parallel computing.

I would be delighted to see a reference to one that refers to a high level language's API as SIMD. Please point one out to me. It's certainly not any of the ones I have available to me.

It is common to differentiate between four classes of computers: SISD, MISD, SIMD, and MIMD machines.

A SISD system is the classical von Neuman machine. A MISD system is a pipelined von Neuman machine, for example the x86 processor.

A SIMD system is one that has one CPU dedicated to control, and a large collection of subordinate ALUs for computation. Each ALU has a small amount of private memory. The IBM Cell processor is the typical SIMD machine.

A special class of SIMD machines are the so-called "vector machines", of which the most famous is the Cray C90. The MMX and SSE instructions in Intel Pentium processors are an example of vector instructions. Some computer scientists regard vector machines a subtype of MISD systems, orthogonal to piplines, because there are no subordinate ALUs with private memory.

MIMD systems multiple independent CPUs. MIMD systems comes in two categories: shared-memory processors (SMP) and distributed-memory machines (also called cluster computers). The dual- and quad-core x86 processors are shared-memory MIMD machines.

Many people associate the word SIMD with SSE due to Intel marketing. But to the extent that vector machines are MISD orthogonal to piplined von Neuman machines, SSE cannot be called SIMD.

That's a fair point, but unrelated to whether or not numpy can be labeled SIMD. These all refer to hardware.

NumPy is a software simulated vector machine, usually executed on MISD hardware. To the extent that vector machines (such as SSE and C90) are SIMD, we must call NumPy an object-oriented SIMD library.

numpy does not "simulate" anything. It is an object-oriented library. If numpy could be said to "simulate" a vector machine, than just about any object-oriented library that overloads operators could. It creates a false equivalence between numpy and software that actually does simulate hardware. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco

Robert Kern skrev:

I would be delighted to see a reference to one that refers to a high level language's API as SIMD. Please point one out to me. It's certainly not any of the ones I have available to me.

Numerical Receipes in Fortran 90, page 964 and 985-986, describes the syntax of Fortran 90 and 95 as SIMD. Peter Pacheco's book on MPI describes the difference between von Neumann machines and vector machines as analogous to the difference between Fortran77 and Fortran 90 (with an example from Fortran90 array slicing). He is ambigous as to whether vector machines really are SIMD, or more related to pipelined von Neumann machines. Grama et al. "Introduction to Parallel Computing" describes SIMD as an "architecture", but it is more or less clear that the mean hardware. They do say the Fortran 90 "where statement" is a primitive used to support selective execution on SIMD processors, as conditional execution (if statements) are detrimental to performance. So at least we here have three books claiming that Fortran is a language with special primities for SIMD processors.

That's a fair point, but unrelated to whether or not numpy can be labeled SIMD. These all refer to hardware.

Actually I don't think the distinction is that important as we are taking about Turing machines. Also, a lot of what we call "hardware" is actually implemented as software on the chip: The most extreme example would be Transmeta, which completely software emulated x86 processors. The vague distinction between hardware and software is why we get patents on software in Europe, although pure software patents are prohibited. One can always argue that the program and the computer together constitutes a physical device; and circumventing patents by moving hardware into software should not be allowed. The distinction between hardware and software is not as clear as programmers tend to believe. Another thing is that performance issues for vector machines and "vector languages" (Fortran 90, Matlab, NumPy) are similar. Precisely the same situations that makes NumPy and Matlab code slow are detrimental on SIMD/vector hardware. That would for example be long for loops with conditional if statements. On the other hand, vectorized operations over arrays, possibly using where/find masks, are fast. So although NumPy is not executed on a vector machine like the Cray C90, it certainly behaves like one performance wise. I'd say that a MIMD machine running NumPy is a Turing machine emulating a SIMD/vector machine. And now I am done with this stupid discussion... Sturla Molden