MATLAB Answers

How can vectorize the below codes?

2 views (last 30 days)
Nx=512;
Ny=Nx
x0 = Nx/2;
y0 = Ny/2;
isolve = 2;
tic
for i=1:512
for j=1:512
ii=(i-1)*Nx+j;
if(isolve == 2)
etas(ii,1)=1.0;
phis(ii,1) =0;
phis(ii,2)=0;
etas(ii,2)=0.0;
else
etas(i,j,1) =1.0;
etas(i,j,2) =0.0;
phis(ii,1) =0;
phis(ii,2)=0;
end
xlength =sqrt((i-x0)^2+(j-y0)^2);
if(xlength <= 14)
if(isolve == 2)
etas(ii,1)=0.0;
etas(ii,2)=1.0;
phis(ii,1) =0;
phis(ii,2)=0;
else
phis(ii,1) =0;
phis(ii,2)=0;
end
end %if
end %j
end %i

  0 Comments

Sign in to comment.

Accepted Answer

the cyclist
the cyclist on 1 Dec 2019
Well, not vectorizing it, but if you preallocate etas and phis with these line:
etas = zeros(262144,2);
phis = zeros(262144,2);
you will speed up the code by a huge factor. (For me, it was a factor of over 1000.)

  5 Comments

Show 2 older comments
the cyclist
the cyclist on 2 Dec 2019
Assuming that the answer to my question above is "yes", then there are massive simplifications possible for your code. If you and I were sitting together, it would be very easy for you to see each step that I did -- preallocating values, deleting lines of code that because unnecessary, consolidating certain statements, etc.
Because we are not sitting together, I think it might be perplexing how I got to such simpler code. I wrote a list of comments for each thing I did, which may help.
Nx = 512;
Ny = Nx;
x0 = Nx/2;
y0 = Ny/2;
isolve = 2;
% 1: preallocate
NxNy = Nx*Ny;
etas = [ones(NxNy,1) zeros(NxNy,1)]; % First column of etas is usually 1. % Second is usually 0;
phis = zeros(NxNy,2); % Both columns of phis are zero.
% 2: since preallocation sets most values, remove all lines that were
% assigning those same values later
% 3: the first "if" statement has the same effect as its else statement, so
% just always do that operation of setting etas(ii,1)=1. But that statement
% gets remove altogether later, because of preallocation of the first column.
% 4: simplified the second if statement, by consolidating two ifs into
% one
% 5: tiny vectorization of an assignment statement for etas
% 6: first column of etas is actually all ones, so preallocate to that
% 7: remove the assignment etas(ii,1) = 1 from the loop, since it has been
% preallocated
tic
for i = 1:Nx
for j = 1:Ny
ii = (i-1)*Nx+j;
xlength = sqrt((i-x0).^2+(j-y0).^2);
if (isolve == 2) && (xlength <= 14)
etas(ii,:) = [0 1];
end %if
end %j
end %i
toc
I tested this code for the values of isolve equal to 1 and 2, and got the same results as your code.
Conceptually, it is now much easier to see what is happening. It seems like:
  • Every row of phis is equal to [0 0]
  • Every row of etas is equal to [1 0], but then some get flipped to [0 1], based on the value of isolve and xlength
Does that seem right?
the cyclist
the cyclist on 2 Dec 2019
The following is a vectorized version. It is not that much faster, and actually more difficult to understand what is going on.
Nx = 512;
Ny = Nx;
x0 = Nx/2;
y0 = Ny/2;
isolve = 2;
NxNy = Nx*Ny;
etas = [ones(NxNy,1) zeros(NxNy,1)];
phis = zeros(NxNy,2);
tic
xvec = 1:Nx;
yvec = 1:Ny;
xlength = sqrt((xvec-x0).^2+(yvec'-y0).^2);
idx = isolve==2 & xlength(:) <= 14;
etas(idx,:) = repmat([0 1],sum(idx),1);
toc
Amir Torabi
Amir Torabi on 2 Dec 2019
very grateful. It really works great!.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!