[SciPy-dev] Re: [SciPy-user] Low-level code integrationdiscussion at scipy'03?

[Fernando Perez]

Perry Greenfield wrote:

> So given that, I (with my admittedly limited imagination :-) can't see
> how syntax such as
> 
> A(i,j,k,l,m,n)
> 
> can work in C (I hope someone can set me straight). Perhaps something
> like
> 
> Index6d(A,i,j,k,l,m,n)
> 
> it the best that can be done. 

Agreed.  In pure C, I don't see any other option either.

> (A[i][j]...[n] would work but would need
> arrays of pointers to pointer to... and I don't think we want to go
> there.)

Yup.  I've done that for 2d arrays, writing a simple function which takes a 
Numpy array and returns the necessary scaffolding of pointers for A[i][j] to 
work.  But it's ugly, wasteful, and doesn't scale cleanly to higher dimensions.

> On the other hand, syntax like that proposed can be provided in C++.
> Still, I'd be very surprised if the blitz people haven't already given
> a great deal of thought to this. They probably can explain why it is
> slower and why it is difficult to do much better than this.

Well, their stuff works really well for pure indexing.  The overhead I 
mentioned was for the magic

U = sum(M(i,j),T(k,l,m,n,j),j);

expression, which is setting all U(i,k,l,m,n) for you by building all the 
nested loops on the fly.  If you directly index the arrays, I think blitz 
performance is indeed very good (at least the benchmarks they provide claim so).

I hope Julian Cummings can make it to the discussion, so we can ask him about 
this directly.

Cheers,

f.