FFTW与Matlab FFT

Question

我在matlab中心发布了这个,但没有得到任何回复所以我想我会转发到这里.

我最近在Matlab中编写了一个简单的例程,它在for循环中使用FFT; FFT主导了计算.我在mex中编写了相同的例程仅用于实验目的,它调用了FFTW 3.3库.事实证明,对于非常大的数组,matlab例程比mex例程运行得更快(大约快两倍).mex例程使用智慧并执行相同的FFT计算.我也知道matlab使用FFTW,但它们的版本是否可能稍微优化一下？我甚至使用了FFTW_EXHAUSTIVE标志,它对大型数组的速度仍然是MATLAB的两倍.此外,我确保我使用的matlab是单线程的"-singleCompThread"标志,我使用的mex文件不在调试模式.只是好奇,如果是这种情况 - 或者如果有一些优化,matlab正在使用我不知道的引擎盖.谢谢.

这是mex部分:

void class_cg_toeplitz::analysis() {
// This method computes CG iterations using FFTs
    // Check for wisdom
    if(fftw_import_wisdom_from_filename("cd.wis") == 0) {
        mexPrintf("wisdom not loaded.\n");
    } else {
        mexPrintf("wisdom loaded.\n");
    }

    // Set FFTW Plan - use interleaved FFTW
    fftw_plan plan_forward_d_buffer;    
    fftw_plan plan_forward_A_vec;       
    fftw_plan plan_backward_Ad_buffer;
    fftw_complex *A_vec_fft;
    fftw_complex *d_buffer_fft;
    A_vec_fft = fftw_alloc_complex(n);
    d_buffer_fft = fftw_alloc_complex(n);

    // CREATE MASTER PLAN - Do this on an empty vector as creating a plane 
    // with FFTW_MEASURE will erase the contents; 
    // Use d_buffer
    // This is somewhat dangerous because Ad_buffer is a vector; but it does not
    // get resized so &Ad_buffer[0] should work
    plan_forward_d_buffer = fftw_plan_dft_r2c_1d(d_buffer.size(),&d_buffer[0],d_buffer_fft,FFTW_EXHAUSTIVE);
    plan_forward_A_vec = fftw_plan_dft_r2c_1d(A_vec.height,A_vec.value,A_vec_fft,FFTW_WISDOM_ONLY);
    // A_vec_fft.*d_buffer_fft will overwrite d_buffer_fft
    plan_backward_Ad_buffer = fftw_plan_dft_c2r_1d(Ad_buffer.size(),d_buffer_fft,&Ad_buffer[0],FFTW_EXHAUSTIVE);

    // Get A_vec_fft
    fftw_execute(plan_forward_A_vec);

    // Find initial direction - this is the initial residual
    for (int i=0;i<n;i++) {
        d_buffer[i] = b.value[i];
        r_buffer[i] = b.value[i];
    }    

    // Start CG iterations
    norm_ro = norm(r_buffer);
    double fft_reduction = (double)Ad_buffer.size(); // Must divide by size of vector because inverse FFT does not do this
    while (norm(r_buffer)/norm_ro > relativeresidual_cutoff) {        
        // Find Ad - use fft
        fftw_execute(plan_forward_d_buffer);    
        // Get A_vec_fft.*fft(d) - A_vec_fft is only real, but d_buffer_fft
        // has complex elements; Overwrite d_buffer_fft        
        for (int i=0;i<n;i++) {
            d_buffer_fft[i][0] = d_buffer_fft[i][0]*A_vec_fft[i][0]/fft_reduction;
            d_buffer_fft[i][1] = d_buffer_fft[i][1]*A_vec_fft[i][0]/fft_reduction;
        }        
        fftw_execute(plan_backward_Ad_buffer); 

        // Calculate r'*r
        rtr_buffer = 0;
        for (int i=0;i<n;i++) {
            rtr_buffer = rtr_buffer + r_buffer[i]*r_buffer[i];
        }    

        // Calculate alpha
        alpha = 0;
        for (int i=0;i<n;i++) {
            alpha = alpha + d_buffer[i]*Ad_buffer[i];
        }    
        alpha = rtr_buffer/alpha;

        // Calculate new x
        for (int i=0;i<n;i++) {
            x[i] = x[i] + alpha*d_buffer[i];
        }   

        // Calculate new residual
        for (int i=0;i<n;i++) {
            r_buffer[i] = r_buffer[i] - alpha*Ad_buffer[i];
        }   

        // Calculate beta
        beta = 0;
        for (int i=0;i<n;i++) {
            beta = beta + r_buffer[i]*r_buffer[i];
        }  
        beta = beta/rtr_buffer;

        // Calculate new direction vector
        for (int i=0;i<n;i++) {
            d_buffer[i] = r_buffer[i] + beta*d_buffer[i];
        }  

        *total_counter = *total_counter+1;
        if(*total_counter >= iteration_cutoff) {
            // Set total_counter to -1, this indicates failure
            *total_counter = -1;
            break;
        }
    }

    // Store Wisdom
    fftw_export_wisdom_to_filename("cd.wis");

    // Free fft alloc'd memory and plans
    fftw_destroy_plan(plan_forward_d_buffer);
    fftw_destroy_plan(plan_forward_A_vec);
    fftw_destroy_plan(plan_backward_Ad_buffer);
    fftw_free(A_vec_fft);
    fftw_free(d_buffer_fft);
};

