How can vectorize the below codes?

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

Answer 1

[the cyclist]

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

[Amir Torabi]

very appreciated. it workd and has a huge effect on reducing time. I have another question!. is this code be able to be written in the form parrallel?

[the cyclist]

I'm not sure if it can easily be vectorized, but it can certainly be simplified in places.

If the code is easy to understand, and sufficiently fast, I would not be too concerned about vectorization. Sometimes that can make the code more difficult to comprehend, which is bad.

I have one question. I notice the lines

etas(i,j,1) =1.0;
etas(i,j,2) =0.0;

Should those actually be

etas(ii,1) =1.0;
etas(ii,2) =0.0;

?

[the cyclist]

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 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]

very grateful. It really works great!.