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.