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 倍。

这是我用于 mex 的代码（应用于 2D fft，名为 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

我在一个具有 2 个四核 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库的编译是否可以优化？有没有办法只使用 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 本身更耗时。理想情况下，您只需创建和销毁一次计划，然后将其多次重复用于相同维度的连续 FFT。

如果您为相同大小的 FFT 重复调用 MEX，那么您可能能够记住计划（例如，保留静态计划变量和维度，并且仅在需要时重新创建计划，即当维度发生变化时）。

或者，您可以拥有三个 MEX 功能 - 一个用于创建计划，一个用于使用给定计划运行 FFT，另一个用于销毁计划。

一旦你解决了上述架构问题，你应该考虑使用FFTW_MEASURE而不是FFTW_ESTIMATE获得更好的性能。

还有一件事：您可能想要添加--enable-sse到您的./configure命令中以在 FFTW 蝴蝶中启用 SIMD 代码生成。

于 2012-04-17T21:38:08.780 回答