嵌套循环优化

Question

我真的很难在 MATLAB 上优化微积分代码。

获得非线性材料的材料属性需要进行繁重的计算。

此计算需要超过 2.4 亿步。 它本身相当简单，因为它包含大量sum。 唯一的问题是数字存储在各种数组和列表中，这有点令人困惑。

这里是原始代码：

Tensor=zeros(3,3,3,3);
for m=1:3
    for n=1:3
        for o=1:3
            for p=1:3
                for x=1:16 
                    for y=1:16
                        for z=1:16 
                            for i=1:3
                                for j=1:3
                                    for k=1:3
                                        for l=1:3
                                            for r=1:3
                                                for s=1:3
                                                    sum=sum+(1/(8*(pi()^2))*P{x,y,z}(i,m)*P{x,y,z}(j,n)*P{x,y,z}(k,o)*P{x,y,z}(l,p)*(TensorC(i,j,k,l)-TensorC0(i,j,r,s))*TensorA(r,s,k,l)*sin(omega(x))*p_omega(x)*p_phi(y)*p_beta(z);
                                                end
                                            end
                                        end                                                    
                                    end
                                end
                            end
                        end
                    end
                end
                Tensor(m,n,o,p)=sum;
            end
        end
    end
end

P{x,y,z}(i,m) 是基础公式的变化（其他相同）：i and m 确定公式的类型并使用x、y 和z 变量计算结果。

求和中的所有其他数字都是在数组和张量中选取的实数。

我试图从最后一个 for 循环中提取尽可能多的变量，以便尽快计算它们并减少操作次数：

Tensor=zeros(3,3,3,3);
CO1=1/(8*(pi()^2));
for m=1:3
    for n=1:3
        for o=1:3
            for p=1:3
            sum=C0_tensor(m,n,o,p);
                for x=1:16
                    CO7=sin(omega(x));
                    CO8=p_omega(x);
                    for y=1:16
                        CO9=p_phi(y);
                        for z=1:16
                            CO10=p_beta(z);
                            for i=1:3
                                CO2=P{x,y,z}(i,m);
                                for j=1:3
                                    CO3=P{x,y,z}(j,n);
                                    for k=1:3
                                        CO4=P{x,y,z}(k,o);
                                        for l=1:3
                                            CO5=P{x,y,z}(l,p);
                                            CO6=TensorC(i,j,k,l);
                                            for r=1:3
                                                for s=1:3
                                                    CO11=TensorC0(i,j,r,s);
                                                    CO12=TensorA(r,s,k,l);
                                                    sum=sum+CO1*CO2*CO3*CO4*CO5*(CO6-CO11)*CO12*CO7*CO8*CO9*CO10;
                                                end
                                            end
                                        end                                                    
                                    end
                                end
                            end
                        end
                    end
                end
                Tensor(m,n,o,p)=sum;
            end
        end
    end
end

不过，计算时间太长了。

我没有看到任何并行化或矢量化计算的方法......

从一个数组或一个矩阵中检索一个特定值的操作似乎很慢......

你认为我应该构建一个包含所有值的巨大张量而不是使用倍数吗？

Answer 1

你不应该使用sum作为变量名，因为你用相同的名字覆盖了有用的函数。

这里只考虑内部循环，您正在为每个 r 和 s 计算一个值，并将其添加到您的输出值中：

  for r=1:3
         for s=1:3
              CO11=TensorC0(i,j,r,s);
              CO12=TensorA(r,s,k,l);
              sum=sum+CO1*CO2*CO3*CO4*CO5*(CO6-CO11)*CO12*CO7*CO8*CO9*CO10;
         end
  end

但是，您可以立即将 C011/C012 作为 3 x 3 矩阵，求和，然后将其添加到您的输出中： （在此处将您的 sum 更改为 out，请注意 .* 而不是 * 在适当的位置）：

C011 = squeeze(TensorCO(i,j,:,:));
C012 = squeeze(TensorCO(:,:,k,l));
s = CO1*CO2*CO3*CO4*CO5*(CO6-CO11).*CO12*CO7*CO8*CO9*CO10;
out = out + sum(s(:));

另外，当你这样做时：

CO7=sin(omega(x));
CO8=p_omega(x);

（在后面的等式中就是 C07*C08）——你不需要为每个 n,m,p. 重新计算 sin(omega(x))，这样就可以完全退出循环。

预计算sin 并乘以p_omega（循环外）：

omega78 = sin(omega).*p_omega;

然后只需在 x 循环中检索值 C78 = omega78(x) 并使用它而不是 C07*C08。