oneMKL与FFTW3对FFT2D加速对比

最新推荐文章于 2024-08-26 13:30:20 发布

C语言狂暴战士

最新推荐文章于 2024-08-26 13:30:20 发布

阅读量461

点赞数

分类专栏： FFT 文章标签： c语言

本文链接：https://blog.csdn.net/m0_65775477/article/details/132679798

版权

FFT 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

Github上某大佬的文章指出，Intel下的onemkL函数库对于大型矩阵的计算很厉害，如下图，在运算性能以及速度上优于诸多方法，本文旨在通过自己亲身实践，学习不同算法函数库用以求解二维实数到复数的快速傅里叶变换，文章最后会附带上全部源代码。（这是本人的第一篇文章，欢迎大家共同交流学习，共勉。）

oneMKL和FFTW3介绍

oneMKL

英特尔 ®oneAPI 数学内核库（ Intel® oneAPI Math Kernel Library ，简称 oneMKL ）是一

个计算数学库，其中包含高度优化的广泛线程例程，适用于需要最高性能的应用程序。该库

提供 Fortran 和 C 编程语言接口。 oneMKLC 语言接口可以从用 C 或 C++ 以及可以引用 C 接

口的任何其他语言编写的应用程序中调用。

oneMKL 在这些主要计算领域提供全面的功能支持：

01、BLAS （1 级、 2 级和 3 级）和 LAPACK 线性代数例程，提供向量、向量矩阵和矩

阵矩阵运算。

02、ScaLAPACK 分布式处理线性代数例程，以及基本线性代数通信子程序 (BLACS) 和

并行基本线性代数子程序 (PBLAS) 。

03、oneMKL PARDISO （基于并行直接稀疏求解器 PARDISO* 的直接稀疏求解器），迭

04、代稀疏求解器，支持用于求解稀疏方程组的稀疏 BLAS （1 、 2 和 3 级）例程，以及

分布式版本提供 oneMKL PARDISO 求解器用于集群。

05、一维、二维或三维中的快速傅立叶变换 (FFT) 函数，支持混合基数（不限于 2 的幂

的大小），以及提供在集群上使用的这些函数的分布式版本。

06、用于优化向量数学运算的向量数学 (VM) 例程。

07、矢量统计 (VS) 例程，为多种概率分布、卷积和相关例程以及汇总统计函数提供高性

能矢量化随机数生成器 (RNG) 。

08、数据拟合库，提供基于样条的函数逼近、函数导数和积分以及搜索功能。

09、扩展特征求解器，基于 Feast 特征值求解器的特征求解器的共享内存编程 (SMP) 版

本。

英特尔 ® oneAPI 数学核心库 (oneMKL) 针对英特尔处理器上的性能进行了优化。

oneMKL 还可以在非 Intel x86 兼容处理器上运行。

oneMKL 可以作为单独的软件包下载：

https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html

也可以作为 Intel® oneAPI Base Toolkit 的一部分下载：

https://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html

进入到下载界面后，选择自己的操作系统以及在线或者离线安装方式之后，按照提示进行安

装即可。

FFTW3

FFTW3（Fastest Fourier Transform in the West）是一个开源的软件库，用于计算快速傅里叶变换（FFT）。它提供了高效的算法和实现，用于在数字信号处理、图像处理、模拟和科学计算等领域进行频率域分析。FFTW3被广泛应用于各种编程语言和平台上。它提供了多种不同的FFT变换算法，以适应不同的数据大小和性能需求。

实验环境配置

处理器：Intel i5

虚拟机：vwmare17，乌班图18.04（推荐使用，个人觉得比较稳定，出了问题网上可查的解决方法也较多，不过也不必纠结于版本，能跑就行）

代码实现

首先，调用oneMKL封装的函数，随机生成2048*2048大小的单精度浮点数：

VSLStreamStatePtr stream;
vslNewStream(&stream,VSL_BRNG_MT19937,1);
vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD,stream,Nx*Nz,input,0.0f,1.0f);
vslDeleteStream(&stream);

使用FFTW3计算real to imag FFT2D：

