2013-12-12 77 views
4

我真的很努力在MATLAB上优化微积分代码。嵌套循环优化

获取非线性材料的材料属性是一项繁重的计算。

此计算需要超过240万步。 它本身相当简单,因为它包含了大量的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的变化(同为其他):i and m确定式的类型并且将结果与xyz变量计算。

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

我试图以尽快计算出它们的减少操作次数,提取从最后的循环尽可能多的变量:

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 

尽管如此,计算实在是太长了。

我看不出并行或矢量化计算的任何方式......

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

待办事项你认为我应该建立一个包含所有值而不是使用多个值的巨大张量?

+1

如果真的很难避免循环,你应该开始考虑使用其他有利于循环的编程语言,比如'C++' – Ray

+1

8o ...好亲切......并祝你好运。你对matlab代码有很大的依赖性吗? – Brian

+0

量化*太长*。 –

回答

1

由于您正在用相同的名称覆盖有用的函数,因此不应将sum用作变量名称。

以刚内环这里,你计算每个rs单个值被添加到您的输出值:

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 ×3个矩阵,做一个总和,然后添加到您的输出: (改变了你的sumout这里,请注意*,而不是*在适当的地方。):

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在以后的公式) - 你不需要((x)的欧米茄),每n,m,p.重新计算的罪,因此可以取出完全的循环。

预计算sin和由p_omega(外侧环)相乘:

omega78 = sin(omega).*p_omega; 

然后,只需检索在x环路和使用而不是C07*C08C78 = omega78(x)

+0

非常感谢您的帮助! 我运用了你的想法,并调整了一些数学和数组/张量结构。 基本上,主循环速度快了10倍。 我会继续努力寻找新的想法! – Empirehell

+1

我优化了张量和数组结构。 现在主循环运行9分钟而不是2小时! – Empirehell