算法优化 for循环问题（已搜索过版内帖）

ice-huya

大家都知道，matlab的for循环计算起来超级慢，尤其是像我现在遇到的需要10000000次循环计算的问题，该怎么解决？而且：
1.计算前不能确定向量的维数，不能进行预分配内存

2.for中有if、continue这样的语句，我不知道能不能转成向量化计算？该怎么转化？

我把程序帖在下面，大家帮忙看看，不胜感激！有什么建议都可以给点，能优化一点也十分感谢！我已被困扰好久，拜托各位了！！

for i=1:cyc-1   %cyc是个上千万的数，所以计算都要几天几夜，还常常算了几十个小时后发现内存不足   
   if (Smax(i)>=Smaxol(i));
      Smaxol(i+1)=Smax(i);   
       ol=i;                                             
   else
      Smaxol(i+1)=Smax(ol);
   end

   R(i)=Smin(i)/Smax(i);
   Y=sqrt((1.0/cos(PI*a(i)/w)));
   Yol=sqrt((1.0/cos(PI*a(ol)/w)));
   Kmax(i)=sqrt(PI*re*(1.0/cos(PI*Smax(i)/2.0/sigmav)+1.0))*(1.0+Y*sqrt(a(i)/2.0/re))*Smax(i);
   Kmaxol(i)=sqrt(PI*re*(1.0/cos(PI*Smaxol(i+1)/2.0/sigmav)+1.0))*(1.0+Yol*sqrt(a(ol)/2.0/re))*Smaxol(i+1);
   deltaK(i)=(1-R(i))*Kmax(i);

   alpha_1=1.0/(1.0-2.0*mu)+(1.0-1.0/(1.0-2.0*mu))/(1.0+0.8861*(t/(Kmax(i)/sigmay)^2.0)^3.2251)^0.75952;
   AA0=(0.825-0.34*alpha_1+0.05*alpha_1^2.0)*(cos(PI*Smax(i)/2.0/sigmaf))^(1.0/alpha_1);
   AA1=(0.415-0.071*alpha_1)*Smax(i)/sigmaf;
   AA3=2.0*AA0+AA1-1.0;
   AA2=1.0-AA0-AA1-AA3;

   if 0<=R(i) && R(i)<1.0;
      fop=max(R(i),AA0+AA1*R(i)+AA2*R(i)^2.0+AA3*R(i)^3.0);
   else
      fop=AA0+AA1*R(i);
   end

   if R(i)>=-2.0 && R(i)<0.0;
      feff(i)=(0.52-0.1*R(i))/(1.0-R(i));
         fth(i)=(1.0-R(i))^0.5;
   end

  if R(i)>=0.0 && R(i)<0.5;
      feff(i)=0.52+0.42*R(i)+0.06*R(i)^2.0;
       fth(i)=(1-R(i))^0.3;
  end

   if R(i)>=0.5 && R(i)<1.0;
      feff(i)=0.52+0.42*R(i)+0.06*R(i)^2.0;
     fth(i)=(1.05-1.4*R(i)+0.6*R(i)^2.0)^0.3;
   end
   

   deltaKeffth(i)=feff(i)*fth(i)*deltaKth0;
   deltaKu=Y*sqrt(PI*a(i))*(Smin(i)-Smin(i+1));
   alpha1=(1.0-1.65*mu)^2.0/5.0-1.0/(20.0*n0)*(1.0-1.65*mu)^(2.0/n0)+(1.0/PI-1.0/(2.2*n0)*(1.0/PI)^(1.0/n0)-((1.0-1.65*mu)^2.0/5.0-1.0/(20.0*n0)*((1.0-1.65*mu)^2.0)^(1.0/n0)))/(1.0+   (t/(Kmax(i)/sigmay)^2.0)/(1.0+1.0/n0))^(1.6+1.0/n0);
   alpha2=(1.0-1.65*mu)^2.0/5.0-1.0/(20.0*n0)*(1.0-1.65*mu)^(2.0/n0)+(1.0/PI-1.0/(2.2*n0)*(1.0/PI)^(1.0/n0)-((1.0-1.65*mu)^2.0/5.0-1.0/(20.0*n0)*((1.0-1.65*mu)^2.0)^(1.0/n0)))/(1.0+(t/(Kmaxol(i)/sigmay)^2.0)/(1.0+1.0/n0))^(1.6+1.0/n0);
   ry(i)=alpha1*(Kmax(i)/sigmay)^2.0;
   rol(i)=alpha2*(Kmaxol(i)/sigmay)^2.0;
   deltar=alpha1*(deltaKu/sigmay)^2.0;

   lamda=(1.0-1.65*mu)^2.0/5.0-1.0/(20.0*n0)*(1.0-1.65*mu)^(2.0/n0)+(1.0/PI-1.0/(2.2*n0)*(1.0/PI)^(1.0/n0)-((1.0-1.65*mu)^2.0/5.0-1.0/(20.0*n0)*((1.0-1.65*mu)^2.0)^(1.0/n0)))/(1.0+(t/(Kmax(i)/sigmay)^2.0)/(1.0+1.0/n0))^(1.6+1.0/n0);
   Kc=(((1.0-2.0*mu)^2.0-sqrt(1.0-mu^2.0))*PI*lamda/((1.0-2.0*mu)^2.0-1.0)/(1.0-2.0*mu)^2.0+(sqrt(1.0-mu^2.0)-1.0)/((1.0-2.0*mu)^2.0-1.0))*KIc;


   if a(i)+ry(i)<a(ol)+rol(i)-deltar;
      mf=((a(ol)+rol(i)-ry(i)-deltar)/a(i))^c;
      deltaKeff(i)=Kmax(i)*(1-mf*fop);
   else
      mf=1.0;
      deltaKeff(i)=Kmax(i)*(1-mf*fop);
      Smaxol(i+1)=Smax(i);
      ol=i;
   end


   if deltaKeff(i)>deltaKeffth(i);
      if Kmax(i)<Kc;
        da(i)=AA*(deltaKeff(i)-deltaKeffth(i))^mm/(1.0-(Kmax(i)/Kc)^nn);
       a(i+1)=a(i)+da(i);
       N(i)=i;
    else                             
        Nf=N(i-1);        af=a(i);
        result=[Nf;af;Kmax(i)];
        save result.dat result -ascii;
        break;   
    end
  else                  
       da(i)=0;
       a(i+1)=a(i)+da(i);
       continue;  
  end  

