行列分解 2D FFT OpenCL 并行化

zhujiang73

搞了大约两周，其间大部分时间用在补习有关理论上了，上学时学的高数忘了不少了。

先把这个算法：

是这样的？  google  到这个：
zhujiang73 发表于 2010-12-29 10:42

改一下：

1.         for(int i = 0; i < mLen / 2; i++)
   
2.         {
   
3.                 float dAngle = -2 * FPI * i / mLen;
   
4.                 W[i].real = (float)cos(dAngle);
   
5.                 W[i].image = (float)sin(dAngle);
   
6.         }
   
7. 
   
8.         reverse(a,b,mLen,M);
   

复制代码

1. void reverse(int *a, int *b, int len, int M)
   
2. {
   
3.     int i,j;
   
4. 
   
5.     for(i=0; i<M; i++)
   
6.     {
   
7.         a[i] = 0;
   
8.     }
   
9. 
   
10.     b[0] = 0;
    
11.     for(i=1; i<len; i++)
    
12.     {
    
13.         j = 0;
    
14.         while(a[j] != 0)
    
15.         {
    
16.             a[j] = 0;
    
17.             j++;
    
18.         }
    
19. 
    
20.         a[j] = 1;
    
21.         b[i] = 0;
    
22.         for(j=0; j<M; j++)
    
23.         {
    
24.             b[i] = b[i]+a[j]*(int)pow(2,M-1-j);
    
25.         }
    
26.     }
    
27. }

复制代码

1. void fft(zcomplex *A, zcomplex *W, int fft_nLen, int fft_M)
   
2. {
   
3.         int p,dist;
   
4.         zcomplex C,B;
   
5. 
   
6.         for(int lev=1; lev<=fft_M; lev++)
   
7.         {
   
8.                 dist = 1<<(lev-1);
   
9.                 for(int t=0; t<dist; t++)
   
10.                 {
    
11.                         p = t*(1<<(fft_M-lev));
    
12. 
    
13.                         for(int i=t; i<fft_nLen; i+=1<<lev)
    
14.                         {
    
15.                                 B = Mul(A[i+dist],W[p]);
    
16.                                 C = Add(A[i],B);
    
17.                                 A[i+dist] = Sub(A[i],B);
    
18.                                 A[i].real = C.real;
    
19.                                 A[i].image = C.image;
    
20.                         }
    
21. 
    
22.                 }
    
23.         }
    
24. }
    

复制代码

然后再用 OpenCL 并行化，同时计算多行的 FFT，计算结果转置，再算一遍再转置回来就行了：

1. __kernel void fft4_nx( __global float4 *A,  __global const float2 *W,  const int nlen, const int fft_M)
   
2. {
   
3.         float4  f4_b;
   
4. 
   
5.         int i;
   
6.         int lev,dist,p,t;
   
7.         
   
8.         float2  w_l;
   
9.         float4  f4_tmp;
   
10.         float4  f4_tmp01;
    
11.         float4  f4_tmp02;
    
12.         
    
13.         float4   a4,b4;
    
14.         
    
15.         int  y = get_global_id(0);
    
16. 
    
17.         for(lev=1; lev<=fft_M; lev++)
    
18.         {
    
19.                 dist = 1<<(lev-1);
    
20.                 for(t=0; t<dist; t++)
    
21.                 {
    
22.                         p = t*(1<<(fft_M-lev));
    
23. 
    
24.                         for(i=t; i<nlen; i+=1<<lev)
    
25.                         {
    
26.                                 w_l = W[p];
    
27.                                 b4 = A[y*nlen+i+dist];
    
28.                                 f4_b.x = b4.x*w_l.x - b4.y*w_l.y;
    
29.                                 f4_b.y = b4.x*w_l.y + b4.y*w_l.x;
    
30. 
    
31.                                 f4_b.z = b4.z*W[p].x - b4.w*W[p].y;
    
32.                                 f4_b.w = b4.z*W[p].y + b4.w*W[p].x;
    
33.                                 
    
34.                                 a4 = A[y*nlen+i];
    
35.                                 
    
36.                                 A[y*nlen+i+dist] = a4 - f4_b;
    
37.                                 
    
38.                                 A[y*nlen+i] = a4 + f4_b;
    
39.                         }
    
40. 
    
41.                 }
    
42.         }
    
43. }
    

复制代码

zhujiang73

回复 1# zhujiang73


    可惜我的入门级 ATI4570 显卡不太给力，用了 GPU 的程序的计算速度和 T6670 CPU 上对应算法单线程程序速度差不多，估计如果是中高端显卡用了 GPU 的程序速度应该快很多。

erlangs

虽然不懂OpenCL,但这绝对是个有趣的东西

erlangs

回复  zhujiang73


    可惜我的入门级 ATI4570 显卡不太给力，用了 GPU 的程序的计算速度和 T6670 CP ...
zhujiang73 发表于 2011-01-13 13:01

   
怎么得知这块运算的瓶颈?

zhujiang73

怎么得知这块运算的瓶颈?
erlangs 发表于 2011-01-13 13:34

    FFT 的瓶颈就是那一大堆的复数运算，把多个 1D FFT 运算并行应该是可以加速的。

zhujiang73

虽然不懂OpenCL,但这绝对是个有趣的东西
erlangs 发表于 2011-01-13 13:27

    OpenCL 是显卡的第二职业，不打游戏时还可以算数学题。

ecjtubaowp

高深！

wsw1wsw2

任何显卡都能openCL?

zhujiang73

回复 1# zhujiang73


    这个 kernel 其实不太好，global 数据访问太多了，但是它比较简单，在我的 ATI4570 上能跑，AMD OpenCL SDK 里的 FFT 例子算法应该更好，但是我的显卡不支持，估计是需要 ATI5＊ 系列以上的显卡。{:3_201:}

zhujiang73

任何显卡都能openCL?
wsw1wsw2 发表于 2011-01-14 09:22

    AMD 系的显卡，推荐 ATI5* 以上的。

    NV  系的我没注意，不过 NV 搞通用计算比较早，支持得应该不错。