intel onemkl api 黑客松---FFT算法的加速与优化

ninggan99

已于 2023-08-31 23:26:27 修改

阅读量261

点赞数 6

文章标签：算法 oneapi c语言

于 2023-08-31 19:48:36 首次发布

本文链接：https://blog.csdn.net/Reading8/article/details/132606755

版权

1 题目描述

使用oneMKL工具，对FFT算法进行加速与优化。

1.点击链接下载和使用最新版本oneMKL

2.调用oneMKL相应API函数，产生2048*2048个随机单精度实数()；

3.根据2产生的随机数据作为输入，实现两维RealtocomplexFFT参考代码；

4.根据2产生的随机数据作为输入，调用oneMKLAPI计算两维RealtocomplexFFT；

5.结果正确性验证，对3和4计算的两维FFT输出数据进行全数据比对（允许适当精度误差），输出“结果正确”或“结果不正确”信息；

6.平均性能数据比对（比如运行1000次），输出FFT参考代码平均运行时间和oneMKLFFT平均运行时间。

2 题目解析

2.1 FFT数学理论

2.1.1傅里叶变换

傅里叶变换（Fourier Transform，简称FT）是一种线性积分变换，用于函数（应用上称作“信号”）在时域和频域之间的变换。因其基本思想首先由法国学者约瑟夫·傅里叶系统地提出，所以以其名字来命名以示纪念。（连续）傅里叶变换的定义如下：

$\widehat f(s)= \int _ {-\infty }^ {\infty } e^ {-i2\pi } \xi xf(x)dx$ （2.1）

其中， $f$ 为原函数， $\widehat f$ 是原函数的傅里叶变换， $x$ 表示时间，其定义域为时域， $\xi$ 表示频率，其定义域为频域， $i$ 是虚数单位。通常情况下， $f$ 是一个实函数，而 $\widehat f$ 则是一个复数值函数。

2.1.2离散傅里叶变换

离散傅里叶变换（Discrete Fourier Transform，简称DFT），是傅里叶变换在时域和频域上都离散的形式，将信号的时域采样变换为频域采样。对于N点序列 $\{{x[n]}\}, 0\leq n< N$ ，它的傅里叶变换（DFT）是一个同样长度的序列：

$\hat{x}[k]=\sum_{n=0}^{N-1}e^{-i2\pi\frac{kn}{N}}x[n]$ (2.2)

其中 $e$ 为自然底数，𝑖是虚数单位。通常以符号 $\mathcal{F}$ 表示这一变换，即：

$\hat{x}=\mathcal{F}x$ (2.3)

快速傅里叶变换（Fast Fourier Transform,FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法。FFT能够将计算DFT的复杂度从 $O(n^2)$ 降低到 $O(nlon(n))$ ，其中𝑛为数据大小。在介绍快速傅里叶变换原理之前，我们先引入N次单位根。我们设 $\omega_{N}^{k}=e^{i2\pi\frac{k}{N}}=cos(2\pi\frac{k}{N})+i\cdot sin(2\pi\frac{k}{N})$ ，其几何意义是复平面上以原点为起点、单位圆的𝑁等分点为终点的向量。特别的， $\omega_{N}^{1}$ 被称作𝑁次单位根。

$\omega_{N}^{k}$ 有以下性质：

$\begin{aligned} &\omega_{N}^{k}=(\omega_{N}^{1})^{k_{\leftarrow}} \\ &\omega_{N}^{k}=\omega_{N}^{k-1}\cdot\omega_{N}^{1}\leftarrow \\ &\omega_{2N}^{2k}=\omega_{N}^{k} \\ &\omega_{N}^{k+\frac{N}{2}}=-\omega_{N|}^{k|_{L}} \\ &\omega_{N}^{k}=\overline{{{\omega_{N}^{-k}}}}_{\leftarrow} \end{aligned}$ (2.4)

若𝑁为偶数，使用 $\omega_{N}^{k}$ 的性质，我们可以进行如下推导

$\begin{aligned} \hat{x}[k]& =\sum_{n=0}^{N-1}e^{-i2\pi\frac{kn}{N}x[n]^{\leftarrow}} \\ &=\sum_{n=0}^{N-1}\omega_{N}^{-nk}x[n]^{\leftarrow} \\ &=\sum_{n=0}^{\frac{N}{2}-1}\omega_N^{-2nk}x[2n]+\sum_{n=0}^{\frac{N}{2}-1}\omega_N^{-(2n+1)k}x[2n+1]^{\leftarrow} \\ &=\sum_{n=0}^{\frac N2-1}\omega_{N}^{-2nk}x[2n]+\omega_{N}^{-k}\sum_{n=0}^{\frac N2-1}\omega_{N}^{-2nk}x[2n+1]^{\leftarrow} \\ &=\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-nk}x[2n]+\omega_{N}^{-k}\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-nk}x[2n+1]^{\leftarrow} \end{aligned}$

