程序优化求指导

system · 2013 年9 月 18 日 02:14

#include <stdlib.h>
 #include <stdio.h> 
 #include <windows.h>
 #include<iostream>
 #include <time.h>
 #include <math.h>
 #include <curand.h>
 #include <curand_kernel.h>
 #include <cuda_gl_interop.h>
 #include <cuda_runtime.h>
 using namespace std;
 #define FeaturesNum1 5
 #define THREDANUM1 32
 #define M_PI1       3.14159265358979323846
 #define PNum1 24
 #define PDim1 17
 __constant__  int pso_lanmark_GPU[FeaturesNum1][3];
 __constant__   double  Featuresinfo_GPU[FeaturesNum1][3];
 __constant__ double alpas_GPU[FeaturesNum1][3][10];
 __constant__ double pose_Xup_GPU[17];
 __constant__   double pose_Xdown_GPU[17];
 __constant__   double Vmax_GPU[17];

 __device__ double myrand1(unsigned int &x,unsigned int &y,unsigned int &z,unsigned int &c,int tidx)
 {  
 unsigned long long t, A = 698769069ULL;  
 x = 69069*x+1234+tidx*1000;  
 y ^= (y<<13+tidx); 
 y ^= (y>>17+tidx); 
 y ^= (y<<5+tidx);  
 t = (A*z + c);
 c = (t >> 32);
 z = t;
 return (double)((x+y+z)/RAND_MAX)/RAND_MAX/4;  
 }

 __global__ void ParticleFly_GPU(double* CudaDat ,double* data_out,int isize,int circulate)
 {

 __shared__ double Fit[PNum1];
 __shared__ int GBestIndex ;
 int tidx = threadIdx.x;  
   unsigned int x; 
   unsigned int y; 
   unsigned int z;
   unsigned int c;  
   x=123456789;
   y=362436000;
   z=521288629;
   c=7654321;
 __shared__ double test[PNum1];
 __shared__ double Wmax ;
 __shared__ double Wmin ;
 __shared__ double W;
 __shared__ double C1;
 __shared__ double C2;
 __shared__ double optical_center[2];
 Wmax = 0.9;
 Wmin = 0.4;
 W = 1;
 C1 = 2;
 C2 = 2;
  optical_center[0] = 172.766;
 optical_center[1] = 306.3046;

   for(int k = 0;k<circulate;k++)
   {
   if(tidx < PNum1)
   {
   if( k <= 1500)
   W = Wmax - k*(Wmax-Wmin)/circulate; 
   else
   W = Wmin;
   for(int j=0; j<PDim1; j++)   //合并访问优化的问题 
   {
   *(CudaDat + (PDim1 + j)*PNum1 + tidx)  = W*(*(CudaDat + (PDim1 + j)*PNum1 + tidx)) +//修改速度
   myrand1(x,y,z,c,tidx)*C1*(*(CudaDat + (2*PDim1 + j)*PNum1 + tidx) -(*(CudaDat + j*PNum1 + tidx)) )+
   myrand1(x,y,z,c,tidx)*C2*(*(CudaDat + (2*PDim1 + j)*PNum1 + GBestIndex) -(*(CudaDat + j*PNum1 + tidx)));
   
   if(*(CudaDat + (PDim1 + j)*PNum1 + tidx) > Vmax_GPU[j]) *(CudaDat + (PDim1 + j)*PNum1 + tidx) = Vmax_GPU[j];
   
   if(CudaDat[(PDim1 + j)*PNum1 + tidx]<-Vmax_GPU[j]) *(CudaDat + (PDim1 + j)*PNum1 + tidx) = -Vmax_GPU[j];
   *(CudaDat +  j*PNum1 + tidx) += (*(CudaDat + (PDim1 + j)*PNum1 + tidx)); //修改坐标
   
   if(CudaDat[j*PNum1 + tidx]>pose_Xup_GPU[j]) 
   CudaDat[j*PNum1 + tidx]=pose_Xup_GPU[j];//保护
   
   if(CudaDat[j*PNum1 + tidx]<pose_Xdown_GPU[j]) 
   CudaDat[j*PNum1 + tidx]=pose_Xdown_GPU[j];

   }
   for(int i = 0;i<17;i++)
   test = CudaDat[i*PNum1];

 //////////////////////////////////////GetFit////////////////////////////////
   double cost = 0;   
   double angle_cos[3],angle_sin[3];

   angle_cos[0] = cos(*(CudaDat + 3*PNum1 + tidx));
   angle_cos[1] = cos(*(CudaDat + 4*PNum1 + tidx));
   angle_cos[2] = cos(*(CudaDat + 5*PNum1 + tidx));
   angle_sin[0] = sin(*(CudaDat + 3*PNum1 + tidx));
   angle_sin[1] = sin(*(CudaDat + 4*PNum1 + tidx));
   angle_sin[2] = sin(*(CudaDat + 5*PNum1 + tidx));

   test[1] = *(CudaDat + 3*PNum1 + tidx);
   for(int i=0; i< FeaturesNum1; i++)
   {
   double alpasx = 0,alpasy =0,alpasz=0;
   for(int j=0; j<10; j++)
   {
   alpasx += *(CudaDat + (j+7)*PNum1 + tidx)*alpas_GPU[0][j];
   alpasy += *(CudaDat + (j+7)*PNum1 + tidx)*alpas_GPU[1][j];
   alpasz += *(CudaDat + (j+7)*PNum1 + tidx)*alpas_GPU[2][j];
   
   }
   double xx = Featuresinfo_GPU[0] + alpasx;
   double xy = Featuresinfo_GPU[1] + alpasy;
   double xz = Featuresinfo_GPU[2] + alpasz;
   /*double xx = Featuresinfo[0];
   double xy = Featuresinfo[1];
   double xz = Featuresinfo[2];

   double wx = xx*angle_cos[2]*angle_cos[1] + xy*(angle_cos[2]*angle_sin[1]*angle_sin[0]-angle_sin[2]*angle_cos[0]) 
   + xz*(angle_cos[2]*angle_sin[1]*angle_cos[0]+angle_sin[2]*angle_sin[0]) +*(CudaDat+ tidx);
   double wy = xx*angle_sin[2]*angle_cos[1] + xy*(angle_sin[2]*angle_sin[1]*angle_sin[0]+angle_cos[2]*angle_cos[0]) 
   + xz*(angle_sin[2]*angle_sin[1]*angle_cos[0]-angle_cos[2]*angle_sin[0]) + *(CudaDat + 1*PNum1 + tidx);
   double wz = -xx*angle_sin[1] + xy*angle_cos[1]*angle_sin[0] +xz*angle_cos[1]*angle_cos[0] + *(CudaDat + 2*PNum1 + tidx);
   double px = optical_center[0] + *(CudaDat + 6*PNum1 + tidx)*(wx/wz) - (double)pso_lanmark_GPU[0];
   double py = optical_center[1] - *(CudaDat + 6*PNum1 + tidx)*(wy/wz) - (double)pso_lanmark_GPU[1];
   cost = cost +  px*px + py*py; 
   //cost = cost + abs(px)+abs(py);
   }
   Fit[tidx] = -1*cost;  //Fit得到全部一样的值

   if( Fit[tidx] >= *(CudaDat + (3*PDim1+1)*PNum1 + tidx))
   {
   *(CudaDat + (3*PDim1+1)*PNum1 + tidx) = Fit[tidx];
   for(int j=0; j<PDim1; j++)
   *(CudaDat + (2*PDim1+j)*PNum1 + tidx) = *(CudaDat + j*PNum1 + tidx);
   }
   }
   __syncthreads();
   if(tidx < PNum1)
   {
   if(tidx == 0)  //用第一个线程计算GBestIndex
   {
   GBestIndex = 0;
   for(int i=0; i<PNum1; i++)
   if(*(CudaDat + (3*PDim1+1)*PNum1 + i) >= *(CudaDat + (3*PDim1+1)*PNum1 + GBestIndex) && i!=GBestIndex) GBestIndex = i;
   }
   }
   }
   __syncthreads();
   if(tidx == 0 )
   {
   for(int i = 0;i<isize;i++)
   {
   //data_out = CudaDat[GBestIndex*isize + i];
   data_out = *(CudaDat + i*PNum1 + GBestIndex);
   }
   }
 }

 extern "C" 
 clock_t ParticleFly_CPU(double* CudaDat,double* data_out,int isize,int circulate,int pso_lanmark1[5][3])
 {
   
   double Featuresinfo1[FeaturesNum1][3]={{-31.1956,34.2044,94.8308},{30.3130,34.2285,94.5579},{0.0072,10.1444,128.9015},{-29.3909,-33.8801, 98.1061},{29.1029,-33.5093,97.8227}};
   double alpas1[FeaturesNum1][3][10]=.......................................

   double pose_Xup1[] = {   450, 450,1000, M_PI1/2, M_PI1/2, M_PI1/2, 900, 0.002653,  0.0016677, 0.0013104, 0.0009396, 0.0008288, 0.0006270, 0.0005781, 0.0005434, 0.0005332, 0.0004848}; //自变量上界    
   double pose_Xdown1[] = {-450,-450,1000,-M_PI1/2,-M_PI1/2,-M_PI1/2,-900,-0.002653, -0.0016677,-0.0013104,-0.0009396,-0.0008288,-0.0006270,-0.0005781,-0.0005434,-0.0005332,-0.0004848}; //自变量下界  
   double Vmax[PDim1];

   for(int i=0; i<PDim1; i++)
   Vmax = (pose_Xup1-pose_Xdown1)*0.1;
   

 double* CudaDat_GPU;
 double* data_out_GPU;
   clock_t start = clock();
   
 cudaMalloc((void**) &CudaDat_GPU,sizeof(double)*PNum1*isize);
 cudaMalloc((void**) &data_out_GPU,sizeof(double)*isize);
 cudaMemset(CudaDat_GPU, 0, sizeof(double)*PNum1*isize);
 cudaMemset(data_out, 0, sizeof(double)*isize);
 clock_t start1 = clock();

 cudaMemcpy(CudaDat_GPU,CudaDat,sizeof(double)*PNum1*isize,cudaMemcpyHostToDevice);

 cudaMemcpyToSymbol( pso_lanmark_GPU, pso_lanmark1, 
   sizeof(int) *5*3); 
 cudaMemcpyToSymbol( Featuresinfo_GPU, Featuresinfo1, 
   sizeof(double) * 5*3); 
 cudaMemcpyToSymbol( alpas_GPU, alpas1, 
   sizeof(int) *5* 3*10); 
 cudaMemcpyToSymbol( pose_Xup_GPU, pose_Xup1, 
   sizeof(double) * 17); 
 cudaMemcpyToSymbol( pose_Xdown_GPU, pose_Xdown1, 
   sizeof(double) * 17); 
 cudaMemcpyToSymbol( Vmax_GPU, Vmax, 
   sizeof(double) * 17); 

   clock_t end1 = clock();

 ParticleFly_GPU<<<1,THREDANUM1>>>(CudaDat_GPU,data_out_GPU,isize,circulate);
 cudaThreadSynchronize();
 clock_t end2 = clock();
 cudaMemcpy(data_out,data_out_GPU,sizeof(double)*isize,cudaMemcpyDeviceToHost);
   clock_t end = clock();

 cudaFree(CudaDat_GPU);
 cudaFree(data_out_GPU);

 cout<<start1 - start<<"  -------  "<<end1 - start1<<"------"<<end2-end1<<"------"<<end-end2<<endl;
 return end - start;
 }

最近在对一个PSO算法进行进行CUDA编程优化，我的配置是pentium dual coreE5300 2.6FHZ，GTX650TI显卡，这个CPU函数调用四次，通过clock()函数测量时间，得到cudaMalloc 部分第一次时间消耗70 毫秒，核函数时间消耗200毫秒左右，对于传入数据，我已经做了一些合并访问的修改，性能提高了2倍左右（以前核函数400毫秒）。
我想问的是-------
对与这个函数，还有哪些地方可以改进的地方，我觉得这个运行速度还有很大的提升空间，自己刚刚学校CUDA不久，跪求指点一二！！！！！

system · 2013 年9 月 18 日 02:24

LZ您好：

论坛仅供讨论CUDA问题，而不能直接为您的科学问题/工程问题的具体实现直接给出代码实现/调试/优化等工作。

因而您直接一句“PSO算法”如何，然后贴几百行几乎无说明的代码，最后来一个“求优化”，这样做是不合适的，请您理解这一点。

如需继续讨论，请您将您的问题抽象为数学问题形式（即不能假设看帖的人和您一样有关于该学科的背景知识），并介绍您当前实现的规划和意图，当前使用的数据结构，当前的线程组织形式和规模，profiler结果如何，哪部分有哪些影响效能的问题有待解决。

如果仅是现在这样干贴几百行代码的话，斑竹无法为您目测出问题。

祝您好运~

system · 2013 年9 月 18 日 02:29

以及大量代码请使用“代码模式”贴出，您这么多代码直接贴出，显示拥挤，缩进混乱，观感极为不友好。

此外，注意到您调用kernel的一句“ParticleFly_GPU<<<1,THREDANUM1>>>(CudaDat_GPU,data_out_GPU,isize,circulate);”您这样只使用一个block的话，将会极大地影响您的算法在GPU上的效能，请您注意这一点。

祝您好运~

system · 2013 年9 月 18 日 02:30

以及，device端计算的耗时请以event计时为准，或者使用profiler统计的时间为准。

system · 2013 年9 月 18 日 02:49

谢谢您的提出的问题，我在最下面做了一些简单的介绍，可能没有说清楚，这是一个PSO算法的基本实现，程序中采用24个点，每个点17维（核函数一个块，每个块32个线程，每个线程处理一个点的运算，每个点需要循环3000次，所有在每个线程里边循环3000次），循环3000，直接上确实有点马虎了。作为新人，还不会使用profiler，真是不好意思

system · 2013 年9 月 18 日 02:50

一看我这编程规范，主机都觉得太囧，请大家原谅，主机以后注意

system · 2013 年9 月 18 日 02:51

楼主提醒得是，我马上整理一下自己的代码

system · 2013 年9 月 18 日 02:54

LZ您好：

您在1#最下方给出的信息和您在5#补充的叙述都不足以告知读者足够清楚的信息，以便分析您的问题的。

以及“核函数一个块，每个块32个线程”，如前所述，这样将极大地影响您的GPU的计算效能，请勿将一个CUDA core看做一个CPU core来直接套用CPU上的实现方法。

祝您好运~

system · 2013 年9 月 18 日 03:03

多谢版主的回复，可能我对CUDA core了解的太肤浅了，这样的话，是不是我应该分配24个块，每个线程处理我的一个点信息呢！

system · 2013 年9 月 18 日 03:13

LZ您好：

您的理解依然是不对的，CUDA一般需要成千上万的线程协同工作才能发挥效能，这一点和CPU是不同的。一般需要将您算法中每一步的处理拿出来，予以充分并行化。

请您适当了解CUDA的基本原理，这将有助于您撰写和优化您的代码。

祝您好运~

system · 2013 年9 月 18 日 05:45

谢谢版主的宝贵意见，我会继续努力的

system · 2013 年9 月 18 日 05:53

我这个程序，数据量，并行度都不是很适合CUDA编写，最后在GPU上的运行速度比CPU上快不了多少，一倍都没有，真是费了很大的功夫，效果还真不怎么样。

system · 2013 年9 月 18 日 05:56

LZ您好：

确实有些算法和计算任务并不适合GPU计算，我们不能为了GPU而GPU。

但是您这个算法究竟如何，还是请您在熟悉CUDA之后，再行评估。

祝您好运~

程序优化 求指导

程序优化求指导