CODEWITHSUNDEEP
performancematlabfor-loopvectorization
noir
noir

Reputation: 185

Vectorize double for loops

I need to evaluate an integral, and my code is

r=0:25;
t=0:250;
Ti=exp(-r.^2);
T=zeros(length(r),length(t));
for n=1:length(t)
  w=1/2/t(n);
  for m=1:length(r)
    T(m,n)=w*trapz(r,Ti.*exp(-(r(m).^2+r.^2)*w/2).*r.*besseli(0,r(m)*r*w));
  end
end

Currently the evaluation is fairly fast, but I wonder if there is a way to vectorize the double for-loop and make it even faster, especially when function trapz is used.

Upvotes: 1

Views: 97

Answers (1)

Oleg
Oleg

Reputation: 10676

You can optimize it by passing matrix argument Y to trapz(A,Y), and using dim = 2, i.e. the loop becomes:

r  = 0:25;
t  = 0:250;
Ti = exp(-r.^2);

tic
T  = zeros(length(r),length(t));
for n = 1:length(t)
    w = 1/2/t(n);
    for m = 1:length(r)
        T(m,n) = w*trapz(r,Ti.*exp(-(r(m).^2+r.^2)*w/2).*r.*besseli(0,r(m)*r*w));
    end
end
toc

tic
T1  = zeros(length(r),length(t));
for n = 1:length(t)
    w = 1/2/t(n);
    Y = bsxfun(@times,Ti.*r, exp(-bsxfun(@plus,r'.^2,r.^2)*w/2).*besseli(0,bsxfun(@times,r',r*w)));
    T1(:,n) = w* trapz(r,Y,2);
end
toc

max(abs(T(:)-T1(:)))

You could probably vectorize it completely, will have a look later.

Upvotes: 2

Related Questions