How can I vectorize two nested loops with different dimensions?

In my code I have two nested loops. How can I vectorize them?

The first nested loop is running on all triangles and all edges in order to verify if satisfy the following condition: and(t(4,i)==2,e(6,j)==1).

e_aux = e;                           % I get t and e from [p,e,t] = initmesh(g), where g is the domain.
for i=1:length(t)                    % Length of the connectivities of t
    for j=1:length(e)                % Length of the connectivities of t
        if and(t(4,i)==2,e(6,j)==1)
            e(1,j) = e_aux(2,j);
            e(2,j) = e_aux(1,j);
            e(6,j)=e_aux(7,j);
            e(7,j)=e_aux(6,j);
        end
    end
end

The second nested loop is running on all nodes of the mesh (Nnodes) and all nodes of the subset domain (Npanels).

for i = 1:2
    for j = 1:Nnodes                                           % Number of nodes in all domain
        for k = 1:Npanels                                      % Number of nodes in a subset domain
            dist = sqrt((xm(k)-x(j))^2 + (ym(k)-y(j))^2);      % Distance between two coordinates (rm and r)
            Green= log(dist)/(2*pi);                           % Green function
            phi(j,i) = phi(j,i) + f(k,i)'*lg(k)*Green/(mur-1); % f(k,i) is a vector with Npanels x 2 dimension
        end
    end
end

In both case, when I have a fine mesh the code gets extremely slow.

[Jan]

Please provide some meaningful input data, e.g. created by rand(). length(t) is fragile: If t is a [4 x 3] matrix for any reason, its length is 4. Use size(t, 2) to determine the length of the 2nd dimension.

In the first code, you overwrite the elements of e repeatedly. Is this really wanted?

Answer 1

[Jan]

The first loop looks strange, because the body does not depend on the loop counter i. Currently this should produce the same output, but the dependency to t seems to be strange:

if any(t(4, :) == 2)
  match = (e(6, :) == 1);
  e([1,2,6,7], match) = e([2,1,7,6], match);
end

Here the innermost loop is vectorized:

for i = 1:2
    for j = 1:Nnodes
        dist     = sqrt((xm - x(j)) .^ 2 + (ym - y(j)) .^2);
        Green    = log(dist) / (2*pi);
        phi(j,i) = sum(f(:, i)' .* lg * Green ./ (mur-1));
    end
end

Without having some input data for testing, the optimization of the code is really hard. So please povide some such that we can test the suggestions before posting them.

[Alan]

Hi Jan,

Thank you for the answer. I am including all data below.

I start with a domain g. Then I create a mesh using initmesh. After obtaining the connectivities I look for the middle points (is it possible generate these middle points vectorized as well?). In addition, I solve a linear system to obtain f. Finally, I have the second nested loop.

% Domain
g = ...
  [2 2 2 4 4 4 4;
   1 0 -1 -19.924858845171276 -0.43352184697054014 19.981205880819925 1.2997507747857209 ...
   ;
   0 -1 1 -0.43352184697054014 19.981205880819925 1.2997507747857209 -19.924858845171276 ...
   ;
   0 1.7320508075688772 0 1.7320508075688732 -19.995300918170731 -0.86683997847771355 ...
   19.957721511320972;
   1.7320508075688772 0 0 -19.995300918170731 -0.86683997847771355 19.957721511320972 ...
   1.7320508075688732;
   2 2 2 1 1 1 1;
   1 1 1 0 0 0 0;
   0 0 0 0 0 0 0;
   0 0 0 0 0 0 0;
   0 0 0 20 20 20 20;
   0 0 0 20 20 20 20;
   0 0 0 0 0 0 0];

% mesh 
[p,e,t] = initmesh(g,'hmax', 25);
for i=1:2
        [p,e,t] = refinemesh(g,p,e,t); % regular mesh
end

% mesh data
x = p(1,:);
y = p(2,:);

connect_bound_obj = find(e(6,:)==2);
e_obj = e(:,connect_bound_obj);
Npanels = length(e_obj); % Number of edges
Nnodes = length(p);      % Number of nodes

% Getting the middle points and the normal
for i=1:Npanels
   xm(i)=0.5*(x(e_obj(1,i))+x(e_obj(2,i))); % x-middle point edge
   ym(i)=0.5*(y(e_obj(1,i))+y(e_obj(2,i))); % y-middle point edge
   lg(i)=sqrt((x(e_obj(2,i))-x(e_obj(1,i)))^2+(y(e_obj(2,i))-y(e_obj(1,i)))^2); % length edge
   nx(i)=(y(e_obj(2,i))-y(e_obj(1,i)))/lg(i); % normal x
   ny(i)=(-x(e_obj(2,i))+x(e_obj(1,i)))/lg(i); % normal y
end


A= rand(Npanels); % lambda*I - K^star
normal = [nx; ny];
f = A\normal'; % vector solution phi


for i = 1:2
    for j = 1:Nnodes                                           % Number of nodes in all domain
        for k = 1:Npanels                                      % Number of nodes in a subset domain
            dist = sqrt((xm(k)-x(j))^2 + (ym(k)-y(j))^2);      % Distance between two coordinates (rm and r)
            Green= log(dist)/(2*pi);                           % Green function
            phi(j,i) = phi(j,i) + f(k,i)'*lg(k)*Green/(mur-1); % f(k,i) is a vector with Npanels x 2 dimension
        end
    end
end

Is it possible understand the code? Let me know if do you understand.