<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
<body text="#000000" bgcolor="#ffffff">
Den 13.02.2011 04:30, skrev Travis Oliphant:
<blockquote
cite="mid:E4B697CE-1621-45D4-97DF-F426483BF083@enthought.com"
type="cite">
<div><br>
</div>
<div>One of the advantages of an open source community is that
different needs get addressed by the people who need them.   So,
if you really need a faster sum and are willing to do the hard
work to make it happen (or convince someone else to do it for
you), then there will be many to say thank you when you are
finished.    You should know though that in my experience It is
harder than you might think at first to "convince someone else
to do it for you."<br>
</div>
</blockquote>
<br>
For things like sum (and may other ufuncs), one could just use
Fortran instrincs and let the Fortran compiler do the rest.
(Obviously NumPy should not depend on a Fortran compiler, unlike
SciPy, but that is another matter.) <br>
<br>
Assuming a wrapper takes care of C and Fortran ordering, we can just
f2py this:<br>
<br>
subroutine fortran_sum_1d( m, arr, out )<br>
use, intrinsic :: sum<br>
integer :: m<br>
real :: arr(m), out<br>
out = sum(arr)<br>
end subroutine<br>
<br>
subroutine fortran_sum_2d( m, n, arr, k, out, dim )<br>
use, intrinsic :: sum<br>
integer :: m, n, k, dim<br>
real :: arr(m,n), out(k)<br>
out = sum(arr, dim+1) ! fortran dims start at 1<br>
end subroutine<br>
<br>
Similar things can be done for other ufuncs. We can also use
Fortran's "elemental" keyword to create our own. The Fortran
compilers from Absoft and Intel (and possibly Portland) can even
take code like this and make it efficient for multicore CPUs. (I
don't know if GNU gfortran can do this too, but I don't expect it
too.)<br>
<br>
One limitation is that the Fortran compiler must know the number of
dimensions, while NumPy arrays are flexible in this respect. Another
limitation is strides, but we can deal with them if they are a
multiple of the dtype size, but not arbitrary number of bytes like
NumPy.<br>
<br>
subroutine fortran_sum_1d_with_strides( m, s, arr, out )<br>
use, intrinsic :: sum<br>
integer :: m, s<br>
real, target :: arr(m)<br>
real, pointer, dimension(:) :: parr<br>
real :: out<br>
parr => arr(::s)<br>
out = sum(parr)<br>
end subroutine<br>
<br>
I think for most people, Fortran is much easier than using NumPy's C
API. Creating fast "ufunc" like functions is very easy with Fortran
95. <br>
<br>
If we e.g. want the L2 norm of a vector, it looks like this in
Fortran:<br>
<br>
subroutine l2norm( m, arr, out )<br>
use, intrinsic :: sum, sqrt<br>
integer :: m<br>
real :: arr(m), out<br>
out = sqrt(sum(arr**2)) ! looks quite like Python...<br>
end subroutine<br>
<br>
Just notice the messyness of C in comparison (there are much worse
examples):<br>
<br>
#include <Python.h><br>
#include <numpy/arrayobject.h<em></em>><br>
#include <math.h><br>
<br>
double l2norm(PyArrayObject *arr)<br>
{<br>
double acc = 0.0, *y;<br>
for (int i=0; i<arr->dimensions; i++) {<br>
acc += y*y;<br>
}<br>
return sqrt(acc);<br>
}<br>
<br>
<br>
<br>
<br>
Sturla<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
</body>
</html>