FFTのチューニングと測定4

No.7 fft,ifft キャシュヒット率向上
　やったことは、小分けにして集中的に処理させるようにした。
　　MB_MAX : 集中的に処理するfftのデータサイズ　128KB　か　256KB　が最速となった。
　　MB_WIN : 上記を分割する最小サイズ　4KB が最速。（ページサイズと同じ）。
　その他改修
　　・sin,cos テーブルを1/2にした。1/4にしたかったが、頭が回らない。
　　　（既存 sin,cos テーブルの1/4しか使わないロジックは組めたが、いざテーブルを縮小しようとすると破綻する。
　　　　それほど難しくは無いはずなのだが、考えが纏まらない。まぁ後にしよう。）
　　・sin,cos テーブルを０～π/4 で計算。精度が若干上がる。
　　・sin,cos テーブル作成を１本化。

/* 処理単位 最小連続領域 fft内では、a[i]とa[m+i]の2領域必要 */ #define MB_WIN (4*1024) /* 処理領域サイズ */ #define MB_MAX (128*1024) ad_t *sc = NULL; ad_t *sc2 = NULL; double *scd = NULL; inline unsigned int bitrv_v(int nx, unsigned int a ); void mk_sincos_table(int nx) { int i,j,k,m,n,nr,nn; double r ; int mx ; if( sc != NULL ){ free(sc); free(sc2); free(scd); sc = NULL ; } if( nx <= 0 ){ return ; } n=1<<nx; mx = MB_MAX/(sizeof(ad_t)) ; /* fft,ifft (mx) 0 ~ pai */ sc=(ad_t *)malloc(sizeof(ad_t)*mx); r = 2.0 * M_PI / (double)mx; for (m=0 ;m<=mx/8;m++){ sc[mx/2+m].r = cos(r*(double)m); sc[mx/2+m].i = sin(r*(double)m); sc[mx/2+mx/4-m].r = sc[mx/2+m].i; sc[mx/2+mx/4-m].i = sc[mx/2+m].r; } for (m=1 ;m<mx/4;m++){ sc[mx-m].r = -sc[mx/2+m].r; sc[mx-m].i = sc[mx/2+m].i; } for( nn=mx/4 ; nn > 0 ; nn/=2 ){ for (m=0 ;m<nn;m++){ sc[nn+m] = sc[nn*2+m*2]; } } /* fft,ifft (n) 0 ~ pai/2 */ scd=(double *)malloc(sizeof(double)*n); r = 2.0 * M_PI / (double)n; scd[n/2+0] = 1.0; scd[n/2+1] = 0.0; for (m=1 ;m<=n/8;m++){ scd[(n/4+m)*2+0] = cos(r*(double)m); scd[(n/4+m)*2+1] = sin(r*(double)m); scd[(n/2-m)*2+0] = scd[(n/4+m)*2+1]; scd[(n/2-m)*2+1] = scd[(n/4+m)*2+0]; } for( nn=n/8 ; nn > 0 ; nn/=2 ){ for (m=0 ;m<nn;m++){ scd[(nn+m)*2+0] = scd[(nn+m)*4+0]; scd[(nn+m)*2+1] = scd[(nn+m)*4+1]; } } /* r_mul (n) 0 ~ pai/4 */ sc2=(ad_t *)malloc(sizeof(ad_t)*n/4); r = M_PI /(double)n; for(j= 2 ; j < n; j *= 2) { i = j / 4; for(k = 0; k < j/2; k+=2) { nr = j + k ; m = bitrv_v(nx,nr); /* wi = cos(r*(double)m)*0.5; */ /* wr = sin(r*(double)m)*0.5; */ sc2[i].i=cos(r*(double)m)*0.5; sc2[i].r=sin(r*(double)m)*0.5; i++; } } printf("sin cos table size %lu MB\n",(sizeof(ad_t)*mx+sizeof(double)*n+sizeof(ad_t)*n/4)/1024/1024); } /* ===== fft ======== */ static inline void xfft_1(int ns,int ne,ad_t a[],int ms,int me,int mw,int md) { int i,j,j1,j2,j3,k,m,cosm,sinm; double cosf,cosw,sinw,tmpr,tmpi; for( j1 = 0 ; j1 < ms/md ; j1 += mw){ for( j2 = ns ; j2 < ne ; j2 += ms*2 ){ for(m = ms; m > me ; m /= 2 ) { for( j3 = j1 ; j3 < m ; j3 += ms/md ){ if( j3 < m/2 ){ cosm = j3*2+m; sinm = j3*2+m+1; cosf = 1.0 ; }else{ cosm = j3*2 + 1; sinm = j3*2; cosf = -1.0 ; } for(j= j2 ; j < j2+ms*2; j += m*2) { for(k = 0 ; k < mw; k++) { i = j + k + j3 ; cosw = scd[k*2+cosm] * cosf; sinw =-scd[k*2+sinm]; tmpr = a[i].r - a[m+i].r; tmpi = a[i].i - a[m+i].i; a[i].r = a[i].r + a[m+i].r; a[i].i = a[i].i + a[m+i].i; a[m+i].r = tmpr * cosw - tmpi * sinw; a[m+i].i = tmpi * cosw + tmpr * sinw; #ifdef FFT_INDEX_CHECK if(i < 0 || m+i >=n) { printf("xfft_1 index error m+i=%d n=%d k=%d j=%d\n",m+i,n,k,j); exit(2) ;}; #endif } } } } } } } static inline void xfft_2(int ns,int ne,ad_t a[],int ms,int me,int mw,int md) { int i,j,j2,k,m; double cosw,sinw,tmpr,tmpi; for( j2 = ns ; j2 < ne ; j2 += ms*2 ){ for(m = ms; m > 0 ; m /= 2 ) { for(j= j2 ; j < j2+ms*2; j += m*2) { for(k = 0 ; k < m; k++) { i = j + k ; cosw = sc[m+k].r; sinw =-sc[m+k].i; tmpr = a[i].r - a[m+i].r; tmpi = a[i].i - a[m+i].i; a[i].r = a[i].r + a[m+i].r; a[i].i = a[i].i + a[m+i].i; a[m+i].r = tmpr * cosw - tmpi * sinw; a[m+i].i = tmpi * cosw + tmpr * sinw; #ifdef FFT_INDEX_CHECK if(i < 0 || m+i >=n) { printf("xfft_2 index error m+i=%d n=%d k=%d j=%d\n",m+i,n,k,j); exit(2) ;}; #endif } } } } } void fft(int nx ,ad_t a[]) { int n; int ms ,me ; int mx ,md ,mw ; n = 1<<nx; ms = n/2; mx = MB_MAX/(sizeof(ad_t)) ; mw = MB_WIN/(sizeof(ad_t)) ; md = mx / (mw*2) ; while ( ms*2 > mx ) { me = ms / md ; if( me*2 < mx ){ me = mx/2 ; } xfft_1( 0, n, a, ms, me, mw, md); ms=me; } xfft_2( 0, n, a, ms, me, mw, md); } /* ===== ifft ======== */ static inline void xifft_1(int ns,int ne,ad_t a[],int ms,int me,int mw,int md) { int i,j,j2,k,m; double cosw,sinw,tmpr,tmpi; for(j2= ns ; j2 < ne; j2 += me) { for(m = 1; m < me ; m *= 2 ) { for(j= j2 ; j < j2+me; j += m*2) { for(k = 0; k < m; k++) { i = j + k ; cosw = sc[m+k].r; sinw = sc[m+k].i; tmpr = cosw * a[m+i].r - sinw * a[m+i].i; tmpi = cosw * a[m+i].i + sinw * a[m+i].r; a[m+i].r = a[i].r - tmpr ; a[m+i].i = a[i].i - tmpi; a[i].r = a[i].r + tmpr ; a[i].i = a[i].i + tmpi; #ifdef FFT_INDEX_CHECK if(i < 0 || m+i >=n) { printf("xifft_1 index error m+i=%d n=%d k=%d j=%d\n",m+i,n,k,j); exit(2) ;}; #endif } } } } } static inline void xifft_2(int ns,int ne,ad_t a[],int ms,int me,int mw,int md) { int i,j,j1,j2,j3,k,m,cosm,sinm; double cosf,cosw,sinw,tmpr,tmpi; for( j1 = 0 ; j1 < me/(md*2) ; j1 += mw){ for( j2 = ns ; j2 < ne ; j2 += me ){ for(m = ms; m < me ; m *= 2 ) { for( j3 = j1 ; j3 < m ; j3 += me/(md*2) ){ if( j3 < m/2 ){ cosm = j3*2+m; sinm = j3*2+m+1; cosf = 1.0 ; }else{ cosm = j3*2 + 1; sinm = j3*2; cosf = -1.0 ; } for(j= j2 ; j < j2+me; j += m*2) { for(k = 0; k < mw; k++) { i = j + k + j3; cosw = scd[k*2+cosm] * cosf; sinw = scd[k*2+sinm]; tmpr = cosw * a[m+i].r - sinw * a[m+i].i; tmpi = cosw * a[m+i].i + sinw * a[m+i].r; a[m+i].r = a[i].r - tmpr ; a[m+i].i = a[i].i - tmpi; a[i].r = a[i].r + tmpr ; a[i].i = a[i].i + tmpi; #ifdef FFT_INDEX_CHECK if(i < 0 || m+i >=n) { printf("xifft_2 index error m+i=%d n=%d k=%d j=%d\n",m+i,n,k,j); exit(2) ;}; #endif } } } } } } } void ifft(int nx,ad_t a[]) { int n; int ms ,me ; int mx ,md ,mw ; long mel ; n = 1<< nx ; mx = MB_MAX/(sizeof(ad_t)) ; mw = MB_WIN/(sizeof(ad_t)) ; md = mx / (mw*2) ; me = mx; ms = 0; if ( n < me ){ me = n ; } xifft_1( 0, n, a, ms, me, mw, md); ms = me ; while ( ms < n ){ mel = (long)ms * md ; if( mel > n ){ me = n ; }else{ me = mel ; } xifft_2( 0, n, a, ms, me, mw, md); ms = me ; } }

