MATLAB Answers

## Matlab for loop vectorization

Asked by shir

### shir (view profile)

on 9 Dec 2012
Accepted Answer by Matt J

### Matt J (view profile)

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;
```

Matt J

### Matt J (view profile)

on 9 Dec 2012

Explain this line,

` VEC(i0,j0)=VEC(i0-1,j0-1);`

Everything is very easy to vectorize apart from that. I know it's supposed to handle the case VEC(i0,j0)=0, but this looks like it might be an improvised way of dealing with it. What is zv(i0,j0) really supposed to be?

Matt J

### Matt J (view profile)

on 9 Dec 2012

Also, is size(X) equal to [x00,y00]? That would simplify things, too.

shir

### shir (view profile)

on 9 Dec 2012

is size(X) equal to [x00,y00]?-yes.

VEC(i0,j0) shouldnt be equal to zero as you can guess from the last line, where it is in denominator.

zv is velocity of the surface

Log in to comment.

## 1 Answer

Answer by Matt J

### Matt J (view profile)

on 9 Dec 2012
Accepted Answer

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);```

shir

### shir (view profile)

on 9 Dec 2012

thank you for your answer. And what if i0,j0 goes throgh arbitrary values like

```      for i0=A-10:1:A+10
for j0=B-10:1:B+10```

without bothering about index out of bounds issues.

Matt J

### Matt J (view profile)

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.

Log in to comment.

### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today