1. ホーム
  2. 質問と答え
  3. FFTWの最適化Matlab FFT

FFTWの最適化Matlab FFT

2012-04-17 14 views
1

MatlabのFFTでは、計算を行っているスレッドの数を選択できません（http://stackoverflow.com/questions/9528833/matlabs-fftn-gets-slower-with-multithreading ）。デフォルトでは、スタンドアローンのMATLAB上のすべてのコアが使用されます。しかし、クラスタでは、各ワーカーはデフォルトで単一のCPUで起動されます。より多くのコアで動作させることができます（maxNumCompThreads関数）。これは代数演算では完全に機能しますが、FFT関数は（奇妙なことに）シングルコアのままです。 私はfftwライブラリ（matlabのように）を使ってmexファイルを書いて、必要なコア数でfftを計算しました。しかし、FFTW_ESTIMATEプランナ（Matlabのデフォルトです）と明確な知恵を使ってコードを比較しようとすると、私のコードはMatlab fftよりも3〜4倍遅くなります。FFTWの最適化Matlab FFT

#include <stdlib.h> 
#include <stdio.h> 
#include <mex.h> 
#include <matrix.h> 
#include <math.h> 
#include </home/nicolas/Code/C/lib/include/fftw3.h>  
void FFTNDSplit(int NumDims, const int N[], double *XReal, double *XImag, double *YReal, double *YImag, int Sign) 
    { 
     fftw_plan Plan; 
     fftw_iodim Dim[NumDims]; 
     int k, NumEl; 
     for(k = 0, NumEl = 1; k < NumDims; k++) 
     { 
     Dim[NumDims - k - 1].n = N[k]; 
     Dim[NumDims - k - 1].is = Dim[NumDims - k - 1].os = (k == 0) ? 1 : (N[k-1] * Dim[NumDims-k].is); 
     NumEl *= N[k]; 
     } 

     //fftw_import_wisdom_from_filename("/home/nicolas/wisdom/wis"); 

     if(!(Plan = fftw_plan_guru_split_dft(NumDims, Dim, 0, NULL, XReal, 
              XImag, YReal, YImag, FFTW_ESTIMATE))) 
     mexErrMsgTxt("FFTW3 failed to create plan."); 

     if(Sign == -1) 
     fftw_execute_split_dft(Plan, XReal, XImag, YReal, YImag); 
     else 
     { 
     fftw_execute_split_dft(Plan, XImag, XReal, YImag, YReal); 
     } 

     //if(!fftw_export_wisdom_to_filename("/home/nicolas/wisdom/wis")) 
     // mexErrMsgTxt("FFTW3 failed to save wisdom."); 

     fftw_destroy_plan(Plan); 
     return; 
    } 


    void mexFunction(int nlhs, mxArray *plhs[], 
        int nrhs, const mxArray *prhs[]) 
    { 

     int i, j,numCPU; 
     int NumDims; 
     const mwSize *N; 

     if (nrhs != 2) { 
      mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
        "Two input argument required."); 
     } 

     if (!mxIsDouble(prhs[0])) { 
      mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
        "Array must be double"); 
     } 

     numCPU = (int) mxGetScalar(prhs[1]); 
     if (numCPU > 8) { 
      mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
        "NumOfThreads < 8 requested"); 
     } 


     /*if (!mxIsComplex(prhs[0])) { 
      mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
        "Array must be complex"); 
     }*/ 

     NumDims = mxGetNumberOfDimensions(prhs[0]); 
     N = mxGetDimensions(prhs[0]); 

     plhs[0] = mxCreateDoubleMatrix(0, 0, mxCOMPLEX); 
     mxSetDimensions(plhs[0], N, NumDims); 
     mxSetData(plhs[0], mxMalloc(sizeof(double) * mxGetNumberOfElements(prhs[0]))); 
     mxSetImagData(plhs[0], mxMalloc(sizeof(double) * mxGetNumberOfElements(prhs[0]))); 

     fftw_init_threads(); 
     fftw_plan_with_nthreads(numCPU); 

     FFTNDSplit(NumDims, N, (double *) mxGetPr(prhs[0]), (double *) mxGetPi(prhs[0]), 
       mxGetPr(plhs[0]), mxGetPi(plhs[0]), -1); 

    }

関連MATLABコード：

function fft2mx(X,NumCPU) 

FFT2mx(X,NumCPU)/sqrt(size(X,1)*size(X,2)); 
return;

Iは、静的を使用してMEXコードをコンパイルここ

私は（FFT2mxという名前の2次元FFTのために適用される、）MEXに使用されるコードであります図書館：

