英特尔 oneAPI 黑客松——使用 oneMKL 对 FFT 算法进行加速与优化

最新推荐文章于 2023-09-05 19:59:49 发布

Miraitowa65

最新推荐文章于 2023-09-05 19:59:49 发布

阅读量541

点赞数

文章标签：算法 oneapi

本文链接：https://blog.csdn.net/A15249/article/details/132426158

版权

英特尔 oneAPI 黑客松——使用 oneMKL 对 FFT 算法进行加速与优化

快速傅里叶变换 FFT
使用 oneMKL 对 FFT 算法进行加速和优化
使用oneMKL的一些感受
完整代码

快速傅里叶变换 FFT

傅里叶变换

傅里叶变换（Fourier transform，简称 FT）是一种线性积分变换，用于函数（应用上称作“信号”）在时域和频域之间的变换。因其基本思想首先由法国学者约瑟夫·傅里叶系统地提出，所以以其名字来命名以示纪念。
傅里叶变换是一种在信号处理、图像处理、通信、物理学等领域中广泛应用的数学工具和技术。它将一个函数（或信号）从时域（时间域）转换为频域（频率域），揭示了信号在不同频率下的成分和特性。傅里叶变换有两种常见的形式：连续傅里叶变换（Continuous Fourier Transform）和离散傅里叶变换（Discrete Fourier Transform）。

（连续）傅里叶变换的定义如下：
连续傅里叶变换
其中，𝑓为原函数，𝑓̂是原函数的傅里叶变换，𝑥表示时间，其定义域为时域，𝜉表示频率，其定义域为频域，𝑖是虚数单位。通常情况下，𝑓是一个实函数，而𝑓̂则是一个复数值函数。

离散傅里叶变换

离散傅里叶变换（Discrete Fourier Transform，简称 DFT），是傅里叶变换在时域和频域上都离散的形式，将信号的时域采样变换为频域采样。
对于𝑁点序列{𝑥[𝑛]}_{0≤𝑛<𝑁}，它的傅里叶变换（DFT）是一个同样长度为𝑁的序列：
离散傅里叶变换
其中𝑒为自然底数，𝑖是虚数单位。

快速傅里叶变换

快速傅里叶变换（Fast Fourier Transform, FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法。FFT 能够将计算 DFT 的复杂度从𝑂(𝑛²)降低到 𝑂(𝑛𝑙𝑜𝑔(𝑛))，其中𝑛为数据大小。
FFT 算法的关键思想是利用信号的对称性和周期性，通过分而治之的策略来加速傅里叶变换的计算，从而大幅提高了计算率。通过将信号分解成多个子信号，然后逐步合并子信号的频谱，从而实现对整体频谱的计算。

使用 oneMKL 对 FFT 算法进行加速和优化

Intel oneMKL（Intel oneAPI Math Kernel Library）是一款高性能数学核心库，旨在为开发人员提供丰富的数学函数和优化算法，以加速各种科学、工程和数据分析应用的性能。oneMKL 是 Intel oneAPI 套件的一部分，适用于多种硬件架构，包括 Intel CPU、GPU 以及其他加速器。

问题简介

调用 oneMKL 相应 API 函数，产生 2048 * 2048 个随机单精度实数()；
根据 1 产生的随机数据作为输入，实现两维 Real to complex FFT 参考代码；
根据 1 产生的随机数据作为输入，调用 oneMKL API 计算两维 Real to complex FFT；
结果正确性验证，对 2 和 3 计算的两维 FFT 输出数据进行全数据比对（允许适当精度误差），输出 “结果正确”或“结果不正确”信息；
平均性能数据比对（比如运行 1000 次），输出 FFT 参考代码平均运行时间和 oneMKL FFT 平均运行时间。

测试准备

