编写了一个程序计算数据，发现一个点的计算时间要2分钟左右，我有五十万个点要计算，耗时太长了，求助大牛帮我优化下 下边的程序，减少下运行时间。
程序的主要功能是，输入一组数据到part[1][3],然后从三个0-90度独立变化的角a，b，c中选择一个组合使得O=misorientation(con[m].C,Mfcc,con[n].C,Mbcc)计算得到的数值最小，在计算每个abc组合的时候需要从 24*24个con[m].C 和 con[n].C的计算结果中选择最小的计算结果和别的abc组合进行比较。由于实际计算中需要很多组数据所以还有一定的 Osum=Osum+Omin; num++;等求平均值得操作，这些在计算单个数据的时候可以不考虑。在计算的时候abc和输入的数据都是转化为矩阵来计算的。


程序如下：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define Pi 3.14159  



void matproduct(float a2[3][3],float b2[3][3],float c2[3][3]);
float misorientation(float con1[3][3],float fcc1[3][3],float con2[3][3],float bcc1[3][3]);
void transpose(float a3[3][3]);
void fccorientation (float par[][3],float fc[3]);

int main()
{
 int i;
 float part[1][3]={2.9976,46.554,37.132},fcc[3];
 fccorientation(part,fcc);
 printf("%f %f %f\n",fcc[0],fcc[1],fcc[2]);
 return 0;
}

void fccorientation (float par[][3],float fc[3])
{float a1,b1,c1,a,b,c,o,p,q,Omin,Oaver,Osum,O;
 float Oavermin=Pi;                            //abc组合当中的取向差的最小值   
 float Ocritical=1;                           //直接输入数据不要输入公式 计算得到的取向差需要符合的外部条件，可以通过计算结果进行调整
 float Mfcc[3][3],Mbcc[3][3];
 int num,i,m,n,x;


 
 struct convert
 {float C[3][3];};
  struct convert con[24]=
    {
        {1,0,0,0,1,0,0,0,1},
        {0,0,-1,0,-1,0,-1,0,0},
        {0,0,-1,0,1,0,1,0,0},
        {-1,0,0,0,1,0,0,0,-1},
        {0,0,1,0,1,0,-1,0,0},
        {1,0,0,0,0,-1,0,1,0},
        {1,0,0,0,-1,0,0,0,-1},
        {1,0,0,0,0,1,0,-1,0},
        {0,-1,0,1,0,0,0,0,1},
        {-1,0,0,0,-1,0,0,0,1},
        {0,1,0,-1,0,0,0,0,1},
        {0,0,1,1,0,0,0,1,0},
        {0,1,0,0,0,1,1,0,0},
        {0,0,-1,-1,0,0,0,1,0},
        {0,-1,0,0,0,1,-1,0,0},
        {0,1,0,0,0,-1,-1,0,0},
        {0,0,-1,1,0,0,0,-1,0},
        {0,0,1,-1,0,0,0,-1,0},
        {0,-1,0,0,0,-1,1,0,0},
        {0,1,0,1,0,0,0,0,-1},
        {-1,0,0,0,0,1,0,1,0},
        {0,0,1,0,-1,0,1,0,0},
        {0,-1,0,-1,0,0,0,0,-1},
        {-1,0,0,0,0,-1,0,-1,0},

    };

 for (a=0.0;a<=90;a=a+1)
        for (b=0.0;b<=90;b=b+1)
            for (c=0.0;c<=90;c=c+1)
            {a1=a/180.0*Pi;
             b1=b/180.0*Pi;
             c1=c/180.0*Pi;
             Mfcc[0][0]=(cos(a1)*cos(c1)-sin(a1)*sin(c1)*cos(b1));
             Mfcc[0][1]=(sin(a1)*cos(c1)+cos(a1)*sin(c1)*cos(b1));
             Mfcc[0][2]=sin(c1)*sin(b1);
             Mfcc[1][0]=(-cos(a1)*sin(c1)-sin(a1)*cos(c1)*cos(b1));
             Mfcc[1][1]=(-sin(a1)*sin(c1)+cos(a1)*cos(c1)*cos(b1));
             Mfcc[1][2]=cos(c1)*sin(b1);
             Mfcc[2][0]=sin(a1)*sin(b1);
             Mfcc[2][1]=-cos(a1)*sin(b1);
             Mfcc[2][2]=cos(b1);


             Osum=0;
             num=0;
             for (i=0;i<1;i++)
             {
                 
              o=par[i][0];
              p=par[i][1];
              q=par[i][2];
              if ((o==0)&&(p==0)&&(q==0))
              {
               Omin=0;
              }                     //对应if条件语句
              else
              {Mbcc[0][0]=cos(o)*cos(q)-sin(o)*sin(q)*cos(p);      //由于初始化在程序的外边所以只能单独赋值
               Mbcc[0][1]=sin(o)*cos(q)+cos(o)*sin(q)*cos(p);
               Mbcc[0][2]=sin(q)*sin(p);
               Mbcc[1][0]=-cos(o)*sin(q)-sin(o)*cos(q)*cos(p);
               Mbcc[1][1]=-sin(o)*sin(q)+cos(o)*cos(q)*cos(p);
               Mbcc[1][2]=cos(q)*sin(p);
               Mbcc[2][0]=sin(o)*sin(p);
               Mbcc[2][1]=-cos(o)*sin(p);
               Mbcc[2][2]=cos(p);
              



               Omin=Pi;
               for (m=0;m<24;m++)
                for (n=0;n<24;n++)
                { O=misorientation(con[m].C,Mfcc,con[n].C,Mbcc);
               
                   if (O<Omin)
                       Omin=O;
                }                   //对应for条件语句
              }                     //对应else 条件语句


              Osum=Osum+Omin;
              num++;
             }                      //对应i层次的for 循环结束
             Oaver=Osum/num;

             if ((Oaver<=Oavermin)&&(Oaver<=Ocritical))   //在a，b，c的变化范围内寻找取向差最小的abc组合，同时最小的取向差还需要满足程序开始设定的标准Ocritical
             {Oavermin=Oaver;
              fc[0]=a;
              fc[1]=b;
              fc[2]=c;
             }
            }//对应最上边角度变化的for语句
            printf("%f",Oavermin);   
}//对应整个函数的结尾