結果

No.1 n=67108864=1<<26 radix=1000 fft 39.964405 , r_mul 5.791993 , ifft 38.881574 , error 0.085938 No.2 n=67108864=1<<26 radix=1000 fft 38.348005 , r_mul 5.744716 , ifft 41.086747 , error 0.085938 No.3 n=67108864=1<<26 radix=1000 fft 3.771606 , r_mul 4.460116 , ifft 4.872133 , error 0.046875 No.4 n=67108864=1<<26 radix=1000 fft 3.780859 , r_mul 3.618985 , ifft 4.872954 , error 0.046875 No.5 n=67108864=1<<26 radix=1000 fft 3.802325 , r_mul 1.972648 , ifft 4.882619 , error 0.046875 No.6 n=67108864=1<<26 radix=1000 fft 3.779546 , r_mul 0.482001 , ifft 4.864007 , error 0.046875 No.7 n=67108864=1<<26 radix=1000 fft 2.680603 , r_mul 0.491500 , ifft 4.124678 ,error 0.046875

fft,ifftは、MB_MAX を境に２分した、速度の向上は、MB_MAX 以下の処理(xfft_2,xifft_1)による。
MB_MAX 以上の処理では、単一スレッドではあまり改善をしていない。
しかし並列に動作させた場合の劣化が少ない。

