Re: Request help speeding up script

[Paul Kienzle]

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