Vectorization of a double for loop
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Geoffrey Dillon
le 9 Nov 2017
Commenté : Geoffrey Dillon
le 13 Nov 2017
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 commentaires
Réponse acceptée
Roger Stafford
le 10 Nov 2017
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
Plus de réponses (1)
Roger Stafford
le 9 Nov 2017
Modifié(e) : Roger Stafford
le 10 Nov 2017
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 commentaires
dpb
le 10 Nov 2017
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
le 10 Nov 2017
@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.
Voir également
Catégories
En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!