1

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

これが私がmexに使用したコードです（FFT2mxという名前の2D 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;

静的ライブラリを使用してmexコードをコンパイルします。

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

すべてがうまく機能します、それはただ遅いです。

FFTWライブラリは、次の引数を使用してコンパイルされています。

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

このコードは、2つのクアッドコアAMD Opteron（tm）を備えた1つのクラスターノードで実行しており、次のコマンドでテストしています。

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

魔女が戻る：

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

mexコードをどのように調整できますか？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の間）、セグメンテーション違反が発生しましたが、その理由がわかりません。ポインタで初期化する別の方法を試しました。また、対応する静的プランを操作するには、プランのポインターが静的である必要があることもどこかで読みました。私が間違っているのを見たことがありますか？

洞察に感謝します。

4

1 に答える 1

2

問題は、各FFTの計画を作成および破棄していることです。計画の作成は通常、FFT自体よりもはるかに時間がかかります。理想的には、プランを1回だけ作成して破棄し、それを同じディメンションの連続するFFTに何度も再利用します。

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

または、3つのMEX関数を使用することもできます。1つはプランの作成用、1つは特定のプランでFFTを実行するため、もう1つはプランを破棄するためです。

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

もう1つ、コマンドに追加--enable-sseして./configure、FFTWバタフライでSIMDコード生成を有効にすることができます。

于 2012-04-17T21:38:08.780 に答える