Why calculating time of FFT is shorter than element-to-element multiplication in intel MKL?

Question

I have vector of 1024*4608 elements (Original_signal) which is stored in one-dimention array.

And I enlarged the Original_signal to Expand_signal by copying every 1024 elements 32 times to 1024*32*4608.

Then I use a Com_array of 1024*32 to do the element-to-element multiplication with the Expand_signal and do the 1024FFT of the After multiplying array.

The core code is like follows:

//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 ");

The time of element-to-element multiplication is 700ms（After removing the memcpy cost）, And the time of FFT is 500ms.

The complex multiplication computation of FFT is N/2log2N And the element-to-element multiplication is N.

In this project N=1024. FFT is 5 times slower than element-to-element multiplication in theory. Why is faster in actual.

Any way to speed up the project?

(notice that Com_array is symmetrical)

Answer 1

In this project N=1024. FFT is 5 times slower than element-to-element multiplication in theory. Why is faster in actual?

As was indicated in comments, the time complexity of the FFT gives you a relative measure for various FFT lengths, up to some constant factor. This factor becomes important when trying to compare with other computations. Also, your analysis assumes that the performance is limited by floating point operations, where in reality the actual performance seems limited by other factors such as special case handling (e.g. NaN, Inf), memory and cache access.

Any way to speed up the project?

Since your performance bottleneck is around the complex element-wise vector multiplication operation, the following will focus on improving the performance of that operation.

I do not have MKL to perform actual benchmarks, but it is probably fair to assume that the vmcMul implementation is both fairly robust to special cases such as NaN and Inf, and fairly optimized in the circumstances.

If you do not need the robustness against the special cases, run on a SSE3 processor, can guarantee that your array sizes are multiples of 2, and that they are 16-bytes aligned, then you may get some performance gains by using a simplified implementation such as the following (based on 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://stackoverflow.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;
  }
}

Once the multiplication is optimized, your implementation could still be improved by getting rid of the extra copies to Temp_signal with memcpy by directly multiplying the Orignal_signal many times with different portions of Com_array, as shown below:

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;
  }
}

This last step would give you another ~20% performance improvement compared to your implementation with vmcMul replaced by packed_vec_mult.

Finally since the loop performs operation on independent blocks, you may be able to get significantly higher throughput (but similar latency) by launching parallel computations over multiple thread, so that the CPU is always kept busy instead of waiting for data in transit to/from memory. My tests suggests somewhere around a factor 2 improvement, but results may differ depending on your specific machine.