mex FFT2mx.cpp /home/nicolas/Code/C/lib/lib/libfftw3.a /home/nicolas/Code/C/lib/lib/libfftw3_threads.a

すべてがうまくいきますが、それはちょっと遅いです。

FFTWライブラリは、次の引数を使用してコンパイルされています：私は2クアッドコアAMD Opteronプロセッサ（TM）で1つのクラスタノード上でこのコードを実行していると私はでテスト

CC="gcc ${BUILD64} -fPIC" CXX="g++ ${BUILD64} -fPIC" \ 
./configure --prefix=/home/nicolas/Code/C/lib --enable-threads && 
make 
make install

：

A = randn([2048 2048])+ i*randn([2048 2048]); 
tic, fft2mx(A,8); toc; 
tic, fftn(A); toc;

魔女リターン：私のMEXコードは

Elapsed time is 0.482021 seconds. 
Elapsed time is 0.151630 seconds.

を調整することができますか？ fftwライブラリのコンパイルは最適化できますか？ ESTIMATEプランナのみを使用してfftwアルゴリズムを高速化する方法はありますか？

私は洞察力を求めています。ありがとうございました。

編集：私は今後、いくつかのセグメンテーションフォールトが発生したのです

# include <string.h> 
# include <stdlib.h> 
# include <stdio.h> 
# include <mex.h> 
# include <matrix.h> 
# include <math.h> 
# include </home/nicolas/Code/C/lib/include/fftw3.h> 

char *Wisfile = NULL; 
char *Wistemplate = "%s/.fftwis"; 
#define WISLEN 8 

void set_wisfile(void) 
{ 
    char *home; 
    if (Wisfile) return; 
    home = getenv("HOME"); 
    Wisfile = (char *)malloc(strlen(home) + WISLEN + 1); 
    sprintf(Wisfile, Wistemplate, home); 
} 

void cleanup(void) { 
    static fftw_plan PlanForward; 
    static int planlen; 
    static double *pr, *pi, *pr2, *pi2; 
    mexPrintf("MEX-file is terminating, destroying array\n"); 
    fftw_destroy_plan(PlanForward); 
    fftw_free(pr2); 
    fftw_free(pi2); 
    fftw_free(pr); 
    fftw_free(pi); 
} 


void mexFunction(int nlhs, mxArray *plhs[], 
       int nrhs, const mxArray *prhs[]) 
{ 

    int i, j, numCPU, NumDims; 
    const mwSize *N; 
    fftw_complex *out, *in1; 
    static double *pr, *pi, *pr2, *pi2; 
    static int planlen = 0; 
    static fftw_plan PlanForward; 
    fftw_iodim Dim[NumDims]; 
    int k, NumEl; 
    FILE *wisdom; 

    if (nrhs != 2) { 
     mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
       "Two input argument required."); 
    } 

    if (!mxIsDouble(prhs[0])) { 
     mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
       "Array must be double"); 
    } 

    numCPU = (int) mxGetScalar(prhs[1]); 
    if (numCPU > 8) { 
     mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
       "NumOfThreads < 8 requested"); 
    } 


    if (!mxIsComplex(prhs[0])) { 
     mexErrMsgIdAndTxt("MATLAB:FFT2mx:invalidNumInputs", 
       "Array must be complex"); 
    } 


    NumDims = mxGetNumberOfDimensions(prhs[0]); 
    N = mxGetDimensions(prhs[0]); 
    for(k = 0, NumEl = 1; k < NumDims; k++) 
    { 
    Dim[NumDims - k - 1].n = N[k]; 
    Dim[NumDims - k - 1].is = Dim[NumDims - k - 1].os = (k == 0) ? 1 : (N[k-1] * Dim[NumDims-k].is); 
    NumEl *= N[k]; 
    } 

/* If different size, free/destroy */ 
    if(N[0] != planlen && planlen > 0) { 
    fftw_free(pr2); 
    fftw_free(pi2); 
    fftw_free(pr); 
    fftw_free(pi); 
    fftw_destroy_plan(PlanForward); 
    planlen = 0; 
    } 
    mexAtExit(cleanup); 


/* Init */ 

fftw_init_threads(); 
// APPROACH 1 
    //pr = (double *) mxGetPr(prhs[0]); 
    //pi = (double *) mxGetPi(prhs[0]); 

