How can I vectorize two nested loops with different dimensions?
1 回表示 (過去 30 日間)
古いコメントを表示
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.
1 件のコメント
Jan
2019 年 6 月 5 日
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?
採用された回答
Jan
2019 年 6 月 5 日
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.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Direct Search についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!