oneAPI GPU 优化指南 - 使用 oneMKL 高效实现傅里叶相关

本文介绍如何使用 oneMKL 高效实现傅里叶相关。交叉相关有诸多应用，传统暴力求和计算量大，傅里叶相关算法更高效，可利用 oneMKL 组合实现。文中给出使用显式缓冲区和统一共享内存（USM）两种实现方式，并对比了性能，还指出避免不必要数据传输可提升性能。

本章节翻译by chenchensmail@163.com 原文：Efficiently Implementing Fourier Correlation Using oneAPI Math Kernel Library (oneMKL) (intel.com)

使用 oneMKL 高效实现傅里叶相关#

现在，我们已经介绍了如何直接使用 oneMKL 内核函数，让我们来看看更复杂的数学运算：交叉相关。交叉相关有许多应用，例如：测量两个 1D 信号的相似性，查找最佳平移以叠加相似图像，体积医学图像分割等。

考虑以下简单信号，表示为 1 和 0 的向量:

Signal 1:  0  0  0  0  0  1  1  0
Signal 2:  0  0  1  1  0  0  0  0

这些信号被视为彼此的循环移位版本，因此将第二个信号相对于第一个信号移动三个元素将给出最大相关得分 2:

Signal 1:  0  0  0  0  0  1  1  0
Signal 2:           0  0  1  1  0  0  0  0

Correlation: (1 * 1) + (1 * 1) = 2

两个或四个元素的移位给出一个相关得分。任何其他移位都会得到零的相关得分。这是如下计算的：

$corr_{\alpha } = \sum_{i=0}^{N-1}sig1_{i} \times sig2_{i+\alpha }$

where N is the number of elements in the signal vectors and $\alpha$ is the shift of sig2 relative to sig1.

实际信号包含更多数据（和噪声），但是无论您是对齐 1D 信号，叠加 2D 图像还是执行 3D 体积图像配准，原则都是相同的。目标是找到最大化相关性的平移。但是，上面显示的暴力求和需每个移位进行： � 乘法和加法。在 1D、2D 和 3D 中，问题分别为 $O(N^{2})$ 、 $O(N^{3})$ 和 $O(N^{4})$ 。

傅里叶相关算法是执行此计算的更有效方法，因为它利用了傅里叶变换的 $O(NlogN)$ :

corr = IDFT(DFT(sig1) * CONJG(DFT(sig2)))

其中 DFT 是离散傅里叶变换， IDFT 是逆变换， CONJG 是复共轭。傅里叶相关算法可以使用 oneMKL 组合，其中包含优化的正向和反向变换以及复共轭乘法函数。因此，整个计算可以在加速器设备上执行。

在许多应用中，只有最终的相关结果很重要，因此所有内容必须从设备传回主机。

在以下示例中，将在设备上创建两个人工信号，就地转换，然后进行相关计算。主机将检索最终结果并报告最佳平移和相关得分。传统智慧认为缓冲区会提供最佳性能，因为它提供了对主机和设备之间数据移动的显式控制。

为了测试这个假设，让我们生成两个输入信号：

  // Create buffers for signal data. This will only be used on the device.
  sycl::buffer<float> sig1_buf{N + 2};
  sycl::buffer<float> sig2_buf{N + 2};

  // Declare container to hold the correlation result (computed on the device,
  // used on the host)
  std::vector<float> corr(N + 2);

