I'm looking to eliminate the following for loop:

2 次查看(过去 30 天)
I'm looking to eliminate the following for loop:
for j=1:96
rn(ijk(j))=rn(ijk(j))+mp(j);
end
Background: rn is a 200x3x3 matrix, ijk is a 1800x1 matrix and np is a 12x8 matrix. ijk describes the mapping between mp and (the much larger) rn. For example, when j=1, mp(j)=0.84 and ijk(1)=601. This loop is embedded in a much larger loop.
There has got to be a way to eliminate the for loop presented here. Any suggestions?

采纳的回答

Sean de Wolski
Sean de Wolski 2012-1-13
Perhaps:
idx = 1:96;
rn(ijk(idx))=rn(ijk(idx))+mp(idx);
  3 个评论
Alex
Alex 2012-1-13
I just played around with a similar setup and I'm getting the same answer using either method.
One way that it could be causing different answers is if the data types are different (ex. doubles vs cells).
You can further debug this to see what is going on by reducing the problem. Start with a smaller idx (even 1) to see if you are getting the same results both ways. If 1 works, make idx bigger until you find the problem.
Sean de Wolski
Sean de Wolski 2012-1-13
Could you try each method in its entirety in two different functions and then compare the results after? Perhaps you're performing the same operation twice or not preallocating the matrices the same way (i.e. different sized results)

请先登录,再进行评论。

更多回答(2 个)

John
John 2012-1-18
Here was the problem: my vector ijk pointed to same index in rn at multiple times. For example, consider
mp=[0.1;0.2;0.3]; ijk=[1,2,1]; rn=zeros(4,1);
The vectorized function gives: idx = 1:3; rn(ijk(idx))=rn(ijk(idx))+mp(idx);
rn =
0.3000
0.2000
0
0
The loop expression gives:
for j=1:3 rn(ijk(j))=rn(ijk(j))+mp(j); end
rn =
0.4000
0.2000
0
0
I'm still not sure how to vectorize the code for this special case, but the answers provided by Sean and Alex are correct to the question.

Sean de Wolski
Sean de Wolski 2012-1-18
%Data:
mp=[0.1;0.2;0.3];
ijk=[1,2,1];
rnsz = [4 1]; %the size you want rn to be:
%Engine:
rn = accumarray(ijk',mp',rnsz); % Build rn; note transpose of subs/vals

类别

Help CenterFile Exchange 中查找有关 Matrix Indexing 的更多信息

产品

Community Treasure Hunt

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

Start Hunting!

Translated by