Double summation with vectorized loops in Matlab

I want to vectorize this double for loop because it is a bottleneck in my code. Since Matlab is a based-one indexing language I have to create an additional term for M = 0.

R,r,lambda are constants

Slm(L,M),Clm(L,M) are matrices 70x70

Plm(L,M) is a matrix 70x71

Cl(L),Pl(L) are vectors 70x1

% function dU_r
  s1 = 0;
  for L = 2:70
     s1 = s1 + ((R/r)^L)*(L+1)*Pl(L)*Cl(L);
     for m = 1:L
        s1 = s1 + ((R/r)^L)*(L+1)*Plm(L,M)*(Clm(L,M)*...
cos(M*lambda) + Slm(L,M)*sin(M*lambda));
     end
  end

  dU_r = -(mu_p/(r^2))*s1;


  % function dU_phi

  s2=0;
  for L = 2:70
      s2 = s2 + ((R/r)^L))*Plm(L,1)*Cl(L);
          for m = 1:l
          s2 = s2 + ((R/r)^L)*(Plm(L,M+1)-M*tan(phi)*Plm(L,M))*...
  (Clm(L,开发者_开发技巧M)*cos(M*lambda) + Slm(L,M)*sin(M*lambda));
          end;
  end;
  dU_phi = (mu_p/r)*s2;

  % function dU_lambda

  s3=0;
      for L=2:70
          for m=1:L
          s3 = s3 + ((R/r)^L)*M*Plm(L,M)*(Slm(L,M)*cos(M*lambda)...
              - Clm(L,M)*sin(M*lambda));
          end;
  dU_lambda = (mu_p/r)*s3;

I think the following code may be a part of the solution (only partially vectorized), for a full vectorization you might want to look at meshgrid, ndgrid, bsxfun and/or repmat, I suppose.

s2     = 0;
L      = 2:70;
PlCl   = Pl.*Cl;
PlmClm = Plm.*Clm;
Rrl    = (R/r).^(L).*(L+1);
COS    = cos((1:70)*lambda);
SIN    = sin((1:70)*lambda);

for l=L
    M  = 1:l;
    s1 = sum(Rrl(l-1)*PlmClm(l,M).*(COS(M) + Slm(l,M).*SIN(M)));
    s2 = s2 + s1;
end;
s2 = s2 + sum(Rrl(L).*PlCl(L));

I haven't tried to run this code, so if it may complain about some dimensions. You should be able to figure that out (just some transposes here or there may do).

Just a note on the side: try to keep away from short variable names (and certainly ambiguous ones like l (that can be misread for a 1, a I or an l). Also, it might be easier to vectorize your code if you have actual (i.e. not coded) expressions to start with.

edit: applied corections suggested by gnovice

For this particular situation, it will likely be more efficient for you to not fully vectorize your solution. As Egon illustrates, it's possible to replace the inner loop through vectorization, but replacing the outer for loop as well could potentially slow-down your solution.

The reason? If you pay attention to how you index the matrices Plm, Clm, and Slm in your loop, you only ever use values from the lower triangular parts of them. All of the values above the main diagonal are ignored. By fully vectorizing your solution, such that it uses matrix operations instead of a loop, you would end up performing unnecessary multiplications with the zeroed-out upper triangular portions of the matrices. This is one of those cases where a for loop may be a better option.

As such, here's a refined and corrected version of Egon's answer that should be close to optimum with regard to the number of mathematical operations performed:

index = 2:70;
coeff = ((R/r).^index).*(index+1);
lambda = lambda.*(1:70).';  %'
cosVector = cos(lambda);
sinVector = sin(lambda);
s2 = coeff*(Pl(index).*Cl(index));
for L = index
  M = 1:L;
  s2 = s2 + coeff(L-1)*((Plm(L,M).*Clm(L,M))*cosVector(M) + ...
                        (Plm(L,M).*Slm(L,M))*sinVector(M));
end

Taking like reference a solution given by Jan Simon 1, I have written the following code for my problem

Rr = R/r;
RrL = Rr;  
cosLambda = cos((1:70)* lambda);
sinLambda = sin((1:70)* lambda);
u1 = uint8(1);
s1 = 0;
s2 = 0;
s3 = 0;
for j = uint8(2):uint8(70)
   RrL = RrL * Rr;
   q1 = RrL * (double(j) + 1);
   t1 = Pl(j) * datos.Cl(j);
   q2 = RrL;
   t2 = Plm(j,1) * Cl(j);
   t3 = 0;
   for m = u1:j
      t1 = t1 + Plm(j,m) * ...
         (Clm(j, m) * cosLambda(m) + ...
          Slm(j, m) * sinLambda(m));
        t2 = t2 + (Plm(j,m+1)-double(m)*tan_phi*Plm(j,m))*...
    (Clm(j,m)*cosLambda(m) + Slm(j,m)*sinLambda(m));
        t3 = t3 + double(m)*Plm(j,m)*(Slm(j,m)*cosLambda(m)...
                - Clm(j,m)*sinLambda(m));
     end
     s1 = s1 + q1 * t1;
     s2 = s2 + q2 * t2;
     s3 = s3 + q2 * t3;
  end
dU_r = -(mu_p/(r^2))*s1;
dU_phi = (mu_p/r)*s2;
dU_lambda = (mu_p/r)*s3