通常会向信号添加随机噪声，以防止神经网络训练期间过度拟合，向图像添加视觉效果，或改善从次优探测器获得的信号的可检测性等。使用 oneMKL 中的简单随机数生成器初始化缓冲区以添加随机噪声：

   // Open new scope to trigger update of correlation result
   {
     sycl::buffer<float> corr_buf(corr);
 
     // Initialize the input signals with artificial data
     std::uint32_t seed = (unsigned)time(NULL); // Get RNG seed value
 7    oneapi::mkl::rng::mcg31m1 engine(Q, seed); // Initialize RNG engine
 8                                               // Set RNG distribution
 9    oneapi::mkl::rng::uniform<float, oneapi::mkl::rng::uniform_method::standard>
10        rng_distribution(-0.00005, 0.00005);
11
12    oneapi::mkl::rng::generate(rng_distribution, engine, N, sig1_buf); // Noise
13    oneapi::mkl::rng::generate(rng_distribution, engine, N, sig2_buf);

请注意，打开了一个新的作用域，并声明了一个缓冲区 corr_buf，用于相关结果。当此缓冲区超出范围时， corr 将在主机上更新。

人工信号放置在每个缓冲区的相对端，类似于上面的小示例：

 1    Q.submit([&](sycl::handler &h) {
 2      sycl::accessor sig1_acc{sig1_buf, h, sycl::write_only};
 3      sycl::accessor sig2_acc{sig2_buf, h, sycl::write_only};
 4      h.single_task<>([=]() {
 5        sig1_acc[N - N / 4 - 1] = 1.0;
 6        sig1_acc[N - N / 4] = 1.0;
 7        sig1_acc[N - N / 4 + 1] = 1.0; // Signal
 8        sig2_acc[N / 4 - 1] = 1.0;
 9        sig2_acc[N / 4] = 1.0;
10        sig2_acc[N / 4 + 1] = 1.0;
11      });
12    }); // End signal initialization

现在信号已经准备好了，让我们使用 oneMKL 中的 DFT 函数进行变换：

1    // Initialize FFT descriptor
2    oneapi::mkl::dft::descriptor<oneapi::mkl::dft::precision::SINGLE,
3                                 oneapi::mkl::dft::domain::REAL>
4        transform_plan(N);
5    transform_plan.commit(Q);
6
7    // Perform forward transforms on real arrays
8    oneapi::mkl::dft::compute_forward(transform_plan, sig1_buf);
9    oneapi::mkl::dft::compute_forward(transform_plan, sig2_buf);

单精度实际到复杂正向变换提交给 SYCL 队列，然后在两个数据上执行就地DFT。现在必须将 $DF T (sig1)$ 的结果乘以 $CONJG(DF T (sig2))$ 。这可以使用手动编码的 kernel 完成:

Q.submit([&](sycl::handler &h)
{
   sycl::accessor sig1_acc{sig1_buf, h, sycl::read_only};
   sycl::accessor sig2_acc{sig2_buf, h, sycl::read_only};
   sycl::accessor corr_acc{corr_buf, h, sycl::write_only};

   h.parallel_for<>(sycl::range<1>{N/2}, [=](auto idx)
   {
      corr_acc[idx*2+0] = sig1_acc[idx*2+0] * sig2_acc[idx*2+0] +
                          sig1_acc[idx*2+1] * sig2_acc[idx*2+1];
      corr_acc[idx*2+1] = sig1_acc[idx*2+1] * sig2_acc[idx*2+0] -
                          sig1_acc[idx*2+0] * sig2_acc[idx*2+1];
   });
});

但是，这种基本实现不太可能提供最佳的跨架构性能。幸运的是， oneMKL 函数 oneapi::mkl::vm::mulbyconj 可以用于此步骤。 mulbyconj 函数需要 ``std::complex<float>``输入，但是缓冲区被初始化为``float`` 数据类型。即使在正向变换后它们包含复杂数据，缓冲区也必须重新转换：

    auto sig1_buf_cplx =
        sig1_buf.template reinterpret<std::complex<float>, 1>((N + 2) / 2);
    auto sig2_buf_cplx =
        sig2_buf.template reinterpret<std::complex<float>, 1>((N + 2) / 2);
    auto corr_buf_cplx =
        corr_buf.template reinterpret<std::complex<float>, 1>((N + 2) / 2);
    oneapi::mkl::vm::mulbyconj(Q, N / 2, sig1_buf_cplx, sig2_buf_cplx,
                               corr_buf_cplx);

IDFT 步骤完成计算:

    // Perform backward transform on complex correlation array
    oneapi::mkl::dft::compute_backward(transform_plan, corr_buf);

当在此示例开始时打开的作用域关闭时，保存相关结果的缓冲区超出范围，导致强制更新主机容器。这是主机和设备之间唯一的数据传输。

使用显式缓冲区的完整傅里叶相关实现如下：

  #include <CL/sycl.hpp>
  #include <iostream>
  #include <mkl.h>
  #include <oneapi/mkl/dfti.hpp>
  #include <oneapi/mkl/rng.hpp>
  #include <oneapi/mkl/vm.hpp>
  
  int main(int argc, char **argv) {
    unsigned int N = (argc == 1) ? 32 : std::stoi(argv[1]);
    if ((N % 2) != 0)
      N++;
    if (N < 32)
      N = 32;
  
    // Initialize SYCL queue
    sycl::queue Q(sycl::default_selector_v);
    std::cout << "Running on: "
              << Q.get_device().get_info<sycl::info::device::name>() << std::endl;
  
    // Create buffers for signal data. This will only be used on the device.
    sycl::buffer<float> sig1_buf{N + 2};
    sycl::buffer<float> sig2_buf{N + 2};
 
    // Declare container to hold the correlation result (computed on the device,
    // used on the host)
    std::vector<float> corr(N + 2);
  
    // Open new scope to trigger update of correlation result
    {
      sycl::buffer<float> corr_buf(corr);
  
      // Initialize the input signals with artificial data
      std::uint32_t seed = (unsigned)time(NULL); // Get RNG seed value
      oneapi::mkl::rng::mcg31m1 engine(Q, seed); // Initialize RNG engine
                                                 // Set RNG distribution
      oneapi::mkl::rng::uniform<float, oneapi::mkl::rng::uniform_method::standard>
          rng_distribution(-0.00005, 0.00005);
  
      oneapi::mkl::rng::generate(rng_distribution, engine, N, sig1_buf); // Noise
      oneapi::mkl::rng::generate(rng_distribution, engine, N, sig2_buf);
  
      Q.submit([&](sycl::handler &h) {
        sycl::accessor sig1_acc{sig1_buf, h, sycl::write_only};
        sycl::accessor sig2_acc{sig2_buf, h, sycl::write_only};
        h.single_task<>([=]() {
          sig1_acc[N - N / 4 - 1] = 1.0;
          sig1_acc[N - N / 4] = 1.0;
          sig1_acc[N - N / 4 + 1] = 1.0; // Signal
          sig2_acc[N / 4 - 1] = 1.0;
          sig2_acc[N / 4] = 1.0;
          sig2_acc[N / 4 + 1] = 1.0;
        });
      }); // End signal initialization
  
      clock_t start_time = clock(); // Start timer
  
      // Initialize FFT descriptor
      oneapi::mkl::dft::descriptor<oneapi::mkl::dft::precision::SINGLE,
                                   oneapi::mkl::dft::domain::REAL>
          transform_plan(N);
      transform_plan.commit(Q);
  
      // Perform forward transforms on real arrays
      oneapi::mkl::dft::compute_forward(transform_plan, sig1_buf);
      oneapi::mkl::dft::compute_forward(transform_plan, sig2_buf);
  
      // Compute: DFT(sig1) * CONJG(DFT(sig2))
      auto sig1_buf_cplx =
          sig1_buf.template reinterpret<std::complex<float>, 1>((N + 2) / 2);
      auto sig2_buf_cplx =
          sig2_buf.template reinterpret<std::complex<float>, 1>((N + 2) / 2);
      auto corr_buf_cplx =
          corr_buf.template reinterpret<std::complex<float>, 1>((N + 2) / 2);
      oneapi::mkl::vm::mulbyconj(Q, N / 2, sig1_buf_cplx, sig2_buf_cplx,
                                 corr_buf_cplx);
  
      // Perform backward transform on complex correlation array
      oneapi::mkl::dft::compute_backward(transform_plan, corr_buf);
  
      clock_t end_time = clock(); // Stop timer
      std::cout << "The 1D correlation (N = " << N << ") took "
                << float(end_time - start_time) / CLOCKS_PER_SEC << " seconds."
                << std::endl;
  
    } // Buffer holding correlation result is now out of scope, forcing update of
      // host container
  
    // Find the shift that gives maximum correlation value
    float max_corr = 0.0;
    int shift = 0;
    for (unsigned int idx = 0; idx < N; idx++) {
      if (corr[idx] > max_corr) {
        max_corr = corr[idx];
        shift = idx;
      }
    }
    int _N = static_cast<int>(N);
    shift =
        (shift > _N / 2) ? shift - _N : shift; // Treat the signals as circularly
                                               // shifted versions of each other.
    std::cout << "Shift the second signal " << shift
              << " elements relative to the first signal to get a maximum, "
                 "normalized correlation score of "
              << max_corr / N << "." << std::endl;
  }

现在将使用统一共享内存(USM) 重新实现傅里叶相关算法，以与显式缓冲区进行比较。只突出两种实现的差异。首先，在 USM 中分配信号和相关数组，然后使用人工数据进行初始化：

   // Initialize signal and correlation arrays
   auto sig1 = sycl::malloc_shared<float>(N + 2, sycl_device, sycl_context);
   auto sig2 = sycl::malloc_shared<float>(N + 2, sycl_device, sycl_context);
   auto corr = sycl::malloc_shared<float>(N + 2, sycl_device, sycl_context);
 
   // Initialize input signals with artificial data
   std::uint32_t seed = (unsigned)time(NULL); // Get RNG seed value
   oneapi::mkl::rng::mcg31m1 engine(Q, seed); // Initialize RNG engine
                                              // Set RNG distribution
   oneapi::mkl::rng::uniform<float, oneapi::mkl::rng::uniform_method::standard>
       rng_distribution(-0.00005, 0.00005);
 
   // Warning: These statements run on the device.
   auto evt1 =
       oneapi::mkl::rng::generate(rng_distribution, engine, N, sig1); // Noise
   auto evt2 = oneapi::mkl::rng::generate(rng_distribution, engine, N, sig2);
   evt1.wait();
   evt2.wait();
 
   // Warning: These statements run on the host, so sig1 and sig2 will have to be
   // updated on the device.
   sig1[N - N / 4 - 1] = 1.0;
   sig1[N - N / 4] = 1.0;
   sig1[N - N / 4 + 1] = 1.0; // Signal
   sig2[N / 4 - 1] = 1.0;
   sig2[N / 4] = 1.0;
   sig2[N / 4 + 1] = 1.0;

其余的实现基本上是相同的，除了指向 USM 的指针被传递给 oneMKL 函数而不是 SYCL 缓冲区：

   // Perform forward transforms on real arrays
   evt1 = oneapi::mkl::dft::compute_forward(transform_plan, sig1);
   evt2 = oneapi::mkl::dft::compute_forward(transform_plan, sig2);
 
   // Compute: DFT(sig1) * CONJG(DFT(sig2))
   oneapi::mkl::vm::mulbyconj(
       Q, N / 2, reinterpret_cast<std::complex<float> *>(sig1),
       reinterpret_cast<std::complex<float> *>(sig2),
       reinterpret_cast<std::complex<float> *>(corr), {evt1, evt2})
       .wait();
 
   // Perform backward transform on complex correlation array
   oneapi::mkl::dft::compute_backward(transform_plan, corr).wait();

还需要释放分配的内存:

  sycl::free(sig1, sycl_context);
  sycl::free(sig2, sycl_context);
  sycl::free(corr, sycl_context);

USM 实现具有更熟悉的语法。它也更简单，因为它依赖于 SYCL runtime 处理的隐式数据传输。但是，程序员错误会影响性能。

请注意前面代码片段中的警告消息。 oneMKL 随机数生成引擎在设备上初始化，因此 sig1 和 sig2 在设备上使用随机噪声进行初始化。不幸的是，添加人工信号的代码在主机上运行，因此 SYCL runtime 必须使主机和设备数据一致。傅里叶相关中使用的信号通常很大，特别是在 3D 成像应用中，因此不必要的数据传输会降低性能。

直接在设备上更新信号数据可以使数据保持一致，从而避免不必要的数据传输：

  Q.single_task<>([=]() {
     sig1[N - N / 4 - 1] = 1.0;
     sig1[N - N / 4] = 1.0;
     sig1[N - N / 4 + 1] = 1.0; // Signal
     sig2[N / 4 - 1] = 1.0;
     sig2[N / 4] = 1.0;
     sig2[N / 4 + 1] = 1.0;
   }).wait();

显式缓冲区和 USM 实现现在具有等效的性能，表明 SYCL runtime 擅长避免不必要的数据传输（前提是程序员注意数据一致性）。

USM 中的完整傅里叶相关实现如下：

 #include <CL/sycl.hpp>
 #include <iostream>
 #include <mkl.h>
 #include <oneapi/mkl/dfti.hpp>
 #include <oneapi/mkl/rng.hpp>
 #include <oneapi/mkl/vm.hpp>
 
 int main(int argc, char **argv) {
   unsigned int N = (argc == 1) ? 32 : std::stoi(argv[1]);
   if ((N % 2) != 0)
     N++;
   if (N < 32)
     N = 32;

   // Initialize SYCL queue
   sycl::queue Q(sycl::default_selector_v);
   auto sycl_device = Q.get_device();
   auto sycl_context = Q.get_context();
   std::cout << "Running on: "
             << Q.get_device().get_info<sycl::info::device::name>() << std::endl;

   // Initialize signal and correlation arrays
   auto sig1 = sycl::malloc_shared<float>(N + 2, sycl_device, sycl_context);
   auto sig2 = sycl::malloc_shared<float>(N + 2, sycl_device, sycl_context);
   auto corr = sycl::malloc_shared<float>(N + 2, sycl_device, sycl_context);

   // Initialize input signals with artificial data
   std::uint32_t seed = (unsigned)time(NULL); // Get RNG seed value
   oneapi::mkl::rng::mcg31m1 engine(Q, seed); // Initialize RNG engine
                                             // Set RNG distribution
   oneapi::mkl::rng::uniform<float, oneapi::mkl::rng::uniform_method::standard>
       rng_distribution(-0.00005, 0.00005);
 
   auto evt1 =
       oneapi::mkl::rng::generate(rng_distribution, engine, N, sig1); // Noise
   auto evt2 = oneapi::mkl::rng::generate(rng_distribution, engine, N, sig2);
   evt1.wait();
   evt2.wait();

   Q.single_task<>([=]() {
      sig1[N - N / 4 - 1] = 1.0;
      sig1[N - N / 4] = 1.0;
      sig1[N - N / 4 + 1] = 1.0; // Signal
      sig2[N / 4 - 1] = 1.0;
      sig2[N / 4] = 1.0;
      sig2[N / 4 + 1] = 1.0;
    }).wait();

   clock_t start_time = clock(); // Start timer
 
   // Initialize FFT descriptor
   oneapi::mkl::dft::descriptor<oneapi::mkl::dft::precision::SINGLE,
                                oneapi::mkl::dft::domain::REAL>
       transform_plan(N);
   transform_plan.commit(Q);
 
   // Perform forward transforms on real arrays
   evt1 = oneapi::mkl::dft::compute_forward(transform_plan, sig1);
   evt2 = oneapi::mkl::dft::compute_forward(transform_plan, sig2);

   // Compute: DFT(sig1) * CONJG(DFT(sig2))
   oneapi::mkl::vm::mulbyconj(
       Q, N / 2, reinterpret_cast<std::complex<float> *>(sig1),
       reinterpret_cast<std::complex<float> *>(sig2),
       reinterpret_cast<std::complex<float> *>(corr), {evt1, evt2})
       .wait();
 
   // Perform backward transform on complex correlation array
   oneapi::mkl::dft::compute_backward(transform_plan, corr).wait();
 
   clock_t end_time = clock(); // Stop timer
   std::cout << "The 1D correlation (N = " << N << ") took "
             << float(end_time - start_time) / CLOCKS_PER_SEC << " seconds."
             << std::endl;
 
   // Find the shift that gives maximum correlation value
   float max_corr = 0.0;
   int shift = 0;
   for (unsigned int idx = 0; idx < N; idx++) {
     if (corr[idx] > max_corr) {
       max_corr = corr[idx];
       shift = idx;
     }
   }
   int _N = static_cast<int>(N);
   shift =
       (shift > _N / 2) ? shift - _N : shift; // Treat the signals as circularly
                                              // shifted versions of each other.
   std::cout << "Shift the second signal " << shift
             << " elements relative to the first signal to get a maximum, "
                "normalized correlation score of "
             << max_corr / N << "." << std::endl;

   // Cleanup
   sycl::free(sig1, sycl_context);
   sycl::free(sig2, sycl_context);
   sycl::free(corr, sycl_context);
 }

请注意，查找最大相关值位置的最后一步在主机上执行。最好在设备上进行此计算，特别是当输入数据很大时。幸运的是， maxloc reduction 是一种常见的并行模式，可以使用 SYCL 实现。这留给读者作为练习，但是 Data Parallel C++ 的图 14-11 提供了一个合适的示例，以帮助你入门。