有没有办法在没有使用intel mkl换位的情况下计算另一个维度的2D FFT的1D FFT？

Question

我想使用mkl来计算存储为1D阵列的2D阵列的1D FFT。例如，

for (int j=0; j<NJ; j++) //rows
{
  for (int i=0; i<NI; i++) //columns
   {
     Pre_2D_array[i+j*NI].x=1.0;
     Pre_2D_array[i+j*NI].y=2.0;
   }
}

我想在行维度上计算Pre_2D_array的1D FFT。我能想到的唯一方法是重塑数组并像这样做FFT，

   for (int i=0; i<NI; i++) //columns
    {
      for (int j=0; j<NJ; j++) //rows
       {
         2D_array[j+i*NJ]=Pre_2D_array[i+j*NI];
       }
    }

DFTI_DESCRIPTOR_HANDLE desc_x = 0;
DftiCreateDescriptor(&desc_x, DFTI_PREC, DFTI_COMPLEX, 1, NJ);
DftiSetValue(desc_x, DFTI_NUMBER_OF_TRANSFORMS, NI);
DftiSetValue(desc_x, DFTI_INPUT_DISTANCE,  NJ);
DftiCommitDescriptor(desc_x);

DftiComputeForward(desc_x, 2D_array);

Althougth这可以得到正确的答案。但是当数组很大时，对原始数组进行转置（重新整形）会浪费太多时间。有没有办法在不重新整形阵列的情况下进行FFT？或者以尽可能快的速度重塑阵列的任何快速方法？

cpuinfo是：

processor   : 0
vendor_id   : GenuineIntel
cpu family  : 6
model       : 79
model name  : Intel(R) Xeon(R) CPU E5-2648L v4 @ 1.80GHz
stepping    : 1
microcode   : 0xb000022
cpu MHz     : 1795.882
cache size  : 35840 KB
physical id : 0
siblings    : 14
core id     : 0
cpu cores   : 14
apicid      : 0
initial apicid  : 0
fpu     : yes
fpu_exception   : yes
cpuid level : 20
wp      : yes
flags       : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc aperfmperf eagerfpu pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 fma cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch arat epb xsaveopt pln pts dtherm tpr_shadow vnmi flexpriority ept vpid fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm rdseed adx smap
bogomips    : 3591.76
clflush size    : 64
cache_alignment : 64
address sizes   : 46 bits physical, 48 bits virtual
power management:

Answer 1

FFTW库在istride等函数中引入了ostride或fftw_plan_many_dft()参数，以避免转置数组。该页面上的最后一个示例是第二维上的DFT。

同样，英特尔数学内核库引入了data layout configuration parameters，如DFTI_INPUT_STRIDES and DFTI_OUTPUT_STRIDES或DFTI_NUMBER_OF_TRANSFORMS。

第二维上的DFT可能看起来像（我没有测试过）：

DftiCreateDescriptor(&desc_x, DFTI_PREC, DFTI_COMPLEX, 1, NJ);
DftiSetValue(desc_x, DFTI_NUMBER_OF_TRANSFORMS, NI);
DftiSetValue(desc_x, DFTI_INPUT_STRIDES, &NI);
DftiSetValue(desc_x, DFTI_OUTPUT_STRIDES, &NI);
DftiSetValue(desc_x, DFTI_INPUT_DISTANCE,  1);
DftiSetValue(desc_x, DFTI_OUTPUT_DISTANCE,  1);
DftiCommitDescriptor(desc_x);

DFTI_OUTPUT_STRIDES被忽略为就地变换（DFTI_PLACEMENT=DFTI_INPLACE）。

Answer 2

据我所知，英特尔MKL不提供在数据之间跨越数据执行FFT的功能。

然而，FFTW确实如此。每4.4.1 Advanced Complex DFTs of the FFTW documentation：

fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany,
                             fftw_complex *in, const int *inembed,
                             int istride, int idist,
                             fftw_complex *out, const int *onembed,
                             int ostride, int odist,
                             int sign, unsigned flags);

该例程规划了多个多维复杂DFT，并扩展了fftw_plan_dft例程（参见复数DFT）以计算多个变换，每个变换具有排名等级和大小n。另外，变换数据不需要是连续的，但是它可以以任意步幅布置在存储器中。为了解释这些可能性，fftw_plan_many_dft添加了新参数howmany，{i,o}nembed，{i,o}stride和{i,o}dist。 FFTW基本接口（参见复杂DFT）提供专用于1,2和3级的例程，但高级接口仅处理一般级别的情况。

howmany是要计算的（非负）变换数。得到的计划计算howmany变换，其中第k个变换的输入位于位置in+k*idist（在C指针算术中），并且其输出位于位置out+k*odist。以这种方式获得的计划通常比为各个变换多次调用FFTW更快。基本的fftw_plan_dft接口对应于howmany=1（在这种情况下，dist参数被忽略）。

每个howmany变换都具有排名等级和大小n，如基本界面。此外，高级接口允许每个变换的输入和输出数组是较大秩秩数组的行主要子阵列，分别由inembed和onembed参数描述。 {i,o}nembed必须是长度为rank的数组，而n应该是元素小于或等于{i,o}nembed。将NULL传递给nembed参数等同于传递n（即，与基本界面中相同的物理和逻辑维度）。

步幅参数表明输入或输出数组的j-th元素分别位于j*istride或j*ostride。 （对于多维数组，j是普通的行主索引。）当与k-th循环中的howmany变换结合时，从上面开始，这意味着第（j，k）个元素位于j*stride+k*dist。 （基本的fftw_plan_dft界面对应于1的步幅。）

对于就地变换，输入和输出步幅和dist参数应该相同;否则，计划者可能会返回NULL。

该页面方便地提供了一个（有点令人困惑的）在2D阵列的列上执行1D FFT的示例：

使用10行和3列转换2d数组的每列：

   int rank = 1; /* not 2: we are computing 1d transforms */
   int n[] = {10}; /* 1d transforms of length 10 */
   int howmany = 3;
   int idist = odist = 1;
   int istride = ostride = 3; /* distance between two elements in 
                                 the same column */
   int *inembed = n, *onembed = n;

另请参阅How do I use fftw_plan_many_dft on a transposed array of data?以获取更多示例。

Answer 3

您不能沿数据的更高维度进行一维FFT。您需要先进行转置，以便FT维度是数据在RAM中连续的维度。

但是，它没有您想象的那么糟糕。在多核机器上，您可以轻松设置一些线程，其唯一的工作是预先/后期安排FT数据。