FFTW到Matlab FFT的优化

Question

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慢3到4倍。

下面是我为mex使用的代码(用于2DFFT，名为FFT2mx)：

#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

我在一个带有两个四核AMD Opteron(tm)的集群节点上运行这段代码，我用：

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库的编译？在仅使用估计规划器的情况下，是否有方法加快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之间)之后，我现在遇到了一些分段错误，我不知道原因。我尝试了不同的方式通过指针初始化。我还在某个地方读到，计划的指针必须是静态的，才能使用相应的静态计划。你看到我做错什么了吗？

再次感谢你的见解。

Answer 1

问题是，您正在为每个FFT创建和破坏一个计划。创建一个计划通常比FFT本身要耗时得多。理想情况下，您只创建和销毁一个计划一次，然后对相同维度的连续FFT多次重复使用它。

如果您正在为相同大小的FFT重复调用MEX，那么您可能能够回溯计划(例如，保持静态计划变量和维度，并且只在需要时重新创建计划，即当维度发生变化时)。

或者，您可以有三个MEX函数-一个用于创建计划，一个用于在给定的计划中运行FFT，另一个用于销毁计划。

一旦解决了上述架构问题，就应该考虑使用FFTW_MEASURE而不是FFTW_ESTIMATE来提高性能。

还有一件事:您可能希望将--enable-sse添加到./configure命令中，以便在FFTW蝴蝶中启用SIMD代码生成。