<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
  <head>
    <meta content="text/html; charset=ISO-8859-1"
      http-equiv="Content-Type">
  </head>
  <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>
         char *addr = arr->data;<br>
         for (int i=0; i<arr->dimensions[0]; i++) {<br>
             y = (double *)addr;<br>
             acc += y*y;<br>
             addr +=  arr->strides[0];<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>