$ (./fft_sp_10 &);(./fft_sp_10 &);(./fft_sp_10 &);(./fft_sp_10 &) $ sin cos table size 768 MB mk_sincos_table* 1.414342 sin cos table size 768 MB mk_sincos_table* 1.444200 sin cos table size 768 MB mk_sincos_table* 1.446548 sin cos table size 768 MB mk_sincos_table* 1.459807 n= 67108864=1<<26 radix= 1000 fft 3.55244350,r_mul 0.64616500,ifft 4.87817400,error 0.046875000 n= 67108864=1<<26 radix= 1000 fft 3.69948500,r_mul 0.82990200,ifft 4.91807500,error 0.046875000 n= 67108864=1<<26 radix= 1000 fft 3.77741450,r_mul 0.81672000,ifft 4.94080800,error 0.046875000 n= 67108864=1<<26 radix= 1000 fft 3.81712500,r_mul 0.81457200,ifft 5.02338200,error 0.046875000

あら？Linuxをアップデートしたせいかマルチスレッド性能が、10%程度劣化してる。（他でも確認）
仕事だったら大事だけど、ほぼ改修前の性能が出ているので良しとしよう。


n=1<<30　で処理したかったが、あきらめたほうが良さそう。
メモリーが厳しいのもあるが、n=1<<29とn=1<<30の計算出来る桁数が,1.5倍であるため
n=1<<30を超える桁数の乗算を行う場合、n=1<<29で行う方が、効率が良いかもしれない。

n=1<<30　の倍の桁を計算する場合。但し便宜上 ｎが倍異なるfftの速度比は倍、fftとifftの速度は同じ。

（１）n=1<<30　による場合
　　fftのデータ領域が、16GB ２つ。
　　・fft 5回とifft４回　計９回
　　・カラツバ方の場合 fft 6回とifft3回　計９回
　　　
（２）n=1<<29　による場合
　　fftのデータ領域が、8GB ５つ。
　　・fft 6回とifft9回　計１５回（n=1<<30相当で７．５回）
　　fftのデータ領域が、8GB ４つ。
　　・fft 8回とifft9回　計１６回（n=1<<30相当で８回）
このようになるとは限らないが、n=1<<29 で計算しても大差は無さそう。

