Matlab for loop vectorization
1 view (last 30 days)
Show older comments
X,Y and z are coordinates representing surface. In order to calculate some quantity, lets call it flow, at point i,j of the surface, i need to calculate contibution from all other points (i0,j0). To do so i need for example to know cos of angles between point i0,j0 and all other points (alpha). Then all contirbutions from i0,j0 must be multiplied on some constants and added. zv0 at every point i,j is final needed result.
I came up with some code written below and it seems to be extremely unappropriate. First of all it slows down rest of the program and seems to use all of the available memory.My system has 4gb physical memory and 12gb swap file and it always runs out of memory, though all of variables sizes are not bigger then 10kb. Please help up with speed up/vectorization and memory problems.
parfor i0=2:1:length(x00);
for j0=2:1:length(y00);
zv=red3dfunc(X0,Y0,f,z0,i0,j0,st,ang,nx,ny,nz);
zv0=zv0+zv;
end
end
function[X,Y,z,zv]=red3dfunc(X,Y,f,z,i0,j0,st,ang,Nx,Ny,Nz)
x1=X(i0,j0);
y1=Y(i0,j0);
z1=z(i0,j0);
alpha=zeros(size(X));
betha=zeros(size(X));
r=zeros(size(X));
XXa=X-x1;
YYa=Y-y1;
ZZa=z-z1;
VEC=((XXa).^2+(YYa).^2+(ZZa).^2).^(1/2);
VEC(i0,j0)=VEC(i0-1,j0-1);
XXa=XXa./VEC;
YYa=YYa./VEC;
ZZa=ZZa./VEC;
alpha=-(Nx(i0,j0).*XXa+Ny(i0,j0).*YYa+Nz(i0,j0).*ZZa);
betha=Nx.*XXa+Ny.*YYa+Nz.*ZZb;
r=VEC;
zv=(1/pi)*st^2*ang.*f.*(alpha).*betha./r.^2;
3 Comments
Accepted Answer
Matt J
on 9 Dec 2012
Here's one idea, though I'm assuming for now that you want to handle VEC=0 situations using the rule 0/0 = 1.
[m,n]=size(X);
i0=2:length(x00);
j0=2:length(y00);
x1=X(i0,j0);
y1=Y(i0,j0);
z1=z(i0,j0);
X=X(:); Y=Y(:); Z=Z(:); %make everything column vector
XXa=bsxfun(@minus,X,x1(:).');
YYa=bsxfun(@minus,Y,y1(:).');
ZZa=bsxfun(@minus,Z,Z1(:).');
VEC=((XXa).^2+(YYa).^2+(ZZa).^2).^(1/2);
XXa=XXa./VEC.^2;
YYa=YYa./VEC.^2;
ZZa=ZZa./VEC.^2;
idx=isnan(XXa);
XXa(idx)=1;
YYa(idx)=1;
ZZa(idx)=1;
betha = bsxfun(@times,Nx(:),XXa) + ...
bsxfun(@times,Ny(:),YYa) + ...
bsxfun(@times,Nz(:),ZZa);
zv=(1/pi)*st^2*ang.*f.*(betha*diag(betha));
zv=reshape(zv,m,n);
2 Comments
Matt J
on 9 Dec 2012
If the indices goes out of bounds, then lines like X(i0,j0) will fail, but making them subsets of 1:m, 1:n shouldn't matter.
More Answers (0)
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!