MATLAB Central Newsreader - urgent: for loops

urgent: for loops

adinda — Tue, 30 Oct 2007 13:11:32 -0400

I have this huge loop in matlab, which of course asks a lot of
computation time. I can't find another way because i need the
differences between the steps in Q, and also in my formula, but I'm
not a genius with matlab. can anyone give me advice on how to remove
my loops?
Thanks

Adinda

The loop:
Hms=zeros(nx,ny,nz,3);
Hmsnk=zeros(nx,ny,nz);
gradHms=zeros(nx,ny,nz,3,6);
%for loop with calculation of the magnetostatic field and its
derivative
for i=1:nx
    for j=1:ny
        for s=1:nz
            for k=1:nx
                for l=1:ny
                    for m=1:nz
                        if (i~=k && j~=l && s~=m)

                            Q=sqrt((k-i)^2+(l-j)^2+(m-s)^2);
                            Hms(i,j,s,1)=Hms(i,j,s,1)+(M(k,l,m,1)/
Q^3-3*(k-i)/Q^5*(M(k,l,m,1)*(k-i)+M(k,l,m,2)*(l-j)+M(k,l,m,3)*(m-s)));

                            Hms(i,j,s,2)=Hms(i,j,s,2)+(M(k,l,m,2)/
Q^3-3*(l-j)/Q^5*(M(k,l,m,1)*(k-i)+M(k,l,m,2)*(l-j)+M(k,l,m,3)*(m-s)));

                            Hms(i,j,s,3)=Hms(i,j,s,3)+(M(k,l,m,3)/
Q^3-3*(m-s)/Q^5*(M(k,l,m,1)*(k-i)+M(k,l,m,2)*(l-j)+M(k,l,m,3)*(m-s)));

                            %gradient van H_ms_x naar de 6 volumes
                            gradHms(i,j,s,1,1)=sum(sum(sum((1/
Q^3-3*(k-
i)^2/Q^5).*difMx1(k,l,m)+(-3*(k-i)*(l-j)/Q^5).*difMy1(k,l,m)...
                                +(-3*(k-i)*(m-s)/
Q^5).*difMz1(k,l,m))));

                            gradHms(i,j,s,1,2)=gradHms(i,j,s,1,2)+(1/
Q^3-3*(k-i)^2/Q^5).*difMx2(k,l,m)+(-3*(k-i)*(l-j)/
Q^5).*difMy2(k,l,m)...
                                +(-3*(k-i)*(m-s)/Q^5).*difMz2(k,l,m);

                            gradHms(i,j,s,1,3)=gradHms(i,j,s,1,3)+(1/
Q^3-3*(k-i)^2/Q^5).*difMx3(k,l,m)+(-3*(k-i)*(l-j)/
Q^5).*difMy3(k,l,m)...
                                +(-3*(k-i)*(m-s)/Q^5).*difMz3(k,l,m);

                            gradHms(i,j,s,1,4)=gradHms(i,j,s,1,4)+(1/
Q^3-3*(k-i)^2/Q^5).*difMx4(k,l,m)+(-3*(k-i)*(l-j)/
Q^5).*difMy4(k,l,m)...
                                +(-3*(k-i)*(m-s)/Q^5).*difMz4(k,l,m);

                            gradHms(i,j,s,1,5)=gradHms(i,j,s,1,5)+(1/
Q^3-3*(k-i)^2/Q^5).*difMx5(k,l,m)+(-3*(k-i)*(l-j)/
Q^5).*difMy5(k,l,m)...
                                +(-3*(k-i)*(m-s)/Q^5).*difMz5(k,l,m);

                            gradHms(i,j,s,1,6)=gradHms(i,j,s,1,6)+(1/
Q^3-3*(k-i)^2/Q^5).*difMx6(k,l,m)+(-3*(k-i)*(l-j)/
Q^5).*difMy6(k,l,m)...
                                +(-3*(k-i)*(m-s)/Q^5).*difMz6(k,l,m);

                            %gradient van H_ms_y naar de 6 volumes
                            gradHms(i,j,s,2,1)=gradHms(i,j,s,2,1)+
(-3*(k-i)*(l-j)/Q^5).*difMx1(k,l,m)+(1/Q^3-3*(l-j)^2/
Q^5).*difMy1(k,l,m)...
                                +(-3*(l-j)*(m-s)/Q^5).*difMz1(k,l,m);

                            gradHms(i,j,s,2,2)=gradHms(i,j,s,2,2)+
(-3*(k-i)*(l-j)/Q^5).*difMx2(k,l,m)+(1/Q^3-3*(l-j)^2/
Q^5).*difMy2(k,l,m)...
                                +(-3*(l-j)*(m-s)/Q^5).*difMz2(k,l,m);

                            gradHms(i,j,s,2,3)=gradHms(i,j,s,2,3)+
