Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Vectorizing code

Subject: Vectorizing code

From: Sven

Date: 25 Oct, 2010 08:31:06

Message: 1 of 3

Dear all,

I am struggling a bit with a piece of my code. This actual piece seems to be the bottleneck in my code (according to profiler). This is mostly due to the for/if loops.
I was wondering if it is possible to vectorize some part (or completely) of it.

The code:

  
for i=1:ngrids-1
        sumk = zeros(ngrids,1);
        for k = 1:i
            Vdif = Vhalf(i+1)-V(k);
            if Vdif < Vhalf(2)
                Vdif = Vhalf(2);
            end
        
            for j=2:length(Vhalf)
                if Vdif <= Vhalf(j) && Vdif > Vhalf(j-1)
                alpha=j;
                end
            end

            int1 = zeros(ngrids,1);
            for j = alpha:ngrids
                int1(j)=f(j)/V(i)*(Vhalf(j+1)-Vhalf(j))*agglom((i*2),(k*2));
                    
            end
            
                int2=f(alpha-1)/(Vhalf(alpha)+Vhalf(i+1)-V(k))*(Vhalf(alpha)-Vhalf(i+1)-V(k))*agglom((((2*alpha)-1)+((2*i)+1)-(k*2))/2);
                
                sumk (k) = (sum(int1)+int2)*(Vhalf(k+1)-Vhalf(k))*f(k);
        end
        Fphalf (i) = sum(sumk);
end

If anyone of you could give me some pointers it would be much appreciated.

Subject: Vectorizing code

From: Roger Stafford

Date: 25 Oct, 2010 14:14:04

Message: 2 of 3

