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

No.4 r_mulのチューニングメモリーのリニアアクセス化
　ソース

/* 実数fft後処理+乗算+実数ifft前処理 */ /* c = a * b ; a=b=c OK */ void *r_mul( int nx,ad_t c[],ad_t a[],ad_t b[] ) { double wr,wi,tmpr,tmpi,r; ad_t anr,amr,bnr,bmr; int m,n,mr,nr ; n = 1<< nx ; anr.r = a[0].r + a[0].i ; anr.i = a[0].r - a[0].i ; anr.r *= b[0].r + b[0].i ; anr.i *= b[0].r - b[0].i ; c[0].r = (anr.r + anr.i)*0.5 ; c[0].i = (anr.r - anr.i)*0.5 ; anr = a[1]; bnr = b[1]; c[1].r = anr.r * bnr.r - anr.i * bnr.i ; c[1].i = anr.r * bnr.i + anr.i * bnr.r ; r = M_PI /(double)n; for( nr = 2 ; nr < n ; nr+=2 ) { m = bitrv_v(nx,nr); mr = bitrv_v(nx,n-m); anr = a[nr]; amr = a[mr]; bnr = b[nr]; bmr = b[mr]; wi = sc[n+m].r*0.5; wr = sc[n+m].i*0.5+0.5; /* fft 後処理 a */ tmpr = wr * (anr.r - amr.r) - wi * (anr.i + amr.i); tmpi = wr * (anr.i + amr.i) + wi * (anr.r - amr.r); anr.r = anr.r - tmpr ; anr.i = anr.i - tmpi; amr.r = amr.r + tmpr ; amr.i = amr.i - tmpi; /* fft 後処理 b */ tmpr = wr * (bnr.r - bmr.r) - wi * (bnr.i + bmr.i); tmpi = wr * (bnr.i + bmr.i) + wi * (bnr.r - bmr.r); bnr.r = bnr.r - tmpr ; bnr.i = bnr.i - tmpi; bmr.r = bmr.r + tmpr ; bmr.i = bmr.i - tmpi; /* 乗算 */ tmpr = anr.r ; anr.r = tmpr * bnr.r - anr.i * bnr.i ; anr.i = tmpr * bnr.i + anr.i * bnr.r ; tmpr = amr.r ; amr.r = tmpr * bmr.r - amr.i * bmr.i ; amr.i = tmpr * bmr.i + amr.i * bmr.r ; /* ifft 前処理 */ tmpr = wr * (anr.r - amr.r) + wi * (anr.i + amr.i); tmpi = wr * (anr.i + amr.i) - wi * (anr.r - amr.r); c[nr].r = anr.r - tmpr ; c[nr].i = anr.i - tmpi; c[mr].r = amr.r + tmpr ; c[mr].i = amr.i - tmpi; } }

この処理でmrもほぼリニアになる。
後で解ったが、ifftとほぼ同じLOOPで処理可能。
ただ m は駄目だ。

結果

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

fftとifftの劇的改善に比べランダムアクセスが残るためか、余り改善しない。
sin,cosテーブルをリニアにアクセス出来るように並べ替えれば高速化可能だが、fftと共用不可となる。


No.5 r_mul でsin,cos関数を使用。
　ソース
　　r_mul用sin,cosテーブル作成部を削除。
　　r_mulをifftのLOOPの記述に合わせる。

