Speeding up / Vectorization of Octave code

There was a really interesting posts on the Help-Octave lists yesterday (25^th) June where
a vectorization paradigm was elucidated. Here it is for the non-Octavers...

Paul Kienzle to Rob Teasdale
date Jun 25, 2007 9:34 PM
subject Re: Request help speeding up script


On Jun 25, 2007, at 5:23 PM, Rob Teasdale wrote:
> Hi all,
>
> I am still learning Octave, and I am still discovering new ways of
> optimising code. I have managed to vectorise a lot of my 'for' loops,
> and I have managed to speed my small script up dramatically, although
> I am having trouble with this last code snippet below. What I am
> trying to achieve is to calculate the area of a quadrilateral by the
> vector cross product of its diagonals. X, Y and Z contain the points
> that describe my mesh. I am sure that there is a way to vectorise at
> least part of this, and I would really appreciate any suggestions,
>
> [R C] = size(X); % all matrices the same size
> row = 0;
> for i = 1:(R-1)
> for j = 1:(C-1)
> P4 = [X(i,j) Y(i,j) Z(i,j)];
> P3 = [X(i,j+1) Y(i,j+1) Z(i,j+1)];
> P2 = [X(i+1,j+1) Y(i+1,j+1) Z(i+1,j+1)];
> P1 = [X(i+1,j) Y(i+1,j) Z(i+1,j)];
> p = P2-P4;
> q = P3-P1;
> row++
> r(row,:) = 0.5 * (cross(p,q)); % calculate the area
> endfor
> endfor



Here is a mechanical translation of your loop, producing 3-column
matrices p and q:

[m,n] = size(X);
I = 1:m-1;
J = 1:n-1;
P4 = [X(I,J)(:), Y(I,J)(:), Z(I,J)(:)];
P3 = [X(I,J+1)(:),Y(I,J+1)(:), Z(I,J+1)(:)];
P2 = [X(I+1,J+1)(:),Y(I+1,J+1)(:),Z(I+1,J+1)(:)];
P1 = [X(I+1,J)(:),Y(I+1,J)(:),Z(I+1,J)(:)];
p = P2-P4;
q = P3-P1;
r = 0.5 * cross(p,q);

Because you traverse by column rather than by row the results will be in a different order than you give. If for some reason you need your particular order, you can transpose X,Y,Z and swap P3 and P1.

- Paul


Now I wish I could do that, and what we call optimization.
-Muthu