"Sven " <sven.molenaar.prive@gmail.com> wrote in message <ia3f8a$od2$1@fred.mathworks.com>...
> Dear all,
>
> I am struggling a bit with a piece of my code. This actual piece seems to be the bottleneck in my code (according to profiler). This is mostly due to the for/if loops.
> I was wondering if it is possible to vectorize some part (or completely) of it.
>
> The code:
>
> for i=1:ngrids-1
> sumk = zeros(ngrids,1);
> for k = 1:i
> Vdif = Vhalf(i+1)-V(k);
> if Vdif < Vhalf(2)
> Vdif = Vhalf(2);
> end
>
> for j=2:length(Vhalf)
> if Vdif <= Vhalf(j) && Vdif > Vhalf(j-1)
> alpha=j;
> end
> end
>
> int1 = zeros(ngrids,1);
> for j = alpha:ngrids
> int1(j)=f(j)/V(i)*(Vhalf(j+1)-Vhalf(j))*agglom((i*2),(k*2));
> end
>
> int2=f(alpha-1)/(Vhalf(alpha)+Vhalf(i+1)-V(k))*(Vhalf(alpha)-Vhalf(i+1)-V(k))*agglom((((2*alpha)-1)+((2*i)+1)-(k*2))/2);
> sumk (k) = (sum(int1)+int2)*(Vhalf(k+1)-Vhalf(k))*f(k);
> end
> Fphalf (i) = sum(sumk);
> end
>
> If anyone of you could give me some pointers it would be much appreciated.
- - - - - - - - - - -
    The two innermost for-loops involving j are presumably crucial to the efficiency of your algorithm. I see a few ways these could be improved.

  As for the first of these, the way it is written you are seeking the last index j for which the given inequalities hold. You should therefore be going backward with j starting from j = length(Vhalf) and moving downwards towards j = 2 and seeking the first j for which the inequalities hold and then doing a 'break'. You can economize further by seeking (while still moving backwards) the first j for which Vdif <= Vhalf(j) and breaking. Then continuing with that j, seek the first j for which Vdif > Vhalf(j-1) and again breaking. You will have arrived at alpha. This should cut down the expected length of array you have to cover and the number of inequalities that need calculating while doing so.

  The second j-for-loop can be greatly improved. In the subsequent code you never use int1 for any other purpose than taking its sum. To avoid adding up the same quantities in this for-loop over and over again, you can establish a cumulative sum of the quantities f(j)*(Vhalf(j+1)-Vhalf(j)) but working backwards from ngrids. That is,

 c = -cumsum(f(ngrids:-1:2).*diff(Vhalf(ngrids+1:-1:2));

can be computed just once before ever entering into any of the for-loops. Then instead of this second for-loop you would just write

 sumint1 = c(ngrids+1-alpha)*agglom((i*2),(k*2))/V(i);
 
which you would use in place of the sum(int1) term in computing sumk(k). There would be no need to allocate space for an int1 array.

  A third improvement might be to avoid the repeated sumk allocation. Just keep a running sum of sumk which you can place in Fphalf(i) when you exit from the "for k = 1:i" loop.

  Your calculation of Vdif can be simplified to

 Vdif = max(Vhalf(i+1)-V(k),Vhalf2);

though this may be only a saving in the number of lines of code.

  Anyway this is a starter for you. Perhaps you will be inspired to find many other ways of improving this code. Note that it isn't necessarily the for-loops that are the source of the excessive computation time but what you do within these loops. Vectorization per se does not necessarily save computation time. It is the algorithm that counts.

Roger Stafford
 

Subject: Vectorizing code

From: Sven

Date: 27 Oct, 2010 13:55:05

Message: 3 of 3

Dear Roger,

Thanks for the starter, never really learned matlab programming in this way. Perhaps thats my own fault. I implemented the suggestions you made and I made a huge improvement. You were also right it has shown me some new possibilities for the rest of my code.

Thanks alot.

Regards,
Sven

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ia43bc$9d5$1@fred.mathworks.com>...
> "Sven " <sven.molenaar.prive@gmail.com> wrote in message <ia3f8a$od2$1@fred.mathworks.com>...
> > Dear all,
> >
> > I am struggling a bit with a piece of my code. This actual piece seems to be the bottleneck in my code (according to profiler). This is mostly due to the for/if loops.
> > I was wondering if it is possible to vectorize some part (or completely) of it.
> >
> > The code:
> >
> > for i=1:ngrids-1
> > sumk = zeros(ngrids,1);
> > for k = 1:i
> > Vdif = Vhalf(i+1)-V(k);
> > if Vdif < Vhalf(2)
> > Vdif = Vhalf(2);
> > end
> >
> > for j=2:length(Vhalf)
> > if Vdif <= Vhalf(j) && Vdif > Vhalf(j-1)
> > alpha=j;
> > end
> > end
> >
> > int1 = zeros(ngrids,1);
> > for j = alpha:ngrids
> > int1(j)=f(j)/V(i)*(Vhalf(j+1)-Vhalf(j))*agglom((i*2),(k*2));
> > end
> >
> > int2=f(alpha-1)/(Vhalf(alpha)+Vhalf(i+1)-V(k))*(Vhalf(alpha)-Vhalf(i+1)-V(k))*agglom((((2*alpha)-1)+((2*i)+1)-(k*2))/2);
> > sumk (k) = (sum(int1)+int2)*(Vhalf(k+1)-Vhalf(k))*f(k);
> > end
> > Fphalf (i) = sum(sumk);
> > end
> >
> > If anyone of you could give me some pointers it would be much appreciated.
> - - - - - - - - - - -
> The two innermost for-loops involving j are presumably crucial to the efficiency of your algorithm. I see a few ways these could be improved.
>
> As for the first of these, the way it is written you are seeking the last index j for which the given inequalities hold. You should therefore be going backward with j starting from j = length(Vhalf) and moving downwards towards j = 2 and seeking the first j for which the inequalities hold and then doing a 'break'. You can economize further by seeking (while still moving backwards) the first j for which Vdif <= Vhalf(j) and breaking. Then continuing with that j, seek the first j for which Vdif > Vhalf(j-1) and again breaking. You will have arrived at alpha. This should cut down the expected length of array you have to cover and the number of inequalities that need calculating while doing so.
>
> The second j-for-loop can be greatly improved. In the subsequent code you never use int1 for any other purpose than taking its sum. To avoid adding up the same quantities in this for-loop over and over again, you can establish a cumulative sum of the quantities f(j)*(Vhalf(j+1)-Vhalf(j)) but working backwards from ngrids. That is,
>
> c = -cumsum(f(ngrids:-1:2).*diff(Vhalf(ngrids+1:-1:2));
>
> can be computed just once before ever entering into any of the for-loops. Then instead of this second for-loop you would just write
>
> sumint1 = c(ngrids+1-alpha)*agglom((i*2),(k*2))/V(i);
>
> which you would use in place of the sum(int1) term in computing sumk(k). There would be no need to allocate space for an int1 array.
>
> A third improvement might be to avoid the repeated sumk allocation. Just keep a running sum of sumk which you can place in Fphalf(i) when you exit from the "for k = 1:i" loop.
>
> Your calculation of Vdif can be simplified to
>
> Vdif = max(Vhalf(i+1)-V(k),Vhalf2);
>
> though this may be only a saving in the number of lines of code.
>
> Anyway this is a starter for you. Perhaps you will be inspired to find many other ways of improving this code. Note that it isn't necessarily the for-loops that are the source of the excessive computation time but what you do within these loops. Vectorization per se does not necessarily save computation time. It is the algorithm that counts.
>
> Roger Stafford
>

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us