float misorientation(float con1[3][3],float fcc1[3][3],float con2[3][3],float bcc1[3][3])
{int i;
 float O1;
 float vfb[3][3]={{0.7416,-0.6667,-0.0749},{0.6498,0.7416,-0.1667},{0.1667,0.0749,0.9832}}; //取向关系的矩阵描述
 float a1[3][3],b1[3][3],c1[3][3],d1[3][3];

 matproduct(vfb,con1,a1);

 matproduct(a1,fcc1,b1);
 
 matproduct(con2,bcc1,c1);

 transpose(c1);



 matproduct(b1,c1,d1);
 //for(i=0;i<3;i++)
     //printf("%f %f %f\n",d1[i][0],d1[i][1],d1[i][2]);
if((d1[0][0]+d1[1][1]+d1[2][2]-1)/2>=1)
O1=0;
else if((d1[0][0]+d1[1][1]+d1[2][2]-1)/2<=-1)
O1=3.14159;
else
 O1=acos((d1[0][0]+d1[1][1]+d1[2][2]-1)/2);

 //printf("%f",O1);

 return(O1);
}


void matproduct(float a2[3][3],float b2[3][3],float c2[3][3])      //对矩阵进行乘法运算，结果存储到c矩阵中
{int i1,j1;
 for (i1=0;i1<3;i1++)
     for (j1=0;j1<3;j1++)
         {c2[i1][j1]=a2[i1][0]*b2[0][j1]+a2[i1][1]*b2[1][j1]+a2[i1][2]*b2[2][j1];}
}


void transpose(float a3[3][3])//对矩阵进行转置运算
{float c3;
 c3=a3[0][1];
 a3[0][1]=a3[1][0];
 a3[1][0]=c3;
 c3=a3[0][2];
 a3[0][2]=a3[2][0];
 a3[2][0]=c3;
 c3=a3[1][2];
 a3[1][2]=a3[2][1];
 a3[2][1]=c3;
}

[ 本帖最后由 zhht87 于 2013-8-16 18:16 编辑 ]

你是要我们看完程序推测你这个程序的功能如何再想怎么优化？
不是我说，脑子呢？

我把程序的功能大概的介绍了一下

你的这个程序 是计算一个点的吗，   怎么一直没有结束运行， 都有十分钟了。

[ 本帖最后由 XiaoXiao_Ren 于 2013-8-16 23:21 编辑 ]

