Will someone please check the following code for efficiency?

I'm working on a rather large matlab project, and I have one function called DoShade that is a huge bottleneck. I've run profile on it, but that didn't give me much useful advice, except what I already knew, which is that the while loop is running some 250,000 times and is taking up most of the time. Even a 0.01s speedup in that while loop will yield huge gains overall. Any advice will be greatly appreciated. Things that don't work:

• Replacing NormSolVector * VectorToOrigin' with dot(,)
• Vectorizing dx, dy, idx, or jdy
• Replacing -realmax with a variable

The typical size of Dem is around 500x500.

function sombra = DoShade(Dem, SolVector, dl)
  %Calculate inverse solar vector
  InvSolVector = -SolVector./max(abs(SolVector(1:2)));
  %Initialize and calculate Normal Solar Vector
  NormSolVector(1:3)=0;
  NormSolVector(3) = sqrt(SolVector(1)^2+SolVector(2)^2);
  NormSolVector(1) = -SolVector(1)*SolVector(3)/NormSolVector(3);
  NormSolVector(2) = -SolVector(2)*SolVector(3)/NormSolVector(3);
  cols = size(Dem,2);
  rows = size(Dem,1);
  %If SolVector(x) is negative, sun is on the west
  if SolVector(1)<=0
    f_i=1;
  else
    f_i=rows; %Default
  end
  %If SolVector(y) is negative, sun is on the north
  if SolVector(2)<=0
    f_j=1;
  else
    f_j=cols; %Default
  end
  %Initialize the shaded variable
  sombra = ones(size(Dem));
  i=f_i;
  for j=1:cols
    m=0;
    zcompare = -realmax;%Initial value lower than any possible projection
    while (1)
      dx=InvSolVector(1)*m;
      dy=InvSolVector(2)*m;
      idx=round(i+dx);
      jdy=round(j+dy);
      if ((idx<1) || (idx>rows) || (jdy<1) || (jdy>cols))
        break;
      end
      VectorToOrigin(1)=dx*dl;
      VectorToOrigin(2)=dy*dl;
      VectorToOrigin(3)=Dem(idx,jdy);
      zprojection = NormSolVector * VectorToOrigin';
      if (zprojection < zcompare)
        sombra(idx,jdy)=0;
      else
        zcompare = zprojection;
      end
      m=m+1;
%     if m>100000
%       error('Limit hit')
%     end
    end
  end
  j=f_j;
  for i=1:rows
    n=0;
    zcompare = -realmax;%Initial value lower than any possible projection
    while (1)
      dx=InvSolVector(1)*n;
      dy=InvSolVector(2)*n;
      idx=round(i+dx);
      jdy=round(j+dy);
      if ((idx<1) || (idx>rows) || (jdy<1) || (jdy>cols))
        break;
      end
      VectorToOrigin(1)=dx*dl;
      VectorToOrigin(2)=dy*dl;
      VectorToOrigin(3)=Dem(idx,jdy);
      zprojection = NormSolVector * VectorToOrigin';
      if (zprojection < zcompare)
        sombra(idx,jdy)=0;
      else
        zcompare = zprojection;
      end
      n=n+1;
%     if n>100000
%       error('Limit hit')
%     end
    end
  end
end

[dpb]

Posting the output of the profiler within the function would likely be helpful since we have no data.

Answer 1

[the cyclist]

I see that you have preallocated memory for sombra.

But does this line

sombra(idx,jdy)=0;

ever grow the size of the array? In other words, are idx or jdy ever larger than the current dimensions of the array? If so, you'll want to preallocate to the largest size sombra can become.

[Hayden Smotherman]

No, the if statement breaks the while loop if idx or jdy ever grow beyond the size of the Dem, which is the same size as sombra is preallocated to.

Answer 2

[Peter O]

A few things to try:

How often does Dem change? Is it common for it to be exactly 496x488, or do its dimensions always change?

If it has the same size on multiple consecutive loops, what happens if you make sombra persistent but sensitive to the size of Dem? You're reallocating ~1MB of memory each time you call. So, in this manner it would only reallocate when needed. I'm not sure if the persistent overhead will cancel it out. Something like:

persistent sombra sz_lastDem

% Other Code Here

sz_Dem = size(Dem)
if isequal(sz_Dem,sz_lastDem)
  sombra(:) = 1;
else 
  sombra = ones(sz_Dem);
  sz_lastDem = sz_Dem;
end

Two other alternatives:

1. Pass multiple Dems, SolVectors, and dl's to it at once by adding a dimension to each and iterating inside the function. You'll eliminate some function call overhead, and if you use the above trick on sombra you don't need to declare them persistent (relaxing memory checks on function entry/exit)

2. This is going to sound gross and go against clean coding practice, but can you inline DoShade? I ask because I had a similar function problem regarding a call to compute a Jacobian ~30k times, and by inlining it and using the above trick on selective allocation of the larger matrix I got around a 30% speedup. All from removing function call overhead and reducing the memory allocation requirements.

[Hayden Smotherman]

Dem never changes size between function calls, at least as far as DoShade is concerned. Neither does dl. Dem and dl are factors changed in a main script, but are pretty much constant. It's the solarvector that is changing every iteration, if that matters.

As for adding a dimension and iterating side of the function, I wanted to do that, but was unsure exactly how I could add a dimension to something like sombra, which is already a matrix. Some structure definition?

EDIT- How can I define sombra to be persistent if it's a function output? Matlab is rightly giving me an error saying I'm trying to define a persistent variable after it has been used. This is when I define "persistent sombra" immediately after the function header.