对于 $x[k+{\frac{N}{2}}]$ ，有

$\begin{gathered} \hat{x}[k+\frac{N}{2}] =\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-n(k+\frac{N}{2})}x[2n]+\omega_{N}^{-(k+\frac{N}{2})}\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-n(k+\frac{N}{2})}x[2n+1]_{\leftrightarrow} \\ \begin{aligned}&=\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-nk-n\frac{N}{2}}x[2n]+\omega_{N}^{-k-\frac{N}{2}}\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-nk-n\frac{N}{2}}x[2n+1\\&=\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-nk}x[2n]-\omega_{N}^{-k}\sum_{n=0}^{\frac{N}{2}-1}\omega_{\frac{N}{2}}^{-nk}x[2n+1]^{\leftarrow}\end{aligned} \end{gathered}$

因此只需要求得 $\sum_{\mathrm{n=0}}^{\frac{N}{2}-1}\omega_{N/2}^{-nk}x[2\mathrm{n}]$ 以及 $\sum_{\mathrm{n=0}}^{\frac{N}{2}-1}\omega_{N/2}^{-nk}x[2\mathrm{n}+1]$ ，即可求得 ${\hat{x}}[k]$ 以及 ${\hat{x}}[k+2]$ 。而求 $\sum_{\mathrm{n=0}}^{\frac{N}{2}-1}\omega_{N/2}^{-nk}x[2\mathrm{n}]|$ 以及 $\sum_{\mathrm{n=0}}^{\frac{\mathrm{N}}{2}-1}\omega_{N/2}^{-nk}x[2\mathrm{n}+1]$ 的过程相当于对𝑥的奇数下标部分和偶数下标部分分别做离散傅里叶变换。通过一次将原本的计算分解为奇数部分和偶数部分的操作，我们将原本离散傅里叶变换中计算所有 ${\hat{x}}[k]$ 所需的 $N^2$ 次乘法分解为2个向量长度为 $\frac{N}{2}$ 的需要 $\frac{N^2}{4}$ 次乘法的更小的离散傅里叶变换，继续地分解下去直到向量的长度为1，可以将总的计算复杂度缩减到 $O(nlon(n))$ 。

2.1.4 二维离散傅里叶变换

一般做影像处理的影像大多不是连续信号，而对于平面上的不连续信号，我们则需使用二维离散傅立叶变换。

假设输入的影像为𝑥，水平方向长度是𝑁，垂直方向长度是𝑀，二维离散傅立叶变换定义为：

$\widehat{x}[u,v]=\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}e^{-i2\pi(\frac{um}{M}+\frac{vn}{N})}x[m,n]$ （2.5）

二维离散傅里叶变换的基本实现方法是先对每一行进行一次一维离散傅里叶变换之后，再对每一列进行一次一维离散傅里叶变换。

2.2 实现细节

1、下载并按照相关引导安装相应版本的Intel®oneAPIBaseToolkitDownload the Intel® oneAPI Base ToolkitSelect your operating system, distribution channel and then download your customized installation of Intel® oneAPI.https://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html2、使用Onemkl数学函数库，需配置Intel®oneAPIBaseToolkit命令行开发环境，在环境中使用，具体操作为打开一个命令行窗口，输入指令

source /opt/intel/oneapi/setvars.sh intel64

setvars.sh前的路径根据实际安装情况调整,例如，setvars.sh安装在路径/home/ganning/intel/oneapi/下则实际指令为

source /home/ganning/intel/oneapi/setvars.sh intel64

(此环境配置方法仅在当前打开的命令行窗口可用，若需在新的命令行窗口使用，请再次配置)，需要注意的是，此步非常重要，是使用onemklapi的首要步骤。

3、这一步为可选选项，可以选择下载Intel®oneAPIDPC++/C++Compiler用intel的编译器来编译源码，在此次比赛中，我们使用了intel提供的的icpx编译器来编译（注意，icpx的使用需激活2、中的环境）。

