Speeding up code: vectorization and others

Question

Alex Kurek 2016-8-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/300444-speeding-up-code-vectorization-and-others

编辑： per isakson 2016-8-28

Hi,

Is there a way to speed this up even more? I did what I could already: everything is precomputed, preallocated etc. and compiled to .mex.

First, just to give you an idea about sizes of continers:

bessels = ones(1201, 1201, 101);    % 1.09 GB
negMcontainer = ones(1201, 1201, 100);
posMcontainer = negMcontainer;
maxN = 100;
levels = maxN + 1;
xElements = 1201;
Aj1 = complex(ones(101, 101);
Aj2 = Aj1;

Code:

parfor i = 1 : xElements
    for j = 1 : xElements
        umn = complex(zeros(levels, levels)); % cleaning
        for n = 0:maxN
            mm = 1;
            for m = -n:2:n
                nn = n + 1; % for indexing
                if m < 0
                    umn(nn, mm) = bessels(i, j, nn) * negMcontainer(i, j, abs(m));
                end
                if m > 0
                    umn(nn, mm) = bessels(i, j, nn) * posMcontainer(i, j, m);
                end
                if m == 0
                    umn(nn, mm) = bessels(i, j, nn);
                end
                mm = mm + 1; % for indexing
            end % m
        end % n
        beta1 = sum(sum(Aj1.*umn));
        betaSumSq1(i, j) = abs(beta1).^2;
        beta2 = sum(sum(Aj2.*umn));
        betaSumSq2(i, j) = abs(beta2).^2;
  end % j
end % i

Best regards, Alex

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Alex Kurek 2016-8-24

编辑：Alex Kurek 2016-8-24

Nobody knows? So far everything I try is not speeding this up.

This if prof Profiler:

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

per isakson 2016-8-24

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/300444-speeding-up-code-vectorization-and-others#answer_232676

编辑：per isakson 2016-8-24

在 MATLAB Online 中打开

Your code choked my computer. R2016a, no parallel toolbox, no mex. However, regarding the innermost loop

for m = -n:2:n
end

Move nn = n + 1; outside the loop. Helps a little bit (tic,toc).
Replace the three if-statements by a if-elseif-elseif-else-end. Again small effect. However, 30% faster when profiling.
Split the loop into two loops one with m<0 and one with m>0 and remove the if-statements. Note: odd and even n. Again small effect. However, when profiling the elapse time was cut in half.

Conclusion: The "JIT/Accelerator" don't need my help in these cases. However, the "JIT/Accelerator" is less effective together with the profiler.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Alex Kurek 2016-8-24

编辑：Alex Kurek 2016-8-24

在 MATLAB Online 中打开

Ok, now it looks like this:

parfor j = 1 : xElements
    for i = 1 : xElements
        umn = complex(zeros(levels, levels)); % cleaning 
        besselIJ = bessels(i, j, :);
        for n = 0:maxN
            nn = n + 1; % for indexing
            currentBessel = besselIJ(nn);
            mm = 1;
            for m = -n:2:n
                  if m < 0
                      umn(nn, mm) = currentBessel * negMcontainer(i, j, abs(m));
                  elseif m == 0
                      umn(nn, mm) = currentBessel;
                  else
                      umn(nn, mm) = currentBessel * posMcontainer(i, j, m) ;
                  end
                  mm = mm + 1;
              end % m
          end % n
          beta1 = sum(sum(Aj1.*umn));
          betaSumSq1(i, j) = abs(beta1).^2;
          beta2 = sum(sum(Aj2.*umn));
          betaSumSq2(i, j) = abs(beta2).^2;
      end % i
  end % j

I have 2016a, mex and parallel and only moving nn = n+1 outside the loop really helped. What helped more, is changing the code co columnwise (first two lines).

One small question: if I split the loop into two to remove if-statements, haw can I efficiently solve the odd/even problem?

Best regards and thanx, Alex

Alex Kurek 2016-8-25

在 MATLAB Online 中打开

Thanx again. I just coded it:

bessels = ones(1201, 1201, 101);    % 1.09 GB
negMcontainer = ones(1201, 1201, 100);
posMcontainer = negMcontainer;
maxN = 100;
levels = maxN + 1;
xElements = 1201;
Aj1 = complex(ones(101, 101));
Aj2 = Aj1;
betaSumSq1 = zeros(xElements, xElements);
betaSumSq2 = betaSumSq1;
for j = 1 : xElements
    umn = complex(zeros(levels, levels)); % cleaning 
    if rem(j, 25) == 0
        disp(j);
    end
    for i = 1 : xElements
          for n = 0 : 2 : maxN
              nn = n + 1;
              mm = 1;
              for m = - n : 2 : -2
                  umn(nn, mm) = bessels(i, j, nn) * negMcontainer(i, j, abs(m));
                  mm = mm + 1;
              end
              umn(nn, mm) = bessels(i, j, nn);
              for m = 2 : 2 : n
                  umn(nn, mm) = bessels(i, j, nn) * posMcontainer(i, j, m) ;
                  mm = mm + 1;
              end
          end
          for n = 1 : 2 : maxN
              nn = n + 1;
              mm = 1;
              for m = - n : 2 : -1
                  umn(nn, mm) = bessels(i, j, nn) * negMcontainer(i, j, abs(m));
                  mm = mm + 1;
              end
              m = 1:2:n;
              numOfEl = ceil(n/2);
              umn(nn, mm:mm+numOfEl-1) = bessels(i, j, nn) * posMcontainer(i, j, m);
  %             for m = 1 : 2 : n
  %                 umn(nn, mm) = bessels(i, j, nn) * posMcontainer(i, j, m);
  %                 mm = mm + 1;
  %             end
          end
          beta1 = sum(sum(Aj1.*umn));
          betaSumSq1(i, j) = abs(beta1).^2;
          beta2 = sum(sum(Aj2.*umn));
          betaSumSq2(i, j) = abs(beta2).^2;
      end % i
  end % j

The speed in .mex and parfor is exectly the same, as before. But it opens a possibility of vectorization. I vectorized last loop (as pasted aboce) and it slowed down a bit. Have I done something wrong?

Best regards, Alex

per isakson 2016-8-27

编辑：per isakson 2016-8-28

I don't know what to say. I would have guessed that your vectorization would be faster than the loop. However, loops are much faster in current releases of Matlab than they were when we learned that loops to should always be avoided.

Lesson learned

It's difficult to be smarter than the "JIT/Acceletator"
The function, profile, cannot always be trusted

Alex Kurek 2016-8-28

在 MATLAB Online 中打开

Thanx. It is possible to speed-up vectorized version by ~17% by doing this:

tic
for j = 1 : xElements
    for i = 1 : xElements
        for n = 1 : 2 : maxN
            nn = n + 1;
            numOfEl = ceil(n/2);
            umn2(nn, 1:numOfEl) = bessels(i, j, nn) * posMcontainer(i, j, 1:2:n);
        end
    end
end
toc % 1.059022 seconds