基数(radix) 1xxx　と 5xxx の処理速度と精度
　全桁 radix-1 と radix-2 の乗算結果
　radix=5?? は、radix=1000 で 500 以上のとき 1000 引いて-500とし次の桁に１繰り越した値を使った場合の精度です。
　２進表現で１bit少なくなるので精度が上がる。
　radixの値は、radix=5?? で誤差 0.125 未満となるようにした。

sin cos table size 12288 MB mk_sincos_table* 19.288233 n= 32=1<< 5 radix=1000000 fft 0.00000144,r_mul 0.00000258,ifft 0.00000260,error 0.007812500 n= 64=1<< 6 radix=1000000 fft 0.00000175,r_mul 0.00000282,ifft 0.00000287,error 0.015625000 n= 128=1<< 7 radix=1000000 fft 0.00000241,r_mul 0.00000315,ifft 0.00000347,error 0.015625000 n= 256=1<< 8 radix=1000000 fft 0.00000399,r_mul 0.00000398,ifft 0.00000498,error 0.062500000 n= 512=1<< 9 radix=1000000 fft 0.00000743,r_mul 0.00000559,ifft 0.00000818,error 0.125000000 n= 1024=1<<10 radix=1000000 fft 0.00001481,r_mul 0.00000870,ifft 0.00001531,error 0.281250000 n= 2048=1<<11 radix= 100000 fft 0.00003047,r_mul 0.00001484,ifft 0.00003201,error 0.007812500 n= 4096=1<<12 radix= 100000 fft 0.00006592,r_mul 0.00002706,ifft 0.00006947,error 0.015625000 n= 8192=1<<13 radix= 100000 fft 0.00014421,r_mul 0.00005723,ifft 0.00013765,error 0.031250000 n= 16384=1<<14 radix= 100000 fft 0.00030691,r_mul 0.00011254,ifft 0.00030151,error 0.062500000 n= 32768=1<<15 radix= 100000 fft 0.00066610,r_mul 0.00022469,ifft 0.00064341,error 0.125000000 n= 65536=1<<16 radix= 100000 fft 0.00146280,r_mul 0.00040211,ifft 0.00136111,error 0.250000000 n= 131072=1<<17 radix= 10000 fft 0.00306610,r_mul 0.00079527,ifft 0.00330002,error 0.007812500 n= 262144=1<<18 radix= 10000 fft 0.00664170,r_mul 0.00161984,ifft 0.00769197,error 0.011718750 n= 524288=1<<19 radix= 10000 fft 0.01429983,r_mul 0.00341949,ifft 0.01750854,error 0.025390625 n= 1048576=1<<20 radix= 10000 fft 0.03075497,r_mul 0.00726500,ifft 0.03936812,error 0.062500000 n= 2097152=1<<21 radix= 10000 fft 0.06450894,r_mul 0.01460850,ifft 0.08666512,error 0.101562500 n= 4194304=1<<22 radix= 10000 fft 0.13735087,r_mul 0.02956750,ifft 0.18965150,error 0.250000000 n= 8388608=1<<23 radix= 1000 fft 0.28949825,r_mul 0.06009750,ifft 0.41275250,error 0.004882812 n= 16777216=1<<24 radix= 1000 fft 0.60363550,r_mul 0.12108900,ifft 0.88746600,error 0.011718750 n= 33554432=1<<25 radix= 1000 fft 1.25817400,r_mul 0.24190800,ifft 1.90104000,error 0.023437500 n= 67108864=1<<26 radix= 1000 fft 2.65220050,r_mul 0.48588900,ifft 4.07407800,error 0.046875000 n=134217728=1<<27 radix= 1000 fft 5.53910800,r_mul 0.97108300,ifft 8.65552200,error 0.109375000 n=268435456=1<<28 radix= 1000 fft 11.47869800,r_mul 1.96056300,ifft 18.41602400,error 0.187500000 n=536870912=1<<29 radix= 1000 fft 23.80241650,r_mul 3.93131700,ifft 38.82842100,error 0.406250000 n=1073741824=1<<30 radix= 100 fft 50.63409150,r_mul 7.84821800,ifft 82.09815700,error 0.007812500 sin cos table size 12288 MB mk_sincos_table* 19.289345 n= 32=1<< 5 radix= 500000 fft 0.00000146,r_mul 0.00000264,ifft 0.00000263,error 0.001953125 n= 64=1<< 6 radix= 500000 fft 0.00000176,r_mul 0.00000284,ifft 0.00000293,error 0.003906250 n= 128=1<< 7 radix= 500000 fft 0.00000245,r_mul 0.00000324,ifft 0.00000353,error 0.007812500 n= 256=1<< 8 radix= 500000 fft 0.00000398,r_mul 0.00000397,ifft 0.00000498,error 0.015625000 n= 512=1<< 9 radix= 500000 fft 0.00000736,r_mul 0.00000561,ifft 0.00000823,error 0.039062500 n= 1024=1<<10 radix= 500000 fft 0.00001491,r_mul 0.00000874,ifft 0.00001519,error 0.062500000 n= 2048=1<<11 radix= 50000 fft 0.00003078,r_mul 0.00001506,ifft 0.00003131,error 0.001953125 n= 4096=1<<12 radix= 50000 fft 0.00006631,r_mul 0.00002742,ifft 0.00006889,error 0.003906250 n= 8192=1<<13 radix= 50000 fft 0.00014308,r_mul 0.00005813,ifft 0.00014106,error 0.007812500 n= 16384=1<<14 radix= 50000 fft 0.00030473,r_mul 0.00010875,ifft 0.00030888,error 0.015625000 n= 32768=1<<15 radix= 50000 fft 0.00065135,r_mul 0.00023616,ifft 0.00065838,error 0.031250000 n= 65536=1<<16 radix= 50000 fft 0.00145519,r_mul 0.00039919,ifft 0.00136547,error 0.062500000 n= 131072=1<<17 radix= 5000 fft 0.00304454,r_mul 0.00079201,ifft 0.00327229,error 0.001586914 n= 262144=1<<18 radix= 5000 fft 0.00663593,r_mul 0.00160174,ifft 0.00765709,error 0.002929688 n= 524288=1<<19 radix= 5000 fft 0.01433247,r_mul 0.00351286,ifft 0.01745560,error 0.006347656 n= 1048576=1<<20 radix= 5000 fft 0.03078732,r_mul 0.00737153,ifft 0.03932924,error 0.011718750 n= 2097152=1<<21 radix= 5000 fft 0.06439238,r_mul 0.01480862,ifft 0.08659762,error 0.023437500 n= 4194304=1<<22 radix= 5000 fft 0.13750100,r_mul 0.02979375,ifft 0.18965050,error 0.062500000 n= 8388608=1<<23 radix= 500 fft 0.29046100,r_mul 0.06008050,ifft 0.41201500,error 0.001220703 n= 16777216=1<<24 radix= 500 fft 0.60163050,r_mul 0.12070600,ifft 0.88503500,error 0.002441406 n= 33554432=1<<25 radix= 500 fft 1.25629450,r_mul 0.24342700,ifft 1.90311300,error 0.005859375 n= 67108864=1<<26 radix= 500 fft 2.65134900,r_mul 0.48445000,ifft 4.07278700,error 0.011718750 n=134217728=1<<27 radix= 500 fft 5.53345550,r_mul 0.97872600,ifft 8.67555700,error 0.031250000 n=268435456=1<<28 radix= 500 fft 11.48626400,r_mul 1.95629100,ifft 18.36818100,error 0.050781250 n=536870912=1<<29 radix= 500 fft 23.78614350,r_mul 3.92355900,ifft 38.80194300,error 0.109375000 n=1073741824=1<<30 radix= 50 fft 50.65816300,r_mul 7.86358500,ifft 82.21162600,error 0.001953125

予定
　・fft,ifft内部マルチスレッド化
　　　fft,ifftをCPU分マルチスレッドで動かせれば良いが、数が不足する場合もある。
　　　単体でマルチスレッド化出来れば、扱い易い。
　・fft,ifftの動作検証
　　　1/3とかだと検算が楽だが今一。
　　　sqrt(2)をニュートン方で計算して検証は、GMPか？
　　　比較的組みやすいガウス＝ルジャンドルのアルゴリズムでπを計算してもよいか。
残件
　・整数化＆ディスクI/O
　　計算とディスクへの保存の都合上 int に10進8桁で保存予定。fftのradixとサイズが異なる。
　　ディスクI/Oは、mmapを使う予定。I/Oのスケジューリングの方が大変そう。
　・整数の乗算
　　fftを使うのは多分1000桁以上。それ以下は、カラツバ法などを使って乗算することになる。
　・n>1<<30 の乗算
　・多倍長割り算
　・他色々