挺有趣的东西，但隔行如隔山啊，我的物理知识只有刚体力学和电磁学的基础知识，取向差，是物理化学领域的概念么？

从代码的角度分析，你的程序存在大量的冗余运算，计算一个点大概需要进行4百多亿次的浮点乘法运算和上千万次的三角运算。

关于三角运算，你实际计算的只是a b c o p q的正余弦，将它们计算后的值代入矩阵式里可以节省很多时间。

对于Mbcc，可以提前计算好再进入abc的穷举循环中参与运算，没必要也不应该放在循环里每次都重新计算一遍。现在你只写了一个点，将来扩充到五十万个点也是一样，这里可以节省很多时间。

关于矩阵运算部分，你的con的意义我不清楚，从我的角度看它们参与矩阵运算只是调整另一个矩阵的行列的顺序及取负。针对这一点可以优化矩阵运算过程，没必要用普通的矩阵运算过程。最终的计算结果只要矩阵的主对角线元素，这里也没必要进行完整的矩阵运算。

不清楚Mfcc的意义（为什么这么计算），如果你能解释清楚整个计算我想将有很大的优化空间。包括a b c这三个循环，这是以穷举的方式寻找目标极值。应该花点时间好好研究一下这个计算过程，0到90度这个区间的正余弦都是单调的，如果能找到某种收敛速度不错的迭代公式将大大缩短计算时间。

以你现在的计算方法计算五十万个点少则三两年多则数十年，如果能花个把月找到更好的算法将计算时间缩短到几天（我的直觉应该可以缩短到几分钟），那这是很值得做的。


如果这个程序对你很重要，那么请从基础讲一讲这个取向差和Mxcc、con的意义，我很有兴趣扩展一下自己的知识面，并帮你完成算法的分析构建。

你好，前几天我看到你回复的帖子说如果计算结果是a,b,c的连续函数的话，可以有方法快速求解，我今天计算了一下 发现在a确定的情况下，计算的结果是b，c的连续函数。但是返回去看你的帖子突然找不到了，请问你能详细的说下，这种连续函数如何求解计算结果的最小值吗？

不愧是楼主啊  你提的建议非常好，先计算三角运算并将部分冗余运算移出循环后计算三个点用时200s，但是计算50万个点也要不到200天。 所以还请你在算法上多帮我分析下，

这个程序主要是用来计算取向的，取向是材料学里边的概念，就是一个个晶体的特征轴的方向，原始数据是角度，计算时转化为矩阵计算，

通过实验可以得到bcc这种晶体的取向数据，而我想要求得fcc这种晶体又和bcc有一定的关系，所以知道bcc 的情况下 可以求出来对应的fcc，在程序中Omin就是理想的fcc 和给出的abc计算得到的fcc 两个取向之间的角度差。

但是复杂的地方在于晶体具有一定的对称性，所以不能直接就用矩阵相乘求解 还要考虑对称性对于数据的影响，就是你说的24个con 矩阵，由于涉及到两种晶体bcc fcc 而这两种晶体都有对称性，所以con矩阵要独立变化两次 就是程序中的m和n 变化的部分。


刚开始编程的时候没有仔细想就用了穷举加打擂台的方法计算角度差的最小值，但是经你提醒之后作图发现Oaver值在a给定的情况下在b和c的变化范围内是连续的，之前论坛上有个朋友说过连续的情况下 可以用迭代收敛的方法求最小值，但没有详细说明算法应该如何设计，所以想请教下版主这种情况下如何设计求最小值的算法。

呵呵，我还当你已经放弃这个问题了。

你的解释还是太简单，我无法根据这点信息建立起一个有效的分析模型。当然，我也理解，也许要解释清楚每一个概念基本相当于自学材料学了。

这样吧，我问两个我最不解的概念吧。

1.fcc这个矩阵的意义是什么，它是怎么得来的？这是我最感兴趣的地方

2.con为什么有24个，为什么不是48个？

我不知道你是怎么作图的，在不考虑你的晶体的对称性的因素的前提下，结果对于a b c的函数是连续的。我没实际分解，但预感这个函数甚至是单调的（在三个维度上）。

但是加入晶体的对称性后，我不确认它是否连续（呵呵，毕竟你这式子太长了，实在没有手工计算的勇气）。顺便说一句，虽然两个矩阵都独立地进行对称变换，但实际上最后的乘积只有最多48种可能，显然24*24次运算中也存在大量冗余。