FFTW 到顶部 Matlab FFT 的优化

Posted 2023-02-16

技术标签:

【中文标题】FFTW 到顶部 Matlab FFT 的优化【英文标题】：Optimization of FFTW to top Matlab FFT 【发布时间】：2012-04-17 21:28:43 【问题描述】：

Matlab 中的 FFT 不允许选择有多少线程正在执行计算 (http://***.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 次之间）后，我现在遇到了一些分段错误，我不知道为什么。我尝试了不同的方法来通过指针初始化。我还在某处读到计划的指针必须是静态的才能与相应的静态计划一起使用。你看到我做错了什么吗？

再次感谢您的见解。

【问题讨论】：

【参考方案1】：

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

如果您为相同大小的 FFT 重复调用 MEX，那么您可以memoize 计划（例如，保留静态计划变量和维度，并且仅根据需要重新创建计划，即当尺寸变化）。

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

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

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

【讨论】：

感谢您的评论。由于我通过一次调用函数来测试代码，因此计划生成不应该是惩罚，对吧？计划生成是的惩罚 - 创建（和销毁）计划需要时间。尝试在 MEX 中循环执行 1000 个 FFT，并将其与 1000 个 MATLAB FFT 进行比较。然后你会看到很大的不同。当应用多个FFT时，是在保存智慧，再次加载后创建一个新计划，这就是记忆的意思。我不确定如何将计划本身传回 matlab，然后传回 mex 函数。但是我仍然不明白为什么当我应用 single FFT 时，我的代码比 maltab 慢得多（它也必须计算其计划）。顺便说一句，FFT 重复 100 次（围绕 fftw_execute），我的 Matlab 得到 11 秒，我的 FFTW 代码得到 13 秒，matlab 得到 56 秒 - FFTW 得到 76 秒，重复 500 次......我仍然缺少一些东西 GPU 在 FFT 上的速度非常快如果您不计算从设备获取数据所需的时间。仔细观察，您似乎仍然在为每次调用创建/破坏一个计划，并且还在分配未释放的内存（即泄漏内存）。您可能希望简化代码并在尝试做太多其他事情（例如线程和智慧等）之前先弄清基本逻辑。

以上是关于FFTW 到顶部 Matlab FFT 的优化的主要内容，如果未能解决你的问题，请参考以下文章

fftw 来自幅度和相位阵列的 C++ 逆二维 FFT

fft 算法的基准测试方法

SSE图像算法优化系列十一：使用FFT变换实现图像卷积。

滚动视图不滚动到顶部

显示滚动到顶部如果顶部！= 0

2D R2C FFT 作为 1D FFT 使用 FFTW