OpenMP实现并行快速傅里叶变换

发表于 17 Jun 2019 总计 4257字 阅读时间 15分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

---

题目

Final Project, Selected topics (not limited to):

1. Parallel Monte Carlo algorithm;
2. Parallellinearsolver;
3. Parallel QR factorization or linear least squares;
4. Parallelbenchmarks;
5. Dynamicloadbalancing;
6. Parallel sorting or selection;
7. Parallelgraphpartitioning;

Implement with MPI or MapReduce or OpenMP or CUDA …

说明

使用OpenMP实现了并行版本的快速傅里叶变换，并与串行版本的傅里叶变换结果作比较，验证其正确性并对比运行时间。

分析

傅里叶变换常用于加速多项式乘法，而常规的快速傅里叶变换（原理略）通常是使用递归实现的，使用并行优化的难度比较高。因此，我在这里实现了非迭代的快速傅里叶版本：先预处理每个位置上元素变换后的位置（每个位置分治后的最终位置为其二进制翻转后得到的位置），然后先将所有元素移到变换后的位置之后直接循环合并。

找变换位置这里其实有一个经典算法叫雷德算法，又被称作蝴蝶变换；不过我没有使用这一算法，因为蝴蝶变换有一定的循环依赖性，很难并行优化。

随后，调整完循环顺序之后，第一层循环变量i表示每一层变换的跨度，第二层循环变量j表示每一层变换的第一个起点，第三层循环遍历k则表示实际变换的位置$k$和$k+i$。在这里，从第二层开始是没有循环依赖的，即对于不同的j，算法不会对同一块地址进行访问（因为访问的下标$k\equiv j\pmod{i}$且$k+i\equiv j\pmod{i}$）。

作为公平起见，用于对比的串行版本快速傅里叶变换是直接在并行版本上删去编译推导#pragma omp for得到的。这是因为递归版本的快速傅里叶变换通常有较大的函数递归开销。

编译

g++ ParallelFastFourierTransform.cpp -oParallelFastFourierTransform -fopenmp

运行

接受两个运行参数，第一个参数是并行部分使用的线程数量，第二个参数是傅里叶变换的长度的指数（例如，当设为10的时候，傅里叶变换的长度就是$2^{10}=1024$；这样做是因为算法本身就要求长度为二的幂次）。

考虑到变换的长度可能很长，我这里直接使用程序自己生成的序列进行变换了（这个程序仅用于验证算法的正确性）。

./ParallelFastFourierTransform <num-of-threads> <power-of-transform-length>

结果

我的软硬件配置是：

• VAIO Z Flip 2016
• Intel(R) Core(TM) i7-6567U CPU @3.30GHZ 3.31GHz
• 8.00GB RAM
• Windows 10.0.18362.175, 64-bit
• Windows Subsystem for Linux [Ubuntu 18.04.2 LTS]：WSL是以软件的形式运行在Windows下的Linux子系统，是近些年微软推出来的新工具，可以在Windows系统上原生运行Linux。
• gcc version 7.3.0 (Ubuntu 7.3.0-27ubuntu1~18.04)：C语言程序编译器，Ubuntu自带。

$ ./ParallelFastFourierTransform 1 22
Serial Time: 4.29194s
Parallel Time: 4.32889s
$ ./ParallelFastFourierTransform 2 22
Serial Time: 4.30011s
Parallel Time: 2.9339s
$ ./ParallelFastFourierTransform 4 22
Serial Time: 4.26301s
Parallel Time: 1.87973s
$ ./ParallelFastFourierTransform 8 22
Serial Time: 4.27979s
Parallel Time: 1.89624s

可以看到，由于我的机器是双核心四线程，在四线程并行加速时得到了最优的加速比。四线程并行优化的傅里叶变换比串行版本快了约56%.

源代码

#include <bits/stdc++.h>
#include <omp.h>
using namespace std;
typedef long long ll;
typedef double lf;
struct Rader : vector<int>
{
	Rader(int n) : vector<int>(1 << int(ceil(log2(n))))
	{
		for (int i = at(0) = 0; i < size(); ++i)
			if (at(i) = at(i >> 1) >> 1, i & 1)
				at(i) += size() >> 1;
	}
};
struct FFT : Rader
{
	vector<complex<lf>> w;
	FFT(int n) : Rader(n), w(size(), polar(1.0, 2 * M_PI / size()))
	{
		w[0] = 1;
		for (int i = 1; i < size(); ++i)
			w[i] *= w[i - 1];
	}
	vector<complex<lf>> pfft(const vector<complex<lf>> &a) const
	{
		vector<complex<lf>> x(size());
#pragma omp parallel for
		for (int i = 0; i < a.size(); ++i)
			x[at(i)] = a[i];
		for (int i = 1; i < size(); i <<= 1)
#pragma omp parallel for
			for (int j = 0; j < i; ++j)
				for (int k = j; k < size(); k += i << 1)
				{
					complex<lf> t = w[size() / (i << 1) * j] * x[k + i];
					x[k + i] = x[k] - t, x[k] += t;
				}
		return x;
	}
	vector<complex<lf>> fft(const vector<complex<lf>> &a) const
	{
		vector<complex<lf>> x(size());
		for (int i = 0; i < a.size(); ++i)
			x[at(i)] = a[i];
		for (int i = 1; i < size(); i <<= 1)
			for (int j = 0; j < i; ++j)
				for (int k = j; k < size(); k += i << 1)
				{
					complex<lf> t = w[size() / (i << 1) * j] * x[k + i];
					x[k + i] = x[k] - t, x[k] += t;
				}
		return x;
	}
	vector<ll> ask(const vector<ll> &a, const vector<ll> &b) const
	{
		vector<complex<lf>> xa(a.begin(), a.end()), xb(b.begin(), b.end());
		xa = fft(xa), xb = fft(xb);
		for (int i = 0; i < size(); ++i)
			xa[i] *= xb[i];
		vector<ll> ans(size());
		xa = fft(xa), ans[0] = xa[0].real() / size() + 0.5;
		for (int i = 1; i < size(); ++i)
			ans[i] = xa[size() - i].real() / size() + 0.5;
		return ans;
	}
};
int main(int argc, char **argv)
{
	if (argc < 3)
		return cerr << "Error: No Enough parameters (" << argv[0] << " <num-of-threads> <power-of-transform-length>).\n", 0;
	omp_set_num_threads(atoi(argv[1]));
	FFT fft(1 << atoi(argv[2]));
	vector<complex<lf>> a(fft.begin(), fft.end());
	double t0 = omp_get_wtime();
	vector<complex<lf>> b = fft.fft(a);
	double t1 = omp_get_wtime();
	cout << "Serial Time: " << t1 - t0 << "s\n";
	vector<complex<lf>> c = fft.pfft(a);
	double t2 = omp_get_wtime();
	cout << "Parallel Time: " << t2 - t1 << "s\n";
	if (b != c)
		cerr << "Error: Parallel result are not equivalent to Serial result.\n";
}