[Numpy-discussion] Cython numerical syntax revisited

[Sturla Molden]

On 3/5/2009 10:11 AM, Dag Sverre Seljebotn wrote:


> Cython can relatively easily transform things like
> 
> cdef int[:,:] a = ..., b = ...
> c = a + b * b

Now you are wandering far into Fortran territory...


> If a and b are declared as contiguous arrays and "restrict", I suppose 
> the C compiler could do the most efficient thing in a lot of cases? 
> (I.e. "cdef restrict int[:,:,"c"]" or similar)

A Fortran compiler can compile a vectorized expression like

   a = b*c(i,:) + sin(k)

into

   do j=1,n
     a(j) = b(j)*c(i,j) + sin(k(j))
   end do

The compiler do this because Fortran has strict rules prohibiting 
aliasing, and because the instrinsic function sin is declared 
'elemental'. On the other hand, if the expression contains functions not 
declared 'elemental' or 'pure', there may be side effects, and temporary 
copies must be made. The same could happen if the expression contained 
variables declared 'pointer', in which case it could contain aliases.

I think in the case of numexpr, it assumes that NumPy ufuncs are 
elemental like Fortran intrinsics.

Matlab's JIT compiler works with the assumption that arrays are 
inherently immutable (everything has copy-on-write semantics). That 
makes life easier.


> However if one has a strided array, numexpr could still give an 
> advantage over such a loop. Or?

Fortran compilers often makes copies of strided arrays. The trick is to 
make sure the working arrays fit in cache. Again, this is safe when the 
expression only contains 'elemental' or 'pure' functions. Fortran also 
often does "copy-in copy-out" if a function is called with a strided 
array as argument.

S.M.