// APPROACH 2 
    pr = (double *) fftw_malloc(sizeof(double) * mxGetNumberOfElements(prhs[0])); 
    pi = (double *) fftw_malloc(sizeof(double) * mxGetNumberOfElements(prhs[0])); 
    tmp1 = (double *) mxGetPr(prhs[0]); 
    tmp2 = (double *) mxGetPi(prhs[0]); 
    for(k=0;k<mxGetNumberOfElements(prhs[0]);k++) 
    { 
    pr[k] = tmp1[k]; 
    pi[k] = tmp2[k]; 
    } 

    plhs[0] = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxCOMPLEX); 
    mxSetDimensions(plhs[0], N, NumDims); 
    mxSetData(plhs[0], (double*) fftw_malloc(sizeof(double) * mxGetNumberOfElements(prhs[0]))); 
    mxSetImagData(plhs[0], (double*) fftw_malloc(sizeof(double) * mxGetNumberOfElements(prhs[0]))); 

    pr2 = mxGetPr(plhs[0]); 
    pi2 = mxGetPi(plhs[0]); 

    fftw_init_threads(); 
    fftw_plan_with_nthreads(numCPU); 

/* Get any accumulated wisdom. */ 

    set_wisfile(); 
    wisdom = fopen(Wisfile, "r"); 
    if (wisdom) { 
    fftw_import_wisdom_from_file(wisdom); 
    fclose(wisdom); 
    } 

/* Compute plan */ 

//printf("%d",planlen); 
    if(planlen == 0) { 

fftw_plan_with_nthreads(numCPU); 
    PlanForward = fftw_plan_guru_split_dft(NumDims, Dim, 0, NULL, pr, pi, pr2, pi2, FFTW_MEASURE); 
    planlen = N[0]; 
    } 

/* Save the wisdom. */ 

    wisdom = fopen(Wisfile, "w"); 
    if (wisdom) { 
    fftw_export_wisdom_to_file(wisdom); 
    fclose(wisdom); 
    } 

/* execute */ 

    fftw_execute_split_dft(PlanForward, pr, pi, pr2, pi2); 
    fftw_cleanup_threads(); 
}

：私はあなたが（知恵と静的な計画を使用して）提案し、この更新されたコードを書いたものを考慮に入れる

いくつかの呼び出し（2から6までの間）を関数に渡すと、なぜ私は理解できません。ポインタで初期化する方法を変えました。私はまた、計画のポインタが対応する静的プランで動作するように静的でなければならないこともどこかに読んでいます。あなたが見ているものは何でも間違っていますか？あなたの洞察力のために再び

感謝。

出典

2012-04-17 Nicolas

答えて

2

問題は、各FFTのための計画を作成し、破壊しているということです。プランの作成は、通常、FFT自体よりもはるかに時間がかかります。理想的にはあなただけ作成し、一度計画を破棄し、同じ寸法（S）の連続したFFTのためにそれを何度も再利用します。

同じサイズのFFTに対してMEXを繰り返し呼び出す場合、は計画をメモすることができます（静的なプラン変数とディメンションを保持し、必要に応じてプランを再作成する、つまりディメンションが変更された場合など） 。

また、プランを作成するためのもの、特定のプランを使用してFFTを実行するもの、プランを破壊するものの3つのMEX機能を持つことができます。

上記のアーキテクチャ上の問題を修正したら、パフォーマンスを向上させるためにFFTW_ESTIMATEの代わりにFFTW_MEASUREを使用することを検討する必要があります。

もう1つ：FFTWの蝶でSIMDコード生成を有効にするには、./configureコマンドに--enable-sseコマンドを追加するとよいでしょう。

出典

2012-04-17 21:38:08

+0

コメントありがとうございます。私は関数を一回呼び出すだけでコードをテストしているので、プランの生成はペナルティではありません。 – Nicolas

+0

プランの生成*はペナルティです - プランを作成（および破壊）するには時間がかかります。 MEX内のループで1000個のFFTを実行し、1000個のMATLAB FFTと比較してみてください。そうすれば、大きな違いが見えます。 –

+0

複数のFFTを適用する場合、知恵を保存していて、再度読み込んだ後に新しい計画を作成することは、暗記という意味です。私はどのように計画自体をmatlabに戻してからmex関数に戻すかわからない。しかし、なぜ私が_単一FFTを適用すると、私のコードがmaltabより遅いのか（なぜなら、それも計画を計算しなければならない）、私はまだ理解していません。ところで、FFT（arround fftw_execute）の100回の繰り返しで、Matlabでは11秒、FFTWでは13秒、matlabでは56秒、500回繰り返されるFFTWでは76秒...私はまだ何かが欠けています – Nicolas

関連する問題

 関連する問題