(-3*(k-i)*(l-j)/Q^5).*difMx3(k,l,m)+(1/Q^3-3*(l-j)^2/
Q^5).*difMy3(k,l,m)...
                                +(-3*(l-j)*(m-s)/Q^5).*difMz3(k,l,m);

                            gradHms(i,j,s,2,4)=gradHms(i,j,s,2,4)+
(-3*(k-i)*(l-j)/Q^5).*difMx4(k,l,m)+(1/Q^3-3*(l-j)^2/
Q^5).*difMy4(k,l,m)...
                                +(-3*(l-j)*(m-s)/Q^5).*difMz4(k,l,m);

                            gradHms(i,j,s,2,5)=gradHms(i,j,s,2,5)+
(-3*(k-i)*(l-j)/Q^5).*difMx5(k,l,m)+(1/Q^3-3*(l-j)^2/
Q^5).*difMy5(k,l,m)...
                                +(-3*(l-j)*(m-s)/Q^5).*difMz5(k,l,m);

                            gradHms(i,j,s,2,6)=gradHms(i,j,s,2,6)+
(-3*(k-i)*(l-j)/Q^5).*difMx6(k,l,m)+(1/Q^3-3*(l-j)^2/
Q^5).*difMy6(k,l,m)...
                                +(-3*(l-j)*(m-s)/Q^5).*difMz6(k,l,m);

                            %gradient van H_ms_z naar de 6 volumes
                            gradHms(i,j,s,3,1)=gradHms(i,j,s,3,1)+
(-3*(k-i)*(m-s)/Q^5).*difMx1(k,l,m)+(-3*(l-j)*(m-s)/
Q^5).*difMy1(k,l,m)...
                                +(1/Q^3-3*(m-s)^2/Q^5).*difMz1(k,l,m);

                            gradHms(i,j,s,3,2)=gradHms(i,j,s,3,2)+
(-3*(k-i)*(m-s)/Q^5).*difMx2(k,l,m)+(-3*(l-j)*(m-s)/
Q^5).*difMy2(k,l,m)...
                                +(1/Q^3-3*(m-s)^2/Q^5).*difMz2(k,l,m);

                            gradHms(i,j,s,3,3)=gradHms(i,j,s,3,3)+
(-3*(k-i)*(m-s)/Q^5).*difMx3(k,l,m)+(-3*(l-j)*(m-s)/
Q^5).*difMy3(k,l,m)...
                                +(1/Q^3-3*(m-s)^2/Q^5).*difMz3(k,l,m);

                            gradHms(i,j,s,3,4)=gradHms(i,j,s,3,4)+
(-3*(k-i)*(m-s)/Q^5).*difMx4(k,l,m)+(-3*(l-j)*(m-s)/
Q^5).*difMy4(k,l,m)...
                                +(1/Q^3-3*(m-s)^2/Q^5).*difMz4(k,l,m);

                            gradHms(i,j,s,3,5)=gradHms(i,j,s,3,5)+
(-3*(k-i)*(m-s)/Q^5).*difMx5(k,l,m)+(-3*(l-j)*(m-s)/
Q^5).*difMy5(k,l,m)...
                                +(1/Q^3-3*(m-s)^2/Q^5).*difMz5(k,l,m);

                            gradHms(i,j,s,3,6)=gradHms(i,j,s,3,6)+
(-3*(k-i)*(m-s)/Q^5).*difMx6(k,l,m)+(-3*(l-j)*(m-s)/
Q^5).*difMy6(k,l,m)...
                                +(1/Q^3-3*(m-s)^2/Q^5).*difMz6(k,l,m);

                        end
                    end
                end
            end
            Hmsnk(i,j,s)=norm((-1/4/pi)*Hms(i,j,s))^2;
        end
    end

end
gradHms=-1/4/pi*gradHms;
Hms=-1/4/pi*Hms;

Re: urgent: for loops

Pete — Thu, 01 Nov 2007 02:12:48 -0400

Wow.
That's a pretty amazing loop.

I'm not nearly clever enough to understand what it is doing, and I don't have
enough of it to test anything (and I'm confident I wouldn't understand it even
then), but I looked at it in awe.

Here are some general thoughts:-

1 - if you haven't already, search for 'matlab vectorization' in google or the
help files to get some advice on writing efficient MATLAB code. Also check
out 'colon' in MATLAB help.

2 - it's easier to see the similarities between lines if you don't break up each
line. I know MATLAB tries to restrict you to about 80 characters per line, but
in this case long lines helped me see that many are almost the same.