fftwf_complex *outputFFTW=(fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*Nx*(Nz/2+1));
fftwf_plan status=fftwf_plan_dft_r2c_2d(Nx,Nz,input,outputFFTW,FFTW_MEASURE);
t0=wallclock_time(); 
fftwf_execute(status);
t1=wallclock_time();

其中，wallclock_time()是自定义的时间函数，如下：

double wallclock_time(void)
{
	struct timeval s_val;
	static struct timeval b_val;
	double time;
	static int base=0;

	gettimeofday(&s_val,0);

	if (!base) {
		b_val = s_val;
		base = 1;
		return 0.0;
	}

	time = (double)(s_val.tv_sec-b_val.tv_sec) + 
		   (double)(1e-6*((double)s_val.tv_usec-(double)b_val.tv_usec));

	return (double)time;
}

使用oneMKL计算real to imag FFT2D：

MKL_LONG rs[3]={0,Nx,1};
MKL_LONG cs[3]={0,Nz/2+1,1};
MKL_Complex8 *outputMKL=(MKL_Complex8 *)malloc(sizeof(MKL_Complex8)*Nx*(Nz/2+1));
DFTI_DESCRIPTOR_HANDLE handle;
MKL_LONG sizeN[2]={Nx,Nz};	
DftiCreateDescriptor(&handle,DFTI_SINGLE,DFTI_REAL,2,sizeN);
DftiSetValue(handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE);
DftiSetValue(handle,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX);
DftiSetValue(handle,DFTI_INPUT_STRIDES,rs);
DftiSetValue(handle,DFTI_OUTPUT_STRIDES,cs);
DftiCommitDescriptor(handle);
t0=wallclock_time();
DftiComputeForward(handle,input,outputMKL);
t1=wallclock_time();

比较残差以及运算性能：

if(Nx==Nz)
		{	
			printf("-------------------------\n");
			diff=0;
			for(i=0;i<Nx*(Nz/2+1);i++)
			{
				diff1=(double)fabs(outputMKL[i].real-outputFFTW[i][0]);
				diff2=(double)fabs(outputMKL[i].imag-outputFFTW[i][1]);
				if(diff1>0.000001 || diff2>0.000001)
				{
					diff=1;
					break;
				}
			}
			if(diff==0) fprintf(stderr," CORRECT (residual < 0.000001) bwetween FFTW3 with oneMKL\n");
			else fprintf(stderr," ERROR (residual >0.000001 ) bwetween FFTW3 with oneMKL\n");
		}

最终结果：

可以发现，oneMKL精度在0.000001数量级下与开源的FFTW3相同，同时性能更加优秀。（不同电脑配置加速率不一样）

全部源代码

按照如下方法编译，这要求你的环境上正确安装编译器。

icx fft_2d_r2c_fftw_and_onemkl.c -o fft_2d_r2c_fftw_and_onemkl -qmkl -std=c99 -lmkl_rt -lm

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<sys/time.h>
#include<stdbool.h>
#include<string.h>
#include<sys/resource.h>
#include<fftw3.h>
#include<mkl.h>
#include"mkl_vsl.h"
#include"mkl_service.h"
#include"mkl_dfti.h"

size_t getCurrentMemoryUsage();
double wallclock_time(void);

