在Matlab中使用矢量化循环进行双加总

我想将此double for循环矢量化，因为它是我的代码中的瓶颈。由于Matlab是基于酮索引语言我要创建为M = 0。在Matlab中使用矢量化循环进行双加总

R，R，λ-附加项是常数

犿（L，M），CLM（L，M）是矩阵70x70

PLM（L，M）是一个矩阵70x71

CL（L），PL（L）是矢量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;

来源

2011-05-30 julian

我认为下面的代码可以是在一个部分解决方案（仅部分矢量化），为一个完全矢量化你可能想看看meshgrid,ndgrid,bsxfun和/或repmat，我想。

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));

我还没有试过运行这段代码，所以如果它可能会抱怨一些尺寸。你应该能够弄清楚（只是在这里或那里可能会有一些转换）。

只需要注意一点：尽量避免使用短变量名称（当然也可以使用l（这可能会误解为1，I或l），这样可能会更容易矢量化你的代码，如果你有实际的（即不编码）表达式开始

编辑：应用通过gnovice

来源

2011-05-31 00:18:55 Egon

有几个错误：1）你需要索引'Rrl'用'1-1'的而不是环'l'。 2）乘以余弦项的系数和“Plm”也乘以正弦项（注意问题中的括号）。 – gnovice 2011-05-31 04:39:08

非常感谢您的帮助。干杯。 – julian 2011-05-31 08:44:30

@gnovice：我很感激你的代码的更正。干杯。 – julian 2011-05-31 08:45:44

建议对于这种特殊情况corections，它可能会更有效的为您不完全矢量化您的解决方案。As Egon illustrates，可以更换内部环路gh矢量化，但替换外部for循环以及可能潜在减速您的解决方案。

原因？如果您注意如何在循环中对Plm,Clm和Slm中的矩阵进行索引，那么您只能使用它们的lower triangular parts中的值。主对角线上方的所有值均被忽略。通过向量化您的解决方案，使其使用矩阵运算而不是循环，您将最终执行矩阵的零化上三角部分的不必要的乘法运算。这是for循环可能是更好的选择的情况之一。

因此，这里有Egon's answer精致而正确的版本，应该是接近最优关于执行的数学运算的数量：以相同的参考扬西蒙1给出解决办法

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

来源

2011-05-31 05:13:37 gnovice

非常感谢您的帮助。干杯。 – julian 2011-05-31 08:43:03

，我写了下面的代码，我的问题

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

来源

2011-06-01 18:37:40 julian

在Matlab中使用矢量化循环进行双加总

回答

相关问题