3 - you do the same calculations quite a few times. You could, for example,
calculate
temp = (1/Q^3-3*(k-i)^2/Q^5);
once and then use 'temp' every time you need that value.

4 - check out the 'image profiler' option if you want to find out exactly what
lines of your code are slow, and that can help decide what to fix.

-------------------------------------------------

I made an attempt at shortening your code myself.
Instead of using difMx1, difMx2, difMx3, difMx4, difMx5, difMx6, use difMx
and add an extra dimension at the end so that,
difMx(:,:,:,1) = difMx1(:,:,:);
difMx(:,:,:,2) = difMx2(:,:,:);
difMx(:,:,:,3) = difMx3(:,:,:);
difMx(:,:,:,4) = difMx4(:,:,:);
difMx(:,:,:,5) = difMx5(:,:,:);
difMx(:,:,:,6) = difMx6(:,:,:);

This means you can do 6 calculations on a single line using the colon, and the
whole of your inner loop can be shortened to

                            Q=sqrt((k-i)^2+(l-j)^2+(m-s)^2);

                            Hms(i,j,s,1)=Hms(i,j,s,1)+(M(k,l,m,1)/Q^3-3*(k-i)/Q^5*(M
(k,l,m,1)*(k-i)+M(k,l,m,2)*(l-j)+M(k,l,m,3)*(m-s)));
                            Hms(i,j,s,2)=Hms(i,j,s,2)+(M(k,l,m,2)/Q^3-3*(l-j)/Q^5*(M
(k,l,m,1)*(k-i)+M(k,l,m,2)*(l-j)+M(k,l,m,3)*(m-s)));
                            Hms(i,j,s,3)=Hms(i,j,s,3)+(M(k,l,m,3)/Q^3-3*(m-s)/Q^5*(M
(k,l,m,1)*(k-i)+M(k,l,m,2)*(l-j)+M(k,l,m,3)*(m-s)));

                            %gradient van H_ms_x naar de 6 volumes
                            gradHms(i,j,s,1,1)=sum(sum(sum((1/Q^3-3*(k-i)^2/
Q^5).*difMx(k,l,m,1)+(-3*(k-i)*(l-j)/Q^5).*difMy(k,l,m,1)+(-3*(k-i)*(m-s)/
Q^5).*difMz(k,l,m,1))));
                            gradHms(i,j,s,1,2:end)=gradHms(i,j,s,1,2:end)+(1/Q^3-3*
(k-i)^2/Q^5).*difMx(k,l,m,2:end)+(-3*(k-i)*(l-j)/Q^5).*difMy(k,l,m,2:end)+
(-3*(k-i)*(m-s)/Q^5).*difMz(k,l,m,2:end);

                            %gradient van H_ms_y naar de 6 volumes
                            gradHms(i,j,s,2,:)=gradHms(i,j,s,2,:)+(-3*(k-i)*(l-j)/
Q^5).*difMx(k,l,m,:)+(1/Q^3-3*(l-j)^2/Q^5).*difMy(k,l,m,:)+(-3*(l-j)*(m-s)/
Q^5).*difMz(k,l,m,:);

                            %gradient van H_ms_z naar de 6 volumes
                            gradHms(i,j,s,3,:)=gradHms(i,j,s,3,:)+(-3*(k-i)*(m-s)/
Q^5).*difMx(k,l,m,:)+(-3*(l-j)*(m-s)/Q^5).*difMy(k,l,m,:)+(1/Q^3-3*(m-s)^2/
Q^5).*difMz(k,l,m,:);

Actually, this is only a start - it could be shortened quite a lot more than this.
But if you want to try it, let me know if this works and then post your updated
code and perhaps it will be easier to work with.

(Note, I'm not a MATLAB genius either so it is quite possible that what I am
telling you is nonsense. But I think it is ok. Certainly I am trying to give a
sensible answer, but you will need to check the results of my code are the
same as yours because I can't test it.)

Re: urgent: for loops

Pete — Thu, 01 Nov 2007 03:39:52 -0400

Oh, I should have added, you'd also need to change difMy and difMz in the same
way as difMx.

Re: urgent: for loops

Pete — Thu, 01 Nov 2007 04:03:16 -0400

Oops.
Don't check out 'image profiler' (point 3 above) - that won't help you at all.
What I meant was just 'profiler' in the 'Desktop' menu, I think.
Sorry....

Re: urgent: for loops

adinda — Mon, 05 Nov 2007 13:25:56 -0500

On 1 nov, 05:03, "Pete " <pete.dot.bankh...@btinternet.dot.com> wrote:
> Oops.
> Don't check out 'image profiler' (point 3 above) - that won't help you at all.
> What I meant was just 'profiler' in the 'Desktop' menu, I think.
> Sorry....

Thank you!