为啥 FFT 的计算时间比 intel MKL 中的元素到元素乘法要短？

Posted 2023-02-16

【中文标题】为啥 FFT 的计算时间比 intel MKL 中的元素到元素乘法要短？

【英文标题】：Why calculating time of FFT is shorter than element-to-element multiplication in intel MKL?为什么 FFT 的计算时间比 intel MKL 中的元素到元素乘法要短？

【发布时间】：2019-03-19 02:51:14

【问题描述】：

我有 1024*4608 个元素的向量（Original_signal），它存储在一维数组中。

我通过将每 1024 个元素复制 32 次到 1024*32*4608 来将 Original_signal 扩大到 Expand_signal。

然后我使用 1024*32 的 Com_array 与 Expand_signal 进行元素到元素的乘法，并在乘法后的数组中进行 1024FFT。

核心代码如下：

//initialize Original_signal
MKL_Complex8 *Original_signal = new MKL_Complex8[1024*4608];
for (int i=0; i<4608; i++)

  for (int j=0; j<1024; j++)
    
      Original_signal[j+i*1024].real=rand();
      Original_signal[j+i*1024].imag=rand();
    
 
//Com_array
MKL_Complex8 *Com_array= new MKL_Complex8[32*1024];
for (int i=0; i<32; i++)
  
    for (int j=0; j<1024; j++)
      
        Com_array[j+i*1024].real=cosf(2*pi*(i-16.0)/10.0*j^2);
        Com_array[j+i*1024].imag=sinf(2*pi*(i-16.0)/10.0*j^2);
      
  


//element-to-element multiplication
MKL_Complex8 *Temp_signal= new MKL_Complex8[1024*32];
MKL_Complex8 *Expand_signal= new MKL_Complex8[1024*32*4608];

gettimeofday(&Bgn_Time, 0);

for (int i=0; i<4608; i++)
  
    for (int j=0; j<32; j++)
      
        memcpy(Temp_signal+j*1024, Original_signal+i*1024, 1024*sizeof(MKL_Complex8));
      
      vmcMul(1024*32, Temp_signal, Com_array, Expand_signal+i*1024*32);
  

gettimeofday(&End_Time, 0);
double time_used = (double)(End_Time.tv_sec-Bgn_Time.tv_sec)*1000000+(double)(End_Time.tv_usec-Bgn_Time.tv_usec);
printf("element-to-element multiplication use time %fus\n, time_used ");


//FFT
DFTI_DESCRIPTOR_HANDLE h_FFT = 0;
DftiCreateDescriptor(&h_FFT, DFTI_SINGLE, DFTI_COMPLEX, 1, 1024);
DftiSetValue(h_FFT, DFTI_NUMBER_OF_TRANSFORMS, 32*4608);
DftiSetValue(h_FFT, DFTI_INPUT_DISTANCE, 1024);
DftiCommitDescriptor(h_FFT);


gettimeofday(&Bgn_Time, 0);

DftiComputeForward(h_FFT,Expand_signal);

gettimeofday(&End_Time, 0);
double time_used = (double)(End_Time.tv_sec-Bgn_Time.tv_sec)*1000000+(double)(End_Time.tv_usec-Bgn_Time.tv_usec);
printf("FFT use time %fus\n, time_used ");

元素到元素相乘的时间是700ms（去掉memcpy成本后），FFT的时间是500ms。

FFT的复数乘法计算为N/2log2N，元素到元素的乘法为N。

在这个项目中 N=1024。理论上，FFT 比元素到元素的乘法慢 5 倍。为什么实际上更快。

有什么方法可以加快项目速度？

（注意 Com_array 是对称的）

【问题讨论】：

您的时序可能包括大量的 I/O。对于元素乘法，您有 2N 次读取。对于 FFT，它是 N 个读数。在 FFT 情况下，函数调用开销也更少。您可能还需要检查 CPU/核心调度计划，以查看多个 FFT 是否并行完成，以及 vcMul 是否也是这种情况。

一般说一种算法的时间复杂度是N，另一种算法的时间复杂度是N log N，并不是说这些数字是可比的。在这两种情况下，都有一个常数因子（C1 * N vs C2 * N * log(N)）将在每种情况下都不同。插入排序是 O(N*N) 而快速排序是 O(N log N) - 但是对于短列表，插入排序通常更快，因为（隐含的）常数更小。