4、生成单精度随机数时，需要使用OneMKL随机数生成器，使用icpx编译器编译链接OneMKL随机数生成器，需要在源文件中包含"mkl_vsl.h"头文件以及在链接时加入正确的链接指令，否则无法使用。例如，源文件为data_generate.cpp，OneMKL相关库文件安装位置在/home/ganning/intel/oneapi/mkl/latest/lib下，随机数生成器需要使用mkl_rt库，我们在编译时需使用指令：

icpx data_generate.cpp -o data_generate -L/home/ganning/intel/oneapi/mkl/latest/lib -lmkl_rt

进行编译链接，成功后使用./data_generate指令执行生成的可执行文件，观察运行结果，可以发现oneMKL随机数生成器已被正确调用。

data_generate.cpp：

#include <stdio.h>
#include "mkl_vsl.h"

int main() {
	//生成随机数个数
	int n = 256 * 256;
	//存储随机数的数组 
	float data[n]; 
    
	VSLStreamStatePtr stream;
	// 初始化随机数生成器流 
	vslNewStream(&stream, VSL_BRNG_MT19937, 777);
	//生成均匀分布的随机数 
	vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD, stream, n, data, 0.0, 1.0);
	// 删除随机数生成器流 
	vslDeleteStream(&stream);

	// 打印生成的随机数（每行 10 个） 
	for (int i = 0; i < n; i++) {
	printf("%f\t", data[i]);
	if ((i + 1) % 10 == 0) {
	printf("\n");
	}
	}
	printf("\n");
    
	return 0;
}

5、使用fftw3库作为验证onemkl中fft类api是否正确调用以及自写代码是否基本正确的基准，首先，我们安装fftw3库，然后根据FFTW编程的方法改写一个二维单精度实数到复数的随机傅立叶映射ground_truth.cpp

从上图可以看出，FFTW3库的使用大致是先用fftw_malloc分配输入输出内存，然后输入数据赋值，然后创建变换方案（fftw_plan），然后执行变换（fftw_execute），销毁方案（fftw_destory_plan），最后释放资源（fftw_free)，值得注意的是，输入输出数据的精度要求使得具体实际操作函数的使用有所不同，但是方法原理都一致。例如我们使用单精度输入，在分配输入输出内存时需要使用fftwf_malloc函数，创建变换方案需要使用fftwf_plan函数，销毁方案，释放资源函数也需换成对应的fftwf_destory_plan()和fftwf_free()函数，另外在使用icpx编译器时也需链接到fftw3f库。

ground_truth.cpp：

#include "fftw3.h"
#include <stdio.h>
#define COLS 4
#define ROWS 4

