<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Jan 8, 2014 at 2:39 PM, Julian Taylor <span dir="ltr"><<a href="mailto:jtaylor.debian@googlemail.com" target="_blank">jtaylor.debian@googlemail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi,<br>
Since AMDs bulldozer and Intels Haswell x86 cpus now also support the<br>
fused-multiply-and-add operation in hardware.<br>
<br>
<a href="http://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation" target="_blank">http://en.wikipedia.org/wiki/Multiply–accumulate_operation</a><br>
<br>
This operation is interesting for numpy for two reasons:<br>
 - Only one rounding is involved in two arithmetic operations, this is<br>
good reducing rounding errors.<br>
 - Two operations are done in one memory pass, so it improves the<br>
performance if ones operations are bound by the memory bandwidth which<br>
is very common in numpy.<br>
<br>
I have done some experiments using a custom ufunc:<br>
<a href="https://github.com/juliantaylor/npufuncs" target="_blank">https://github.com/juliantaylor/npufuncs</a><br>
<br>
See the README.md on how to try it out. It requires a recent GCC<br>
compiler, at least 4.7 I think.<br>
<br>
It contains SSE2, AVX, FMA3 (AVX2), FMA4 and software emulation<br>
variants. Edit the file to select which one to use. Note if the cpu does<br>
not support the instruction it will just crash.<br>
Only the latter three are real FMA operations, the SSE2 and AVX variants<br>
just perform two regular operations in one loop.<br>
<br>
My current machine only supports SSE2, so here are the timings for it:<br>
<br>
In [25]: a = np.arange(500000.)<br>
In [26]: b = np.arange(500000.)<br>
In [27]: c = np.arange(500000.)<br>
<br>
In [28]: %timeit npfma.fma(a, b, c)<br>
100 loops, best of 3: 2.49 ms per loop<br>
<br>
In [30]: def pure_numpy_fma(a,b,c):<br>
   ....:     return a * b + c<br>
<br>
In [31]: %timeit pure_numpy_fma(a, b, c)<br>
100 loops, best of 3: 7.36 ms per loop<br>
<br>
<br>
In [32]: def pure_numpy_fma2(a,b,c):<br>
   ....:     tmp = a *b<br>
   ....:     tmp += c<br>
   ....:     return tmp<br>
<br>
In [33]: %timeit pure_numpy_fma2(a, b, c)<br>
100 loops, best of 3: 3.47 ms per loop<br>
<br>
<br>
As you can see even without real hardware support it is about 30% faster<br>
than inplace unblocked numpy due better use of memory bandwidth. Its<br>
even more than two times faster than unoptimized numpy.<br>
<br>
If you have a machine capable of fma instructions give it a spin to see<br>
if you get similar or better results. Please verify the assembly<br>
(objdump -d fma-<suffix>.o) to check if the compiler properly used the<br>
machine fma.<br>
<br>
An issue is software emulation of real fma. This can be enabled in the<br>
test ufunc with npfma.set_type("libc").<br>
This is unfortunately incredibly slow about a factor 300 on my machine<br>
without hardware fma.<br>
This means we either have a function that is fast on some platforms and<br>
slow on others but always gives the same result or we have a fast<br>
function that gives better results on some platforms.<br>
Given that we are not worth that what numpy currently provides I favor<br>
the latter.<br>
<br>
Any opinions on whether this should go into numpy or maybe stay a third<br>
party ufunc?<br></blockquote><div><br></div><div>Multiply and add is a standard function that I think would be good to have in numpy. Not only does it save on memory accesses, it saves on temporary arrays.<br><br></div><div>
Another function that could be useful is a |a|**2 function, abs2 perhaps.<br><br></div><div>Chuck <br></div></div></div></div>