这是matlab部分:

% Take FFT of A_vec.
A_vec_fft = fft(A_vec); % Take fft once

% Find initial direction - this is the initial residual 
x = zeros(n,1); % search direction
r = zeros(n,1); % residual
d = zeros(n+(n-2),1); % search direction; pad to allow FFT
for i = 1:n
    d(i) = b(i); 
    r(i) = b(i); 
end

% Enter CG iterations
total_counter = 0;
rtr_buffer = 0;
alpha = 0;
beta = 0;
Ad_buffer = zeros(n+(n-2),1); % This holds the product of A*d - calculate this once per iteration and using FFT; only 1:n is used
norm_ro = norm(r);

while(norm(r)/norm_ro > 10^-6)
    % Find Ad - use fft
    Ad_buffer = ifft(A_vec_fft.*fft(d)); 

    % Calculate rtr_buffer
    rtr_buffer = r'*r;

    % Calculate alpha    
    alpha = rtr_buffer/(d(1:n)'*Ad_buffer(1:n));

    % Calculate new x
    x = x + alpha*d(1:n);

    % Calculate new residual
    r = r - alpha*Ad_buffer(1:n);

    % Calculate beta
    beta = r'*r/(rtr_buffer);

    % Calculate new direction vector
    d(1:n) = r + beta*d(1:n);      

    % Update counter
    total_counter = total_counter+1; 
end

在时间方面,对于N = 50000和b = 1:n,使用mex需要大约10.5秒,使用matlab需要4.4秒.我正在使用R2011b.谢谢

Answer 1

由于我不知道MATLAB FFT实现的任何细节,所以有一些观察而不是确定的答案:

• 根据您的代码,我可以看到两种速度差异的解释:
  • 速度差异由FFT优化级别的差异来解释
  • MATLAB中的while循环执行次数要少得多

我将假设你已经研究过第二个问题,并且迭代次数是可比的.(如果不是,这很可能是一些准确性问题,值得进一步调查.)

现在,关于FFT速度比较:

• 是的,理论上说FFTW比其他高级FFT实现更快,但只要你比较苹果和苹果,它就是相关的:在这里你要比较下一级的实现,在汇编级别,不仅是选择算法但其对特定处理器的实际优化以及具有不同技能的软件开发人员正在发挥作用
• 我在一年中在许多处理器上组装优化或审查优化的FFT(我在基准测试行业),优秀的算法只是故事的一部分.有些注意事项对您编码的体系结构非常具体(考虑延迟,指令调度,寄存器使用优化,内存中数据排列,分支采用/不采用延迟等)以及产生差异的因素与选择算法同样重要.
• 在N = 500000的情况下,我们还讨论了大型内存缓冲区:更多优化的另一扇门,它可以快速地为您运行代码的平台提供特定的功能:您如何设法避免缓存未命中将不会受到算法,以及数据流如何以及软件开发人员可能使用哪些优化来有效地将数据输入和输出内存.
• 虽然我不知道MATLAB FFT实现的细节,但我很确定一大批DSP工程师已经(并且仍然)在优化它,因为它是这么多设计的关键.这很可能意味着MATLAB拥有正确的开发人员组合,可以产生更快的FFT.

Answer 2

由于低级和体系结构特定的优化,这是经典的性能提升.

Matlab使用来自英特尔MKL(数学核心库)二进制文件(mkl.dll)的FFT .这些是英特尔为英特尔处理器优化(在汇编级别)的例程.即便在AMD上,它似乎也提供了不错的性能提升.

FFTW看起来像普通的c库,没有经过优化.因此使用MKL的性能提升.