Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: urgent: for loops
Date: Thu, 1 Nov 2007 02:12:48 +0000 (UTC)
Organization: Queen's University
Lines: 77
Message-ID: <fgbcn0$jap$1@fred.mathworks.com>
References: <1193749892.996365.194660@o38g2000hse.googlegroups.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1193883168 19801 172.30.248.38 (1 Nov 2007 02:12:48 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 1 Nov 2007 02:12:48 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 943450
Xref: news.mathworks.com comp.soft-sys.matlab:435543



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.)