How can vectorize the below codes?

Question

Amir Torabi 2019-12-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494208-how-can-vectorize-the-below-codes

评论： Amir Torabi 2019-12-2

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

the cyclist 2019-12-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494208-how-can-vectorize-the-below-codes#answer_404165

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

the cyclist 2019-12-2

在 MATLAB Online 中打开

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 2019-12-2

在 MATLAB Online 中打开

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