Vectorization of a double for loop
显示 更早的评论
Hello all, I am trying to speed up the assembly of a matrix with a complicated structure. Right now, I have several sets of nested loops such as
for i=M+2:M+N-2
for j=i+1:M+N-1
BL1(j,i) = gamm/(gbp1*h(i)^(1-beta))*...
((j-i+3/2)^beta-2*(j-i+1/2)^beta+(j-i-1/2)^beta);
end
end
Here gamm, beta, are constants, h is a vector of length M+N and BL1 is a (dense) matrix of size (M+N-1) x (M+N-1). Is it possible to vectorize this? I've tried several things, but none seem to work.
Thanks for taking a look.
3 个评论
dpb
2017-11-9
gbpi also a constant?
Geoffrey Dillon
2017-11-9
Jan
2017-11-10
"gbpi"? It is gbp1.
采纳的回答
Roger Stafford
2017-11-10
The following code should correct the error in my first answer. This version also avoids the repetition involved in your code for the expression
((j-i+3/2)^beta-2*(j-i+1/2)^beta+(j-i-1/2)^beta)
and hopefully it should therefore take less execution time.
BL1 = zeros(M+N-1);
T = (1:N-3)';
X = (T+3/2).^beta - 2*(T+1/2).^beta + (T-1/2).^beta;
for i = M+2:M+N-2
BL1(T(1:M+N-1-i)+i,i) = gamm/(gbp1*h(i)^(1-beta))*X(1:M+N-1-i);
end
更多回答（1 个）
Roger Stafford
2017-11-9
编辑：Roger Stafford
2017-11-10
You can save unneeded repetition by computing one of the factors outside the loop.
BL1 = zeros(M+N-1);
T = (1:M+N-1)';
X = (T+3/2).^beta - 2*(T+1/2).^beta + (T-1/2).^beta;
for i=M+2:M+N-2
BL1(T+i,i) = gamm/(gbp1*h(i)^(1-beta)) * X;
end
This answer is incorrect. Read my correction in the revised answer below. RAS
5 个评论
Jan
2017-11-10
EDITED: The quote was not ', but a similar character. But not ´ . Strange.
Roger Stafford
2017-11-10
@Jan: Thanks for the editing. I thought I typed in an ordinary, garden-variety single quote. Perhaps my computer is playing tricks on me.
Geoffrey Dillon
2017-11-10
dpb
2017-11-10
As written for the above M,N --> BL1 will be M+N-1 x M+N-2 --> 4127x4126 as those are the upper limits for the i,j loops respectively and the indices for BL1 are those indices identically. Every row below M+2 (=4098) will be identically zero and all columns diagonally to the left will also be zero beyond that row.
Is that the intent? But if it is to be square and 4127x4127, then the i upper index is wrong or the column subscript expression isn't correct.
Did you preallocate BL1 before timing the loop solution? Awaiting the answer to the above observation I've not looked at whether there's an efficient vectorization to be done or not; often the looping solution can be pretty effective if do the obvious.
Roger Stafford
2017-11-10
@Geoffrey: My apologies. I mistakenly read the line
for j=i+1:M+N-1
as though it had the parentheses:
for j=i+(1:M+N-1)
I have written a second “answer” which should correct this error.
类别
在 帮助中心 和 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 
