本章节翻译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
两个或四个元素的移位给出一个相关得分。 任何其他移位都会得到零的相关得分。 这是如下计算的:
where N is the number of elements in the signal vectors and is the shift of sig2 relative to sig1.
实际信号包含更多数据(和噪声), 但是无论您是对齐 1D 信号,叠加 2D 图像 还是执行 3D 体积图像配准,原则都是相同的。 目标是找到最大化相关性的平移。但是,上面显示的暴力求和需每个移位进行: � 乘法和加法。在 1D、2D 和 3D 中,问题分别为 、 和 。
傅里叶相关算法是执行此计算的更有效方法, 因为它利用了傅里叶变换的 :
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。 现在必须将 的结果乘以 。 这可以使用手动编码的 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 提供了 一个合适的示例,以帮助你入门。