Optimizing a Vectorized Matlab Function

Question

When i run profiler it tell me that the most time consuming code is the function vdist. Its a program that measures distance between two points on earth considering earth as an ellipsoid. The code looks standard and i don't know where and how it can be improved upon. The initial comments say, it has already been vectorized. Is there a counterpart to it in some other language which can be used as a MEX file. All i want is improvement in terms of time efficiency. Here is a link to the code from Matlab FEX.

http://www.mathworks.com/matlabcentral/fileexchange/8607-vectorized-geodetic-distance-and-azimuth-on-the-wgs84-earth-ellipsoid/content/vdist.m

The function is called from within a loop as- (You can find the function as its the most time consuming line here)

              109 for i=1:polySize 
              110    % find the two vectors needed
       11755  111    if i~=1 
0.02   11503  112        if i<polySize 
0.02   11251  113         p0=Polygon(i,:); p1=Polygon(i-1,:); p2=Polygon(i+1,:);    
         252  114        else 
         252  115         p0=Polygon(i,:); p1=Polygon(i-1,:); p2=Polygon(1,:); %special case for i=polySize 
         252  116        end 
         252  117    else 
         252  118         p0=Polygon(i,:); p1=Polygon(polySize,:); p2=Polygon(i+1,:); %special case for i=1 
         252  119    end 
 0.02  11755  120    Vector1=(p0-p1); Vector2=(p0-p2); 
 0.06  11755  121    if ~(isequal(Vector1,Vector2) || isequal(Vector1,ZeroVec) || isequal(Vector2,ZeroVec)); 
              122        %determine normals and normalise and
 0.17  11755  123        NV1=rotateVector(Vector1, pi./2); NV2=rotateVector(Vector2, -pi./2); 
 0.21  11755  124        NormV1=normaliseVector(NV1); NormV2=normaliseVector(NV2); 
              125        %determine rotation by means of the atan2 (because sign matters!)
       11755  126        totalRotation = vectorAngle(NormV2, NormV1); % Bestimme den Winkel totalRotation zwischen den normierten Vektoren 
       11755  127      if totalRotation<10 
       11755  128          totalRotation=totalRotation*50; 
       11755  129      end 
0.01   11755  130      for res=1:6 
0.07   70530  131         U_neu=p0+NV1; 
17.01  70530  132         [pos,a12] = vdist(p0(:,2),p0(:,1),U_neu(:,2),U_neu(:,1)); 
0.02   70530  133         a12=a12+1/6.*res*totalRotation; 
       70530  134         ddist=1852*safety_distance; 
4.88   70530  135         [lat2,lon2] = vreckon(p0(:,2),p0(:,1),ddist, a12); 
0.15   70530  136         extendedPoly(f,:)=[lon2,lat2];f=f+1; 
< 0.01 70530  137      end 
       11755  138    end 
       11755  139 end

What are the lines within that function that are acting as a bottleneck? Also can you show us an example of how you call this function. If you're calling it from within a loop and with scalars then you're not going to get any benefit from vectorized code. — slayton
Make sure you give vectorized input, it will likely not be easy to improve the general performance of this algorithm but perhaps you can cut out some pieces that you don't use? — Dennis Jaheruddin
@slayton- i've included the calling loop above. As far as the bottleneck goes, its not like a couple of lines are taking all the time, but its kinda distributed. Yeah, some lines are taking more time than others, but the highest for a line is 10% of the whole function. — Vikram

High Performance Mark High Performance Mark · Accepted Answer · 2012-10-26T09:38:29

No matter how hard I study the code that's been posted, I can't see why the call to vdist is made inside the loop.

When I'm trying to optimise a block of code inside a loop one of the things I look for are statements which are invariant, that is which are the same at each call, and which can therefore be lifted out of the loop.

Looking at

130      for res=1:6 
131         U_neu=p0+NV1; 
132         [pos,a12] = vdist(p0(:,2),p0(:,1),U_neu(:,2),U_neu(:,1)); 
133         a12=a12+1/6.*res*totalRotation; 
134         ddist=1852*safety_distance; 
135         [lat2,lon2] = vreckon(p0(:,2),p0(:,1),ddist, a12); 
136         extendedPoly(f,:)=[lon2,lat2];f=f+1; 
137      end

I see

in l131 the variables p0, NV1 appear only on the rhs, and they only appear on the rhs elsewhere inside the loop, so this statement is loop-invariant and can be lifted out of the loop; only a small time saving perhaps;
in l134 again, I see another loop-invariant statement, which can again be lifted out of the loop for another small time saving;
but then I started to look very closely, and I can't see why l132, where the call to vdist is made, is inside the loop either. None of the values on the rhs of that assignment are modified in the loop (other than U_neu but I've already lifted that out of the loop).

Tidying up what was left a bit, this is what I ended up with:

U_neu=p0+NV1; 
[pos,a12] = vdist(p0(:,2),p0(:,1),U_neu(:,2),U_neu(:,1)); 
ddist=1852*safety_distance; 
for res=1:6 
   extendedPoly(f,:) = vreckon(p0(:,2),p0(:,1),ddist, a12+1/6.*res*totalRotation); 
   f=f+1; 
end

Optimizing a Vectorized Matlab Function

2 Answers