[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Request help speeding up script
From: |
Paul Kienzle |
Subject: |
Re: Request help speeding up script |
Date: |
Mon, 25 Jun 2007 22:34:05 -0400 |
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