Eliminate for loop: can this code be vectorized?

1 次查看（过去 30 天）
Quinten Rensen
Quinten Rensen 2018-7-27
评论： Quinten Rensen 2018-7-28
Hi all,
I am trying to optimize my code, my question is if it is possible to vectorize the following code:
for r = 1:2:(N-1)
U(r,1)=U(r,1)+ (1./K(r,r)).*(F(r,1)-K(r,:)*U(:,1));
end
where U and F are Nx1 vectors and K is a NxN matrix. This part of the code takes a considerable amount of time when N becomes large. Maybe someone has a suggestion on how to vectorize this part?
Thank you!

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

采纳的回答

KSSV
KSSV 2018-7-27
编辑：KSSV 2018-7-27
r = 1:2:(N-1) ;
U(r,1)=U(r,1)+ (1./K(r,r))*(F(r,1)-K(r,:)*U(:,1));
unchecked check the values before using.
  6 个评论
显示 4更早的评论
Quinten Rensen
Quinten Rensen 2018-7-27
编辑：Quinten Rensen 2018-7-27
Hi Jan, thank you for your quick replies. Maybe it is useful to note that I have defined K as a sparse matrix, so K = [(indices) , value ]
Here is the original code:
function [U,res]=GS(U,F,nelx,nely,sweeps,K)
omega = 1;
count =1;
N = 2.*(nelx+1).*(nely+1);
while count <= sweeps %Number of sweeps;
for r = 1:2:(N-1)
U(r,1)=U(r,1)+omega.*(1./K(r,r)).*(F(r,1)-K(r,:)*U(:,1));
end
for c = 2:2:N
U(c,1)=U(c,1)+omega.*(1./K(c,c)).*(F(c,1)-K(c,:)*U(:,1));
end
count = count + 1;
end
%%Calculate final residual at all points
res = zeros(2.*(nelx+1).*(nely+1),1); %Preallocate
res = F - K*U; %All residuals;
end%End of GS function
And this is the code after I implemented the suggestions:
function [U,res]=GS(U,F,nelx,nely,sweeps,K)
omega = 1;
count =1;
N = 2.*(nelx+1).*(nely+1);
while count <= sweeps %Number of sweeps;
r = 1:2:(N-1)
U(r,1)=U(r,1)+omega.*(1./K(sub2ind(size(K), r, r))).*(F(r,1)-K(sub2ind(size(K), r, :))*U(:,1));
c = 2:2:N
U(c,1)=U(c,1)+omega.*(1./K(sub2ind(size(K), c, c))).*(F(c,1)-K(sub2ind(size(K), c, :))*U(:,1));
count = count + 1;
end
%%Calculate final residual at all points
res = zeros(2.*(nelx+1).*(nely+1),1); %Preallocate
res = F - K*U; %All residuals;
end%End of GS function
Quinten Rensen
Quinten Rensen 2018-7-28
Ok, I found what was wrong, I had to use the transpose of (1./K(sub2ind(size(K), r, r))). Then the code becomes
r = 1:2:(N-1)
U(r,1)=U(r,1)+omega.*(1./K(sub2ind(size(K), r, r)))'.*(F(r,1)-K(r,:)*U(:,1));
Thank you for your help!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by