How can i optimize my code without a for-loop ?

Hello, I've partially written a script which uses a for loop. But the problem is taht it slows too much the program, could you look at it and tell me if it is possible to use only matrix calculus ?

Here is the code:

m = [0 0 m0]; %m0 is a constant
R_rec = repmat(r_rec,size(r_source,1),1); %the dimension is 78000...
R = r_source - R_rec;

for n=1:size(R)
    r = R(n,:);
    r1 = r/sqrt(r*r');
    B0(n,:) = 1/(sqrt(r*r'))^(3/2).*((3*dot(m,r1).*r1) - m);
end

S = sum(B0,1);

Bx = S(1);
By = S(2);
Bz = S(3);
end

The puropose of this script is to calculate the magnetic field of magnet (r_source) at a point (r_sec) r_source is 78000*3 and r_sec = 1*3 that's why I use repmat() to get the position vector r used in the formula. It takes 2-3 minutes to compute the field about +/- 100 points, which is too much...

How can I optimize it ?

Thanks a lot, Paul.

Ps: sorry for my bad English, it's not my mother tongue.

Answer 1

[Ced]

How about this? With 100k rows, I get old version: 0.90 s new version: 0.01 s

m = [0 0 m0]; %m0 is a constant
R_rec = repmat(r_rec,size(r_source,1),1); %the dimension is 78000...
R = r_source - R_rec;

r_norm = sqrt(sum(R.*R,2));
R_norm = R./repmat(r_norm,1,3);
v = 3*m0*R_norm(:,3);
B0 = R_norm.*repmat(v,1,3);
B0(:,3) = B0(:,3)-m0;
B0 = B0./(repmat(r_norm.^(3/2),1,3));

S = sum(B0,1);

Bx = S(1);
By = S(2);
Bz = S(3);

Note that this version will only work if m(1) = m(2) = 0. But I think you won't have any problem in adjusting it if that is not the case.

Cheers

PS: Shouldn't it be sqrt(r*r'))^3 instead of sqrt(r*r'))^(3/2) at the end? But maybe I'm missing something.

[Paul Eluard]

Yes, thanks a lot ! It takes about 0.011 sec on my computer so you're the fastest and yes, it is (r*r')^(3/2) and not with the square root...

Answer 2

[Andrei Bobrov]

R = bsxfun(@minus, r_source, r_rec);

sqr1 = 1./sqrt(sum(R.^2,2));
r1 = bsxfun(@times,R,sqr1);
r2 = 3*r1(:,3)*m0;
r3 = bsxfun(@times,r2,r1);
r3(:,3) = r3(:,3) - m0;

S = sum(bsxfun(@times,sqr1.^1.5,r3));

[Paul Eluard]

It also works but i got a faster answer; yours takes about 0.013 but is very clear and easy to understand, thank you