ad_t *sc = NULL; void *mk_sincos_table(int nx) { int n ,m,nn; double r ; if( sc != NULL ){ free(sc); sc = NULL ; } if( nx <= 0 ){ return ; } n=1<<nx; sc=(ad_t *)malloc(sizeof(ad_t)*n); r = 2.0 * M_PI / (double)n; for ( m=0;m<n;m++){ sc[m].r = 0.0; sc[m].i = 0.0; } /* fft,ifft (nx) 0 ~ pai */ for (m=0 ;m<n/4;m++){ sc[n/2+m].r = cos(r*(double)m); sc[n/2+m].i = sin(r*(double)m); } sc[n/2+n/4].r = 0.0; sc[n/2+n/4].i = 1.0; for (m=1 ;m<n/4;m++){ sc[n-m].r = -sc[n/2+m].r; sc[n-m].i = sc[n/2+m].i; } for( nn=n/4 ; nn > 0 ; nn/=2 ){ for (m=0 ;m<nn;m++){ sc[nn+m] = sc[nn*2+m*2]; } } } /* 実数fft後処理+乗算+実数ifft前処理 */ /* c = a * b ; a=b=c OK */ void *r_mul( int nx,ad_t c[],ad_t a[],ad_t b[] ) { double wr,wi,tmpr,tmpi,r; ad_t anr,amr,bnr,bmr; int m,n,mr,nr,j,k ; n = 1<< nx ; anr.r = a[0].r + a[0].i ; anr.i = a[0].r - a[0].i ; anr.r *= b[0].r + b[0].i ; anr.i *= b[0].r - b[0].i ; c[0].r = (anr.r + anr.i)*0.5 ; c[0].i = (anr.r - anr.i)*0.5 ; anr = a[1]; bnr = b[1]; c[1].r = anr.r * bnr.r - anr.i * bnr.i ; c[1].i = anr.r * bnr.i + anr.i * bnr.r ; r = M_PI /(double)n; m = 0; for(j= 2 ; j < n; j *= 2) { for(k = 0; k < j; k+=2) { nr = j + k ; mr = j*2-1-k ; m = inc_rvbit(nx-1,m); #ifdef FFT_INDEX_CHECK if(nr < 0 ) { printf("r_mul index under nr\n"); exit ;}; if(mr< 0 ) { printf("r_mul index under mr\n"); exit ;}; if(m< 0 ) { printf("r_mul index under m\n"); exit ;}; if(nr >=n ) { printf("r_mul index over nr\n"); exit ;}; if(mr >=n ) { printf("r_mul index over mr\n"); exit ;}; if(m >=n ) { printf("r_mul index over m\n"); exit ;}; #endif anr = a[nr]; amr = a[mr]; bnr = b[nr]; bmr = b[mr]; wi = cos(r*(double)m)*0.5; wr = sin(r*(double)m)*0.5+0.5; /* fft 後処理 a */ tmpr = wr * (anr.r - amr.r) - wi * (anr.i + amr.i); tmpi = wr * (anr.i + amr.i) + wi * (anr.r - amr.r); anr.r = anr.r - tmpr ; anr.i = anr.i - tmpi; amr.r = amr.r + tmpr ; amr.i = amr.i - tmpi; /* fft 後処理 b */ tmpr = wr * (bnr.r - bmr.r) - wi * (bnr.i + bmr.i); tmpi = wr * (bnr.i + bmr.i) + wi * (bnr.r - bmr.r); bnr.r = bnr.r - tmpr ; bnr.i = bnr.i - tmpi; bmr.r = bmr.r + tmpr ; bmr.i = bmr.i - tmpi; /* 乗算 */ tmpr = anr.r ; anr.r = tmpr * bnr.r - anr.i * bnr.i ; anr.i = tmpr * bnr.i + anr.i * bnr.r ; tmpr = amr.r ; amr.r = tmpr * bmr.r - amr.i * bmr.i ; amr.i = tmpr * bmr.i + amr.i * bmr.r ; /* ifft 前処理 */ tmpr = wr * (anr.r - amr.r) + wi * (anr.i + amr.i); tmpi = wr * (anr.i + amr.i) - wi * (anr.r - amr.r); c[nr].r = anr.r - tmpr ; c[nr].i = anr.i - tmpi; c[mr].r = amr.r + tmpr ; c[mr].i = amr.i - tmpi; } } }

結果

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

一言『CPUキャッシュヒットミスによるペナルティは、三角関数より重い』
r_mul専用 sin,cosテーブルは、fftのデータサイズの1/2程度で作成が可能だが、データの配置がかなり面倒
未確認だが、nx 別に作成する必要がなさそう。この場合fftのデータサイズの1/4程度になる。
少し考えよう。


●その他の確認した項目

最大サイズのFFT結果

n=536870912=1<<29 radix=1000 fft 33.229422 , r_mul 15.924549 , ifft 39.121055 , error 0.437500 n=536870912=1<<29 radix=100 fft 33.393006 , r_mul 15.857790 , ifft 45.415891 , error 0.003906 n=1073741824=1<<30 radix=100 fft 69.092204 , r_mul 31.908605 , ifft 69.060423 , error 0.008301

n=1<<29　で radix=1000は、やはり無理か。
n=1<<30は、VMのメモリー52GBでLinuxが無応答になり原因の調査すら出来ない。メモリーリークは無いようだが何故だ？
VMのメモリー56GBでは正常。
※ Linuxを最新の状態に更新したらVMのメモリー52GBで問題無く動作。バグを踏んだらしい。
n=1<<30でのメモリー使用量は、fft領域で16GB*2 sin,cosテーブルで現状16GB。計48GB
sin,cosテーブルは、半分以下に切り詰めれるが、加減算やディスク保存に使う領域等で合計50GBは必要。
最大 n=1<<30で処理したいが、VM環境では無理そう。


多重実行性能
　fft_sp_04　は、No.5のモジュール

$ (./fft_sp_04 &) ;(./fft_sp_04 &) ;(./fft_sp_04 &) ;(./fft_sp_04 &) ; $ mk_sincos_table 1.296122 mk_sincos_table 1.310993 mk_sincos_table 1.338312 mk_sincos_table 1.342632 n=67108864=1<<26 radix=1000 fft 8.507075 , r_mul 2.069274 , ifft 8.423515 , error 0.046875 n=67108864=1<<26 radix=1000 fft 8.480474 , r_mul 2.064933 , ifft 8.439081 , error 0.046875 n=67108864=1<<26 radix=1000 fft 8.513428 , r_mul 2.074727 , ifft 8.491021 , error 0.046875 n=67108864=1<<26 radix=1000 fft 8.557486 , r_mul 2.067748 , ifft 8.503273 , error 0.046875

r_mulは、４多重で実行時間はほぼ同じ。これは並列実行数に比例して性能向上が期待できる。
処理時間の殆どがsin,cos関数によるものなので、sin,cosテーブルをリニアにアクセス出来るようにした場合、４倍以上の性能向上が期待出来る。
fft,ifftは、４多重で実行時間は約２倍。たぶんメモリー帯域を食潰して半分程度の性能が出ていないと思われる。
CPUキャッシュのヒット率を向上させれば良い。
つまりCPUキャッシュに乗る程度の小さな領域に小分けにし集中的に処理させるようにする。