在 Intel 官网上下载安装 oneAPI 后，将工具包集成到 VS2022 中。创建空项目，进入 VS2022 中的项目属性页面，在 VC++目录中配置可执行文件目录，包含目录以及库目录三个路径。具体配置路径如下：
可执行目录：C:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\bin\intel64
包含目录：C:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\include
库目录：C:\Program Files (x86)\Intel\oneAPI\compiler\2023.2.0\windows\compiler\lib\intel64_win；
C:\Program Files(x86)\Intel\oneAPI\mkl\2023.2.0\lib\intel64
在这里插入图片描述
然后在 win64 环境下配置 oneMKL 库，需要在项目中链接添加以下库作为附加依赖项：

	mkl_intel_ilp64.lib
	mkl_intel_thread.lib
	mkl_core.lib
	libiomp5md.lib

在 win64 环境下配置 fftw3 库，需要在项目中链接添加以下库作为附加依赖项：

	libfftw3-3.lib
	libfftw3f-3.lib
	libfftw3l-3.lib

由于 fftw3 接口是作为 oneMKL 的包装器开源提供的，fftw3 接口集成在 oneMKL 中。所以在 oneMKL 安装完成后可以直接调用，不需要添加上述的附加依赖项。
在这里插入图片描述

设置好附加依赖项后，需要在 VS2022 中集成的 oneAPI 工具包中将是否使用 oneMKL 选项设置成 Sequential。

具体测试过程

首先需要产生 2048 * 2048 个随机单精度实数，通过调用 oneMKL 中的随机数生成器 vslNewStream 创建和初始化数据流，然后将数据流作为输入即可。具体代码如下：

// 分配输入数组
float* inputData = (float*)malloc(N * M * sizeof(float));

// 生成随机单精度实数数据作为输入
VSLStreamStatePtr stream;
vslNewStream(&stream, VSL_BRNG_MT19937, 1);
vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD, stream, N * M, inputData, 0.0f, 1.0f);

然后分别调用 fftwf_plan_dft_r2c_2d 函数和 DftiCreateDescriptor 函数创建 FFTW 计划和 oneMKL 描述符，并分别为异地存储的输出分配空间。

// 创建 FFTW 计划
fftwf_complex* fftwOutputData = (fftwf_complex*)fftwf_malloc((N / 2 + 1) * M * sizeof(fftwf_complex));
fftwf_plan fftwPlan = fftwf_plan_dft_r2c_2d(N, M, inputData, fftwOutputData, FFTW_ESTIMATE);

// 创建 OneMKL 描述符
MKL_LONG rs[3] = { 0, N, 1 };
MKL_LONG cs[3] = { 0, N / 2 + 1, 1 };