int main(int argc,char *argv[])
{
	int i,Nx,Nz;
	double t0,t1,t,k,diff1,diff2,diff;
	size_t memoryUsage;
	Nx=2048;
	Nz=2048;
	//N2=16;
	k=log(Nz)/log(2)+1;
	
	while(Nz<=Nx)
	{
		float *input=(float *)malloc(sizeof(float)*Nx*Nz);
		/*Random numbers*/
    		VSLStreamStatePtr stream;
    		vslNewStream(&stream,VSL_BRNG_MT19937,1);
    		vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD,stream,Nx*Nz,input,0.0f,1.0f);
		vslDeleteStream(&stream);

		/*Use FFTW3 to do FFT2d*/
    		fftwf_complex *outputFFTW=(fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*Nx*(Nz/2+1));
    		fftwf_plan status=fftwf_plan_dft_r2c_2d(Nx,Nz,input,outputFFTW,FFTW_MEASURE);
    		t0=wallclock_time(); 
    		fftwf_execute(status);
    		t1=wallclock_time();
    		t=t1-t0;
    		fprintf(stderr,"   FFTW   Plan N=(%-5d  *  %-5d)   time is %f s\n",Nx,Nz,t);
		if(Nx==Nz)
		{
			t0=wallclock_time(); 
    			for(i=0;i<1000;i++) fftwf_execute(status);
    			t1=wallclock_time();
    			t=t1-t0;
    			fprintf(stderr,"   1000 times FFTW   Plan N=(%-5d  *  %-5d)   time is %f s\n",Nx,Nz,t);	
		}

		/*Use oneMKL to do FFT2d*/
    		MKL_LONG rs[3]={0,Nx,1};
    		MKL_LONG cs[3]={0,Nz/2+1,1};
    		MKL_Complex8 *outputMKL=(MKL_Complex8 *)malloc(sizeof(MKL_Complex8)*Nx*(Nz/2+1));
    		DFTI_DESCRIPTOR_HANDLE handle;
    		MKL_LONG sizeN[2]={Nx,Nz};	
    		DftiCreateDescriptor(&handle,DFTI_SINGLE,DFTI_REAL,2,sizeN);
    		DftiSetValue(handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE);
    		DftiSetValue(handle,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX);
    		DftiSetValue(handle,DFTI_INPUT_STRIDES,rs);
    		DftiSetValue(handle,DFTI_OUTPUT_STRIDES,cs);
    		DftiCommitDescriptor(handle);
    		t0=wallclock_time(); 
		DftiComputeForward(handle,input,outputMKL);
    		t1=wallclock_time();
    		t=t1-t0;
    		fprintf(stderr,"   oneMKL Plan N=(%-5d  *  %-5d)   time is %f s\n",Nx,Nz,t);
		if(Nx==Nz)
		{
			t0=wallclock_time(); 
    			for(i=0;i<1000;i++) DftiComputeForward(handle,input,outputMKL);
    			t1=wallclock_time();
    			t=t1-t0;
    			fprintf(stderr,"    1000 times oneMKL Plan N=(%-5d  *  %-5d)   time is %f s\n",Nx,Nz,t);
		}
		
		/*Compare oneMKL and FFTW3*/
		if(Nx==Nz)
		{	
			printf("-------------------------\n");
			diff=0;
			for(i=0;i<Nx*(Nz/2+1);i++)
			{
				diff1=(double)fabs(outputMKL[i].real-outputFFTW[i][0]);
				diff2=(double)fabs(outputMKL[i].imag-outputFFTW[i][1]);
				if(diff1>0.000001 || diff2>0.000001)
				{
					diff=1;
					break;
				}
			}
			if(diff==0) fprintf(stderr," CORRECT (residual < 0.000001) bwetween FFTW3 with oneMKL\n");
			else fprintf(stderr," ERROR (residual >0.000001 ) bwetween FFTW3 with oneMKL\n");
		}

		free(input);
		free(outputMKL);
		DftiFreeDescriptor(&handle);
		//fftwf_destory_plan(status);
		fftwf_free(outputFFTW);

		Nz=pow(2.0,k);
		k+=1.0;
	}

	memoryUsage = getCurrentMemoryUsage();
	fprintf(stderr,"  CPU usage size : %zu KB \n",memoryUsage);

	return 0;
}

size_t getCurrentMemoryUsage(){
	struct rusage usage;
	getrusage(RUSAGE_SELF, &usage);
	return usage.ru_maxrss;
}

double wallclock_time(void)
{
	struct timeval s_val;
	static struct timeval b_val;
	double time;
	static int base=0;

	gettimeofday(&s_val,0);

	if (!base) {
		b_val = s_val;
		base = 1;
		return 0.0;
	}

	time = (double)(s_val.tv_sec-b_val.tv_sec) + 
		   (double)(1e-6*((double)s_val.tv_usec-(double)b_val.tv_usec));

	return (double)time;
}