Ban the Asterisk!

QUESTION: I'm trying to calculate shear and I have to do it a lot. Does anyone know of a way to execute the following code as fast as possible? In the code, delta is an integer array.

    FOR i=0,n-1 DO array[*,i] = SHIFT( array[*,i], delta[i] )

ANSWER: The following discussion occurred between Wayne Landsman, Martin Downing, and Ken Bowman. We hear first from Wayne.

Well, I think a loop over rows is necessary, but the speed can be vastly improved by using the folllowing IDL maxim for efficient programming, which gets my vote in the IDL Maxim Most Deserving of Wider Recognition category.

***Avoid the use of an asterisk on the left-hand side of an assignment statement.

In this case, it means rewriting the assignment here as:

    FOR i=0,n-1 DO array[0,i] = SHIFT( array[*,i], delta[i] )

For a (2048,2048) floating point array on my Solaris machine running V5.5, I find a factor of 15 improvement in speed.

I believe that if one uses an asterisk on the left hand side, IDL first creates a temporary variable, places the contents of the right hand side into this temporary variable, and then places the temporary variable back into array. By specifying array[0,i] one directly gives the starting location for placing the contents, which greatly increases the speed of the operation.

The original questioner was impressed with this answer. He had this to say.

WOW! Wayne, that's magic. I get 10-15X speed-up. The world has to be told about this. One drawback - the code does not strictly make sense... But, hell, who cares!?

Martin seconded his surprise.

I second that WOW! I hadn't realised specifying a single cell of an image on the left hand side allows intelligent access to the elements mapped by the dimensional shape of the right hand side. Offset assignments in 2D or 3D images become a breeze.

For example, to set image[100:109,100:129] = 255 all you have to do it this.

    image[100,100] = Replicate(255, 10, 30)

I had stupidly assumed that the assignment would follow C language rules and fill in a raster ordering. Thanks for the tip. I will be using this method from now on.

Before we could all get too giddy, though, Ken threw a bit of reality on the fire.

I thought using an asterisk on the LHS meant that IDL created a temporary index array, something like this:

   FOR i = 0, n-1 DO array[[0,1,2,3,...,m-1],i] = SHIFT(array[*,i],delta[i])

where m is the size of the first dimension. The slow down is due to the the indirect array indexing.

The fast method (e.g., array[0,i] = ...) specifies where to start storing the RHS, so it only works for the first dimension (that is, in memory order).

Google
 
Web Coyote's Guide to IDL Programming