How can vectorize the below codes?

2 views (last 30 days)
Amir Torabi on 1 Dec 2019
Commented: Amir Torabi on 2 Dec 2019
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

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.)

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 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 on 2 Dec 2019
very grateful. It really works great!.