Vectorizing/Profiling Code with double for loops etc

Hi,

Need help increasing the efficiency of the code below (getting rid of double for-loops etc.) I saw some topics relating to bsxfun, but im not too sure if thats applicable here. I appreciate all help in advance. Thank you.

EDIT: Tried changing some other things to vectorize this, but could not make it happen. END Edit

Here's the code:

if true
  % code
end
clear;clc;close all
tic
cylgrid = ones(201,502);
init = 100:-1:-100; 
i=1:201;
for j=1:502
cylgrid(:,j)=init;
end
%initialize circle
j=500;
for ii = -100:100 %fluctuate i
  for jj= 1:502 %fluctuate j
      if sqrt( (ii ).^2 +(jj - j./5).^2) < 25
  cylgrid(ii+101,jj) = 0;%psi = 0 at surface and inside cylinder
      end
  end
end
psi = cylgrid;%just for notation
psinew= cylgrid;
error = 0.001;
t=0;dt=1;
while error < max(max(psinew-psi)) || t==0
  psi=psinew;
for ii = 2:200
  for jj = 2:500
      psinew(ii,jj) = psi(ii,jj) +(1/4).*(psi(ii-1,jj)+psi(ii+1,jj)+psi(ii,jj-1)+psi(ii,jj+1)-4.*psi(ii,jj));
  end
end
%making sure we begin with our circle again...
for ii = -100:100 %fluctuate i
  for jj= 1:502 %fluctuate j
      if sqrt( (ii ).^2 +(jj - j./5).^2) < 25
  psinew(ii+101,jj) = 0;%psi = 0 at surface and inside cylinder
      end
  end
end
t=t+dt;
end
uprofile= abs(psinew);
for iii = 100:-1:1
  for jjj=1:502
      uprofile(101-iii,jjj) = uprofile(101-iii,jjj)./iii;
      uprofile(iii+101,jjj) = uprofile(iii+101,jjj)./iii;
  end
end

figure(1); title('Streamline Plot');

contourf(psinew,100);set(gca,'DataAspectRatio',[1 1 1]);

figure(2);title('Velocity Profile');

contourf(uprofile,10);set(gca,'DataAspectRatio',[1 1 1]);

toc

[E]

Help? If something is unclear, let me know and I can clarify.

Answer 1

[Geoff Hayes]

Hi E,

Within your while loop, there is the following chunk of code:

%making sure we begin with our circle again...
for ii = -100:100 %fluctuate i
    for jj= 1:502 %fluctuate j
        if sqrt( (ii ).^2 +(jj - j./5).^2) < 25
            psinew(ii+101,jj) = 0;%psi = 0 at surface and inside cylinder
        end
    end
end

Unless I'm misreading something, the ii and jj that satisfy the condition sqrt( (ii ).^2 +(jj - j./5).^2) < 25 will always be the same (and are no different than those that were used in the similar block when the circle is initialized). If that is the case, then why not just keep track of these indices and re-use them, rather than doing this double loop again and again for the same set? When the circle is initialized, just add the following matrix to keep track of the indices that are being set to zero:

% initialize matrix to keep track of indices that should be reset to zero
indcsToReset = cast(zeros(size(cylgrid)),'logical');

%initialize circle
j=500;
for ii = -100:100 %fluctuate i
  for jj= 1:502 %fluctuate j
      if sqrt( (ii ).^2 +(jj - j./5).^2) < 25
          cylgrid(ii+101,jj) = 0;%psi = 0 at surface and inside cylinder

            % reset this index
            indcsToReset(ii+101,jj) = true; 
        end
    end
  end

Then you can replace the first chunk above with the comment making sure we begin with our circle again with

psinew(indcsToReset)=0;

Hopefully the above is a start to improving performance…

Geoff

[E]

Thanks! I was wondering how to do something like that. Never understood what the 'cast' command does. I'll check it out now.

Any ideas to get rid of the double for-loop below? I don't think bsxfun will do it in this case. But not sure.

if true
  % code
end
for ii = 2:200
for jj = 2:500
    psinew(ii,jj) = psi(ii,jj) +(1/4).*(psi(ii-1,jj)+psi(ii+1,jj)+psi(ii,jj-1)+psi(ii,jj+1)-4.*psi(ii,jj));
end
end

[Geoff Hayes]

A quick glance - you should remove the periods (i.e. (1/4).*(psi…) and in the last operation since you are just doing scalar multiplication (no matrix element-wise multiplication going on).

and doesn't psi(ii,jj)* cancel out because you are subtracting it again:

psinew(ii,jj) = psi(ii,jj) + (1/4)*(….. - 4*psi(ii,jj));

If that is the case, your equation here can be reduced to:

0.25*(psi(ii-1,jj)+psi(ii+1,jj)+psi(ii,jj-1)+psi(ii,jj+1));

[E]

Yes it does. But in the later problem when I introduce viscosity, I will use a relaxation parameter, so instead of 1/4, it will be something like 1.5/4, where 1 or 1.5 is the relaxation parameter.

The first comment helped cut my run time in half! So if this second double for-loop can be consolidated in some way, then the result can be obtain even faster!

[Geoff Hayes]

Cool. Looking at that line of code, it seems to be just summing the neighbours (above, below, left, right) of the current cell given by the index (ii,jj) then subtracts 4 times itself before adding it all to itself. If that is the case, and this is always true then it could be replaced with the following:

% get the number of rows and columns in psi (r==201,c==502)
[r,c] = size(psi);

while error < max(max(psinew-psi)) || t==0
    psi=psinew;

      % avoid the double loop altogether:
      psinew(2:r-1,2:c-2) = psi(2:r-1,2:c-2) + 0.25*(psi(1:r-2,2:c-2) + …
                            psi(3:r,2:c-2) + psi(2:r-1,3:c-1) +         …
                            psi(2:r-1,1:c-3) - 4*psi(2:r-1,2:c-2)); 

      % making sure we begin with our circle again…
      % etc.
end

So in the above, shifted copies of psi are added to one another (hopefully) simulating what the double for loop was doing.

Question - why does the inner for loop iterate from 2 to 500 and not 501 since the number of columns in psinew is 502?

[E]

Thanks! I'll try this out right now. Response to your question..I goofed, but it won't change much about the fluid flow.

[E]

WOW ! this was ridiculously good! From my original code, which took about 1 hour to run, this code now took 3.5 seconds!

I am amazed!... Kudos to you and thank you very much!

[E]

I'd been trying to figure out how to do that "psinew" line of code for a long time now. I kept getting errors. I really appreciate this!

[Geoff Hayes]

Glad to have been able to help!

Answer 2

[Jan]

A general hint: SQRT() is expensive, so try to avoid it:

sqrt( (ii ).^2 +(jj - j./5).^2) < 25

Faster:

ii.^2 + (jj - j ./ 5).^2 < 625

[E]

Did not know that. Thanks! Will implement that now.