英特尔oneAPI黑客松 基于oneMKL的FFT加速

Windows平台环境配置过程

1. 前往Visual Studio官网安装Visual Studio Community，安装时选择使用C++进行开发。
2. 前往oneMKL下载页面下载适合的安装器并安装。在已经安装好VS的情况下，运行安装器后应该会出现Integrate with IDE字样，如下图
3. 在VS中新建C++控制台项目。进入项目后在项目属性页中启用oneMKL，如下图。

二维Real to complex FFT的实现

基于递归的简单FFT2D实现如下

void fft(complex<float>* a, int n) {
    if (n <= 1) {
        return;
    }
    complex<float>* even = new complex<float>[n / 2];
    complex<float>* odd = new complex<float>[n / 2];
    for (int i = 0; 2 * i < n; i++) {
        even[i] = a[2 * i];
        odd[i] = a[2 * i + 1];
    }
    fft(even, n / 2);
    fft(odd, n / 2);
    double angle = -2 * pi / n;
    complex<float> w(1), wn(cos(angle), sin(angle));
    for (int i = 0; 2 * i < n; i++) {
        a[i] = even[i] + w * odd[i];
        a[i + n / 2] = even[i] - w * odd[i];
        w *= wn;
    }
    delete[] even;
    delete[] odd;
}

void fft2D(complex<float>* fft_result, float* r2c_data, int rows, int cols) {
    for (int i = 0; i < num_points; i++) {
        fft_result[i] = complex<float>(r2c_data[i], 0.0);
    }
    for (int i = 0; i < rows; i++) {
        fft(fft_result + i * cols, cols);
    }
    complex<float>* column = new complex<float>[rows];
    for (int i = 0; i < cols; i++) {
        for (int j = 0; j < rows; j++) {
            column[j] = fft_result[j * cols + i];
        }
        fft(column, rows);
        for (int j = 0; j < rows; j++) {
            fft_result[j * cols + i] = column[j];
        }
    }
    delete[] column;
}

基于oneAPI的随机数生成以及FFT2D计算

相较于官方提供的实例代码，为了方便后续速度比较，这里使用了

DftiSetValue(mkl_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)

来避免生成的随机数数组被覆盖。由于上文所实现的fft2d函数数据存储格式为未经优化的复数-复数存储，为了方便后续正确性验证，这里使用了下述函数

DftiSetValue(mkl_handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)

来指定oneMKL FFT2D计算结果的存储格式。

void random_data_gen(float* r2c_data) {
    VSLStreamStatePtr stream;
    vslNewStream(&stream, VSL_BRNG_MT19937, 114514);
    vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD, stream, num_points, r2c_data, 0.0, 1.0);
}

void one_fft(float* r2c_data, complex<float>* mkl_result, int dims[]) {
    DFTI_DESCRIPTOR_HANDLE mkl_handle = NULL;
    DftiCreateDescriptor(&mkl_handle, DFTI_SINGLE, DFTI_REAL, 2, dims);
    DftiSetValue(mkl_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
    DftiSetValue(mkl_handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
    DftiCommitDescriptor(mkl_handle);
    DftiComputeForward(mkl_handle, r2c_data, mkl_result);
}

结果比对

经过上述设置后，oneMKL FFT2D结果存储格式如下图有色色块部分（8 * 8）。故对于一个N * N的矩阵而言，比对范围为N * (N / 2 + 1)。

函数如下

void comparison(complex<float>* fft_result, complex<float>* mkl_result) {
    int col = NUM /2 + 1, row = NUM, index = 0;
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < col; j++) {
            index = i * NUM + j;
            float dist = abs(fft_result[index] - mkl_result[index]);
            if (dist > 1e-4) {
                printf("结果不正确\n");
                return;
            }
        }
    }
    printf("结果正确\n");
}

结果展示

结果显示，oneMKL的时间优势是十分显著的。

结语

在比赛过程中，我了解了oneAPI的威力以及oneMKL工具的巨大潜力，也十分期待在未来的工作中继续进行oneAPI的应用。这次实现FFT2D算法的过程加深了我对傅里叶变换的记忆，而与oneMKL的对比结果也让我感受到了优化和加速工程问题的重要性。