DFTI_DESCRIPTOR_HANDLE hand = NULL;
DftiCreateDescriptor(&hand, DFTI_SINGLE, DFTI_REAL, 2, NI);
DftiSetValue(hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
DftiSetValue(hand, DFTI_INPUT_STRIDES, rs);
DftiSetValue(hand, DFTI_OUTPUT_STRIDES, cs);
DftiCommitDescriptor(hand);

// 分配 OneMKL 输出数组
MKL_Complex8* mklOutputData = (MKL_Complex8*)malloc((N / 2 + 1) * M * sizeof(MKL_Complex8));

分别执行 FFTW 计划和 oneMKL 描述符计算 FFT 变换的结果，并对每次的执行时间进行计时。

// FFTW FFT
LARGE_INTEGER fftwStart, fftwEnd;
QueryPerformanceCounter(&fftwStart);
fftwf_execute(fftwPlan);
QueryPerformanceCounter(&fftwEnd);
fftwTotalTime += (double)(fftwEnd.QuadPart - fftwStart.QuadPart);

// OneMKL FFT
LARGE_INTEGER mklStart, mklEnd;
QueryPerformanceCounter(&mklStart);
DftiComputeForward(hand, inputData, mklOutputData);
QueryPerformanceCounter(&mklEnd);
mklTotalTime += (double)(mklEnd.QuadPart - mklStart.QuadPart);

每次执行完成后，对 FFTW 和 oneMKL的输出进行比对，验证输出结果的正确性。具体比对方式是分别比较二者进行 FFT 计算后输出的复数的实部和虚部差异，预先设定精度误差允许范围为1e-6，超过该精度误差则认为输出结果不正确。统计输出的结果中不同的数据个数。

// 比较两个库的输出数据
int mismatchCount = 0;
for (int i = 0; i < (N / 2 + 1) * M; i++) {
    float diff_real = fabs(mklOutputData[i].real - fftwOutputData[i][0]);
    float diff_imag = fabs(mklOutputData[i].imag - fftwOutputData[i][1]);
    if (diff_real > 1e-6 || diff_imag > 1e-6) {
        mismatchCount++;
    }
}

设定程序循环1000次，每次循环都重新运行整个测试过程。即重新生成新的数据作为输入，并对每次循环的输出结果进行比对，统计是否存在结果不正确的次数，以达到排除结果正确的偶然性的目的。同时计算比较 FFTW 和 oneMKL 的平均运行时间。

此外，为了可视化1000次循环的循环进程，额外添加了一个进度条用于记录循环完成进度。

// 进度条可视化1000次循环进度
void printProgressBar(int progress, int total) {
    float percentage = (float)progress / total * 100;
    int barWidth = 50;
    int numBars = (int)(percentage / 100 * barWidth);

    printf("[");
    for (int i = 0; i < barWidth; i++) {
        if (i < numBars) {
            printf("=");
        }
        else {
            printf(" ");
        }
    }
    printf("] %.2f%%\r", percentage);
    fflush(stdout);
}

最终测试结果见下图：
在这里插入图片描述
可以看到 oneMKL 中计算的两维 FFT 的结果是正确的，同时oneMKL 库的计算速度明显比 fftw3 快。而且 oneMKL 中集成的 fftw3是经过 oneMKL 加速的，所以实际运行中 oneMKL 的 FFT 计算速度会比开源的 fftw3 库快很多。

使用oneMKL的一些感受

在实际使用中，oneMKL 提供了简洁且易于理解的接口，以便用户能够迅速上手。只需几行代码，就可以将一系列数据从时域转换到频域。这让用户能够更专注于应用程序的核心逻辑，而不必为数学计算的复杂性而烦恼。
另一个引人注目的方面是性能提升。使用 oneMKL 后，对大数据量的计算速度明显加快了。这是因为 oneMKL 利用了多线程和硬件优化，以提供更高的并行性和计算效率。
另外，oneMKL 提供的帮助文档和示例代码清晰实用。能够帮助用户更好地理解如何使用库中的函数，搭建正确的代码逻辑。

总的来说，使用 Intel oneMKL 进行 FFT 变换让我感受到了现代数学核心库的优越性。它能够帮助用户节省时间，提升了计算性能，并且为用户提供了更高效的开发体验。无论是在科学计算、数据分析还是其他领域，oneMKL 都是一个强大的工具，让用户能够更好地将精力集中在解决问题和创造价值上。

完整代码

#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
#include <mkl.h>
#include <math.h>
#include <windows.h>

#define NUM_ITERATIONS 1000

// 进度条可视化1000次循环进度
void printProgressBar(int progress, int total) {
    float percentage = (float)progress / total * 100;
    int barWidth = 50;
    int numBars = (int)(percentage / 100 * barWidth);

    printf("[");
    for (int i = 0; i < barWidth; i++) {
        if (i < numBars) {
            printf("=");
        }
        else {
            printf(" ");
        }
    }
    printf("] %.2f%%\r", percentage);
    fflush(stdout);
}

int main() {
    double fftwTotalTime = 0.0;
    double mklTotalTime = 0.0;
    int N = 2048, M = 2048;
    MKL_LONG NI[2] = { M, N };
    int flag = 0;

    for (int iter = 0; iter < NUM_ITERATIONS; iter++) {
        // 分配输入数组
        float* inputData = (float*)malloc(N * M * sizeof(float));

        // 生成随机单精度实数数据作为输入
        VSLStreamStatePtr stream;
        vslNewStream(&stream, VSL_BRNG_MT19937, 1);
        vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD, stream, N * M, inputData, 0.0f, 1.0f);

        // 创建 FFTW 计划
        fftwf_complex* fftwOutputData = (fftwf_complex*)fftwf_malloc((N / 2 + 1) * M * sizeof(fftwf_complex));
        fftwf_plan fftwPlan = fftwf_plan_dft_r2c_2d(N, M, inputData, fftwOutputData, FFTW_ESTIMATE);

        // 创建 OneMKL 描述符
        MKL_LONG rs[3] = { 0, N, 1 };
        MKL_LONG cs[3] = { 0, N / 2 + 1, 1 };

        DFTI_DESCRIPTOR_HANDLE hand = NULL;
        DftiCreateDescriptor(&hand, DFTI_SINGLE, DFTI_REAL, 2, NI);
        DftiSetValue(hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
        DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
        DftiSetValue(hand, DFTI_INPUT_STRIDES, rs);
        DftiSetValue(hand, DFTI_OUTPUT_STRIDES, cs);
        DftiCommitDescriptor(hand);

        // 分配 OneMKL 输出数组
        MKL_Complex8* mklOutputData = (MKL_Complex8*)malloc((N / 2 + 1) * M * sizeof(MKL_Complex8));

        // FFTW FFT
        LARGE_INTEGER fftwStart, fftwEnd;
        QueryPerformanceCounter(&fftwStart);
        fftwf_execute(fftwPlan);
        QueryPerformanceCounter(&fftwEnd);
        fftwTotalTime += (double)(fftwEnd.QuadPart - fftwStart.QuadPart);

        // OneMKL FFT
        LARGE_INTEGER mklStart, mklEnd;
        QueryPerformanceCounter(&mklStart);
        DftiComputeForward(hand, inputData, mklOutputData);
        QueryPerformanceCounter(&mklEnd);
        mklTotalTime += (double)(mklEnd.QuadPart - mklStart.QuadPart);

        // 比较两个库的输出数据
        int mismatchCount = 0;
        for (int i = 0; i < (N / 2 + 1) * M; i++) {
            float diff_real = fabs(mklOutputData[i].real - fftwOutputData[i][0]);
            float diff_imag = fabs(mklOutputData[i].imag - fftwOutputData[i][1]);
            if (diff_real > 1e-6 || diff_imag > 1e-6) {
                mismatchCount++;
            }
        }

        if (mismatchCount > 0) {
            printf("Iteration %d - 结果不正确！\n", iter + 1);
            flag += 1;
        }

        // 释放资源
        fftwf_destroy_plan(fftwPlan);
        fftwf_free(fftwOutputData);
        DftiFreeDescriptor(&hand);
        free(inputData);
        free(mklOutputData);

        printProgressBar(iter + 1, NUM_ITERATIONS);
    }

    printf("\n");

    if (flag == 0) {
        printf("1000次运行结果全部正确！\n");
    }

    // 输出比对结果和平均运行时间
    if (fftwTotalTime > 0.0 && mklTotalTime > 0.0) {
        double fftwAverageTime = fftwTotalTime / NUM_ITERATIONS;
        double mklAverageTime = mklTotalTime / NUM_ITERATIONS;

        printf("OneMKL 中集成的 FFTW3 FFT 平均运行时间: %f 微秒\n", fftwAverageTime);
        printf("OneMKL FFT 平均运行时间: %f 微秒\n", mklAverageTime);
    }

    return 0;
}

Miraitowa65

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
1
评论
英特尔 oneAPI 黑客松——使用 oneMKL 对 FFT 算法进行加速与优化

调用 oneMKL 相应 API 函数，产生 2048 * 2048 个随机单精度实数()；根据 1 产生的随机数据作为输入，实现两维 Real to complex FFT 参考代码；根据 1 产生的随机数据作为输入，调用 oneMKL API 计算两维 Real to complex FFT；结果正确性验证，对 2 和 3 计算的两维 FFT 输出数据进行全数据比对（允许适当精度误差），输出 “结果正确”或“结果不正确”信息；
复制链接

扫一扫