end

ice-huya

可能程序贴的有点乱，再补充说明一下我的疑问：

1. Smax(i)、Smin(i)是两个n*1的向量，维数相等，n的数量级在10000000左右，cyc是Smax的长度
2. 循环中除了有i在不断增加外，还有个ol是在某些条件下随i变化，我不知道这样的话，是不是就不能用向量化运算？
3. 最后的一个if中，又有continue，是不是也不能用向量计算代替了？

程序贴的有点乱了，给大家带来不便，还是希望能给我点指导，非常非常感谢！

pasuka

ice-huya 发表于 2011-12-22 14:54
可能程序贴的有点乱，再补充说明一下我的疑问：

1. Smax(i)、Smin(i)是两个n*1的向量，维数相等，n的数量 ...

先用profile，看看程序热点在哪里，然后有针对性的修改代码
又，啥注释也没有，理论公式也没有的一段冗长代码，我很怀疑是否有人愿意无偿替lz优化代码

rocwoods

      基本的注释没有，基本的用于测试的数据没有。还有一点需要说明的是MATLAB的for循环计算起来慢，那是很久之前的事。对于高版本的MATLAB，for循环一点也不慢。楼主可以测试下，MATLAB中10亿次空循环的时间，和C++中的对比。然后看下10亿次1+1这样的操作和C++的对比，看看差距明显吗？再看看一次计算sin(1:1e6)和C++中循环计算，是不是超越C++了？MATLAB for循环表现的慢其实本身不是MATLAB循环机制的慢，而是函数解析的慢。
    在MATLAB中，built-in函数是调用效率最高的，单次调用的开销和C++的相当，但是普通的m文件相对来说就慢了不少，匿名函数和普通文件差不多，稍慢些。至于inline这样的结构，单次调用开销是built-in函数的上万倍。
    所以很多时候，MATLAB表现的慢是浪费在了函数开销以及内存不断分配上了，真正的正事没干多少。这就好比从北京平谷到天津蓟县，本来几十公里的路程，非得坐大巴几小时到首都机场，然后等几小时飞机，然后10分钟飞到天津，然后再几个小时到蓟县。一共一天的时间，就10分钟花在正事上。
    向量化的本质就是“摊薄”函数调用的开销，使得本身“有用计算”的开销在单次函数调用中明显大于函数解析的开销，只要抓住这个原则，MATLAB的程序不会慢的。
    用profile看看你程序瓶颈在哪里，看看哪些函数调用次数最多，然后看看那个函数能否向量化，而不是一味避免循环和continue。
    循环1千万次不算多，我做过的项目，单幅图像将近25G那么大，要逐一判断每个像素点(近100亿个)是否在一个15万边的复杂多边形内。如果在MATLAB中用最笨的逐一循环去做，粗略估算了下时间要30年左右。而合理利用循环和向量化，外加算法角度的优化，七、八个小时即可完成。

ice-huya

rocwoods 发表于 2011-12-22 22:41
基本的注释没有，基本的用于测试的数据没有。还有一点需要说明的是MATLAB的for循环计算起来慢，那是 ...

非常感谢版主的回复，又给我了希望:loveliness:

我会试一下profile再进行重点优化。

这个程序是算裂纹扩展的，我怕代码带长贴不下才把能省的注释省去了。敢问一下，如果可以的话，能不能给我传一份您的循环次数很多的、但优化后的程序，让我学习一下。我的邮箱ice-huya@163.com

再次感谢！还希望以后多和您交流。

ice-huya

pasuka 发表于 2011-12-22 19:37
先用profile，看看程序热点在哪里，然后有针对性的修改代码
又，啥注释也没有，理论公式也没有的一段冗长 ...

谢谢你！:)

我先试试看profile，如果还是解决不了，再贴上详细的注释请教大家