int main() {
    float (*in)[2];
    float (*out)[2];
    fftwf_plan p;

    int N = ROWS * COLS;

    in = (float(*)[2]) fftwf_malloc(sizeof(float) * 2 * N);
    out = (float(*)[2]) fftwf_malloc(sizeof(float) * 2 * N);

    float r2c_data[ROWS][COLS];
    printf("Before Transform:\n");
    for (int i = 0; i < ROWS; i++) {
        for (int j = 0; j < COLS; j++) {
            in[i * ROWS + j][0] = i + j; 
            in[i * ROWS + j][1] = 0.0;  
            r2c_data[i][j] = in[i * ROWS + j][0];
            printf("%7.2f\t", r2c_data[i][j]);
        }
        printf("\n");
    }


    p = fftwf_plan_dft_2d(ROWS, COLS, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftwf_execute(p);

    printf("After Transform:\n");
    for (int i = 0; i < N; i++) {
        printf("%7.2f + %7.2fi\t", out[i][0], out[i][1]);
        if ((i + 1) % COLS == 0) printf("\n");
    }
    printf("\n");

    fftwf_destroy_plan(p);
    fftwf_free(in);
    fftwf_free(out);

    return 0;
}

同样的，使用fftw3函数库在编译时需加入与使用Onemkl随机数生成器类似的编译链接指令

icpx ground_truth.cpp -o ground_truth -lfftw3f

执行生成可执行文件，执行结果如下，fftw3库被正常调用，且改写成功。

6、使用onemklfft二维实数到复数傅立叶变换，并将变换结果打印输出。

mkl_real2c.cpp：

#include <stdio.h>
#include <mkl.h>
#include "mkl_dfti.h"

#define COLS  4
#define ROWS  4

int main() {
    printf("Before Transform:\n");
    float *in = (float *)mkl_malloc(sizeof(float) * ROWS * COLS, 32);
    for (int i = 0; i < ROWS; i++) {
        for (int j = 0; j < COLS; j++) {
            in[i * ROWS + j] = i + j;
            printf("%7.2f",in[i * ROWS + j]);
        }
        printf("\n");
    }

    DFTI_DESCRIPTOR_HANDLE my_desc_handle = NULL;
    MKL_LONG status;
    MKL_LONG dim_sizes[2] = {ROWS, COLS};
    status = DftiCreateDescriptor(&my_desc_handle, DFTI_SINGLE, DFTI_REAL, 2, dim_sizes);
    DftiSetValue(my_desc_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
    DftiCommitDescriptor(my_desc_handle);
    MKL_Complex8 *out = (MKL_Complex8 *)mkl_malloc(sizeof(MKL_Complex8) * ROWS * COLS, 32);
    for(int i=0;i<ROWS*COLS;i++){
    	out[i].real=0;
    	out[i].imag=0;
    }
    DftiComputeForward(my_desc_handle, in, out);
    printf("After Transform:\n");
    for (int i = 0; i < ROWS; i++) {
        for (int j = 0; j < COLS; j++) {
            printf("%7.2f + %7.2fi\t", out[i * ROWS + j].real, out[i * ROWS + j].imag);
        }
        printf("\n");
    }

    // 释放内存和描述符
    mkl_free(in);
    mkl_free(out);
    DftiFreeDescriptor(&my_desc_handle);

    return 0;
}

经过编译链接执行后的输出结果为：

对比步骤4中fftw3结果，可以发现onemkl中的fft接口被已被正确调用，但是，值得注意的是，onemklfft为了节约空间开销，计算的结果利用傅立叶共轭对称的性质，只保存了完整结果的一半，因此在后面对比流程中，只对比了一半的结果。

7、两维RealtocomplexFFT参考代码实现

具体实现细节原理可参考realtocomplex_fft2d.cpp，经过编译链接运行，与fftw3结果简单对比，可以正确实现二维实数的FFT变换。

将以上各步整合汇总，加上性能对比，结果验证，可形成最终的解答。

3 题目解答

本次开发所用硬件配置为Intel® Core™ i5-8300H CPU @ 2.30GHz × 8, 16 GB RAM，所用操作系统是ubuntu 20.04 LTS，测试方法是使用参数卡设置相关参数，输出对比结果。

3.1 程序运行步骤

为了便于操作和观察结果，我们创建了一个名为parameter.txt的文件

参数卡parameter.txt界面：

--random data begin
  0.0
--random data end
  1.0
--iteration for comparing running time  
  1000
--comflag for chosing methods to compare(comflag=0 means comparision between RealtoComplexfft and FFTW3,comflag=1 means comparision between RealtoComplexfft and intelonemkl,comflag=2 means comparision between FFTW3 and intelonemkl)
  1
--precision allowed
  0.001

在这个文件中，我们定义了生成随机数范围的参数n0,n1,运行次数niter，需要对比的方法参数comflag以及允许的最大误差precision。可以在此文件中修改代码运行的相关参数，例如我们需要生成0.0到1之间的随机数，对比迭代1000次FFTW3与intelonemkl的性能，若它们之间结果最大误差不超过0.001则输出True，我们可以按上图所示定义参数。

运行此代码需要执行以下步骤：

1.定义parameter.txt文件，将其与onemkl_main.cpp放置在同一路径下。

2.<source /home/ganning/intel/oneapi/setvars.sh intel64>

3. <icpx onemkl_main.cpp -o onemkl_main -L/home/ganning/intel/oneapi/mkl/latest/lib -lmkl_rt>

4.<./onemkl_main>

3.5 相关代码附件

data_generate.cpp：随机数生成实例

ground_truth.cpp：fftw3计算二维实数fft转换实例

mkl_real2c.cpp：intelonemkl计算二维实数fft转换实例

realtocomplex_fft2d.cpp：Real to complex fft 参考代码

onemkl_main.cpp：汇总及性能对比代码

parameter.txt：定义onemkl_main.cpp参数文件

oneMkl.pdf:英特尔oneAPI黑客松---FFT算法的加速与优化说明文档

参考文献

ninggan99

关注

6
点赞
踩
4

收藏

觉得还不错? 一键收藏
2
评论
intel onemkl api 黑客松---FFT算法的加速与优化

基于Intel® Core™ i5-8300H CPU @ 2.30GHz × 8, 16 GB RAM，Ubuntu 20.04 LTS系统的intel onemkl api 下载安装及随机数生成、二维实数fft变换使用，对比了intelonemkl与fftw3以及自写fft的二维实数fft变换性能。
复制链接

扫一扫