欢迎来到 ***！您是否对 mode of multiplication 使用特定值？其中一些模式执行错误检查或提高准确性。这些特征可能会导致计算速度变慢。长度为 N 的 dft 的 flop 计数约为 5Nlog_2(N) (//Cooley Tukey 算法)。 fftw.org/fftw-paper.pdf 对于 N=1024，每个值大约对应 50 次翻牌。对于实际信号，它可以除以 2。它确实比 1 乘法大得多！

@francis 我在 vmcMul 中专门使用了 VML_EP 模式，它可以加速乘法并降低精度。但时间成本仍然大于 FFT。

@SleuthEye 我将此线程绑定在 CPU 的一个核心中。所以我认为他们在同一个环境中。

【参考方案1】：

在这个项目中 N=1024。理论上，FFT 比元素到元素的乘法慢 5 倍。为什么实际更快？

正如 cmets 中所指出的，FFT 的时间复杂度为您提供了各种 FFT 长度的相对度量，直至某个常数因子。当试图与其他计算进行比较时，这个因素变得很重要。此外，您的分析假设性能受到浮点运算的限制，而实际上实际性能似乎受到其他因素的限制，例如特殊情况处理（例如NaN、Inf）、内存和缓存访问。

有什么方法可以加快项目速度？

由于您的性能瓶颈在于复杂的逐元素向量乘法运算，以下将重点关注提高该运算的性能。

我没有 MKL 来执行实际的基准测试，但可以公平地假设 vmcMul 实现对于 NaN 和 Inf 等特殊情况都相当稳健，并且在这种情况下得到了相当优化。

如果您不需要针对特殊情况的鲁棒性，在 SSE3 处理器上运行，可以保证您的数组大小是 2 的倍数，并且它们是 16 字节对齐的，那么您可以通过使用获得一些性能提升一个简化的实现，如下所示（基于Sebastien's answer to another post）：

#include <pmmintrin.h>
#include <xmmintrin.h>

// Computes and element-by-element multiplication of complex vectors "a" and "b" and
// stores the results in "c".
// Vectors "a", "b" and "c" must be:
//   - vectors of even length N
//   - 16-bytes aligned
// Special cases such as NaN and Inf are not handled.
//
// based on https://***.com/questions/3211346/complex-mul-and-div-using-sse-instructions#4884057
void packed_vec_mult(int N, MKL_Complex8* a, MKL_Complex8* b, MKL_Complex8* c)

  int M = N/2;

  __m128* aptr = reinterpret_cast<__m128*>(a);
  __m128* bptr = reinterpret_cast<__m128*>(b);
  __m128* cptr = reinterpret_cast<__m128*>(c);
  for (int i = 0; i < M; i++)
  
    __m128 t0 = _mm_moveldup_ps(*aptr);
    __m128 t1 = *bptr;
    __m128 t2 = _mm_mul_ps(t0, t1);
    __m128 t3 = _mm_shuffle_ps(t1, t1, 0xb1);
    __m128 t4 = _mm_movehdup_ps(*aptr);
    __m128 t5 = _mm_mul_ps(t4, t3);
    *cptr = _mm_addsub_ps(t2, t5);

    ++aptr;
    ++bptr;
    ++cptr;
  

一旦优化了乘法，您的实现仍然可以通过直接将Orignal_signal 与Com_array 的不同部分多次相乘来消除Temp_signal 和memcpy 的额外副本来改进您的实现，如下所示：

MKL_Complex8* outptr = Expand_signal;
for (int i=0; i<4608; i++)

  for (int j=0; j<32; j++)
  
    packed_vec_mult(1024, Original_signal+i*1024, Com_array+j*1024, outptr);
    outptr += 1024;
  

与将vmcMul 替换为packed_vec_mult 的实现相比，最后一步将使您的性能再提高约20%。

最后，由于循环在独立的块上执行操作，您可以通过在多个线程上启动并行计算来获得显着更高的吞吐量（但相似的延迟），以便 CPU 始终保持忙碌而不是等待传输中的数据到/从内存。我的测试表明大约有 2 倍的改进，但结果可能会因您的特定机器而异。