[学习]多线程

N体问题

模拟一个太空中有N个星体，各个星体有三个维度的速度和位置，相互之间受到万有引力影响运动的效果。模拟最小时间分位点为0.005秒，本题目中万有引力常数设定为1，星体的初始数据在“nbody.txt”文本文件中。要求用单线程和多线程分别完成一份，将迭代20次后的数据输出到文件中，并且比较串行程序和多线程之间的计算结果差异。

一、实现方案

模拟方案涉及大量计算，对于算法的考察较小，考研对数据结构和计算方式设计，以及对与多线程openmp和pthread库的运用。

• 实验环境：ubuntu20.04（Hype-V）
• 编译器版本：g++ 9.3.0

1、输入数据

一个txt文件，其中包括1024个星体的7项数据，一行为一个星体，分别是星体质量，星体在x、y、z上的位置，星体在x、y、z上的速度。

对于一个星体需要存放的数据如下，主要是x、y、z空间位置，sx、sy、sz速度和mass星体质量，还有计算时一些临时变量，如计算1号星体和其他星体相互作用力总和的变量，引入这个变量还有一个作用，可以减少误差（主要是指除法的误差累积）。星体结构体数据如下：

struct planet // 星体结构体
{
	double mass;	// body质量
	double sx, sy, sz;	// body三维方向上的速度
	double x, y, z;		// body三维方向上的位置
	double fx, fy, fz;	// body三维方向上收到的力的总和
	planet()			// 把所有数据初始化为0
	{
		mass = sx = sy = sz = x = y = z = 0;
		fx = fy = fz = 0;
	}
};

2、计算引力（单线程）

计算一个星体收到其他星体力的总和，无非两两之间计算力然后更新。 例如1号星体和2号星体、1号和3号、1号和4号……1号和1024号，2号和3号、2号和4号……等等（注意，这里并不需要求2号和1号，因为在计算1号和2号时已经计算了力大小，只需要取反）。因为这个模块也会在后来的多线程中被调用，因此抽象成为一个单独的函数。为了减小误差，将每次计算的力存放到临时变量中，并且在更新完所有的合力之后，再更新星体的速度和位置。 不要一上来就写并行化，拿单线程的程序验证算法的正确性是一个好习惯。

void speed_change(planet &a, planet &b){...}
void series_stimulus(unsigned int count, planet *array)				// 单线程模块，count表示迭代多少次，array是计算使用的元数据。
{
	while (count--)
	{
		for (int i = 0; i < max_size; i++)
		{
			for (int j = i + 1; j < max_size; j++)
			{
				speed_change(array[i], array[j]);
			}
		}
		for (int i = 0; i < max_size; i++)
		{
			postion(array[i]);
		}
	}
}

3、多线程优化

a、openmp多线程

有了上面单线程的计算方法，多线程就是把单线程的for，while等循环体并行化而已，openmp！！！简直没有比这个更加快捷的方法。逐个讲解pragma omp parallel带的参数。

• for : 表明并行化循环
• num_threads(n) : 并行线程数量，n可以是任意整数，其实还是看配置
• private(a,b,…,z) : 表明哪些变量是私有变量，不同线程之间不会相互干扰
• shared(a,b,…,z) : 表明哪些变量是公有变量，不同线程之间共享，例如本题中星体数据就是共享的，并且能够更新

void parallel_stimulus1(unsigned int count, planet *array)	// openmp骞惰绋嬪簭
{
	int i, j;
	while (count--)
	{
#pragma omp parallel for num_threads(4) private(i, j) shared(array)
		for (i = 0; i < max_size; i++)
		{
			for (j = i + 1; j < max_size; j++)
			{
				speed_change(array[i], array[j]);
			}
		}
		for (i = 0; i < max_size; i++)
		{
			postion(array[i]);
		}
	}
}

当然编译带有openmp的cpp文件需要在编译指令加入 -fopenmp 的参数，比如说你的文件叫做“a.cpp”那么编译参数就是：

g++ a.cpp -o a.out -fopenmp

b、pthread多线程

如果说上面的openmp是自动多线程，那么pthread就是手动多线程。使用pthread创建新线程之前必须要为新线程分配空间，thread_handle就是这样一个容纳新线程的容器，然后使用new为它分配空间（这里因为笔者笔记本只有4个核心，因此只创建4线程）。有了容器之后要装东西，不然一个空的线程容器有什么用呢。调用pthread_create函数创建新线程：

• 第一个参数是容器位置（也就是我们希望把这个新线程放在容器哪个位置）
• 第二个参数在这里没有用指定为NULL
• 第三个参数是函数指针（即希望在这个容器内运行什么），必须是void*类型
• 第四个参数是给函数使用的参数（这个容器内的函数总不可能完成无米之炊吧，所以米在哪里在这个位置告诉它）

pthread不会帮我们实现线程同步，所以需要手动实现线程同步。这里使用pthread_join函数，将之前创建的多个线程合并，统一更新星体的速度和位置。

void parallel_stimulus2(unsigned int count, planet *array)		// pthread主线程函数
{
	int num_of_thread = 4;
	thread_set *set = new thread_set[num_of_thread];
	pthread_t *thread_handle = new pthread_t[num_of_thread];
	while (count--)
	{
		for (long i = 0; i < num_of_thread; i++)
		{
			set[i].rank = i;
			set[i].o = array;
			pthread_create(&thread_handle[i], nullptr, thread_pthread, (void *)&set[i]);
		}
		for (long i = num_of_thread-1; i >= 0; i--)
		{
			pthread_join(thread_handle[i], nullptr);
		}
		
		for (int i = 0; i < max_size; i++)
		{
			postion(array[i]);
		}
	}
	delete[] thread_handle;
	delete[] set;
	return;
}

然后是容器内容物要怎么办，因为要手动并行化，所以这个线程做什么也需要自己计算。本处举个🌰是4线程，每个线程计算1024个星体的四分之一。为了传递参数方便，使用结构体thread_set，其一是方便、其二是因为pthread_create第四个参数是一个指针（而不是一个指针的指针），不能传递多个参数。

struct thread_set // pthread参数表
{
	int rank;		// 指定的线程号
	planet *o;		// 计算的元数据
};
void *thread_pthread(void *argv)							// pthread子线程函数
{
	thread_set set = *(thread_set *)argv;
	long rank = set.rank;
	int n = 0;
	
	if (rank == 3)
	{
		n = max_size - 1;
	}
	else
	{
		n = (rank + 1) * max_size / 4;				// 如果修改线程数，这里一并修改。
	}
	for (int i = rank * max_size / 4; i < n; i++)
	{
		for (int j = i + 1; j <= max_size; j++)
		{
			speed_change(set.o[i], set.o[j]);
		}
	}
	return nullptr;
}

为了能够使用pthread，需要在编译参数上增加“-lpthread”是“l”不是“i”。🌰同上。

g++ a.cpp -o a.out -lpthread

4、输入输出

需要从文件中读入计算数据，然后将迭代结果输出到文件。此处比较推荐的是ofstream头文件提供的文件流功能。能够像cout、cin一样使用。

a、文件读入

使用前建立一个fstream文件流指针，调用open方法打开需要读入的文件。然后像使用cin一样使用fp即可。

#include <fstream>

void read(planet *array)				
{
	fstream fp;
	fp.open("nbody.txt", ios::in);
	if (fp.bad())					// 检查文件是否存在
	{
		cout << "error to open file";
	}
	for (int i = 0; i < max_size; i++)
	{
		fp >> array[i].mass >> array[i].x >> array[i].y >> array[i].z >> array[i].sx >> array[i].sy >> array[i].sz;
	}
	fp.close();				// 读写结束后关闭文件
}

b、文件输出

一样，fstream也可以像cout一样使用，向文件写内容。

void write(planet *array, const char *file)			
{
	fstream fp;
	fp.open(file, ios::out);
	if (fp.bad())
	{
		cout << "output open error";
	}
	for (int i = 0; i < max_size; i++)
	{
		fp << setprecision(15) << setiosflags(ios::left) << array[i].mass << " " << array[i].x << ' '
		   << array[i].y << ' ' << array[i].z << ' '
		   << array[i].sx << ' ' << array[i].sy << ' ' << array[i].sz << endl;
	}
	fp.close();
}

5、main函数

int main()
{
	double start, end;
	planet array[max_size + 16];
	planet copy_planet[max_size + 16];
	read(array);

	unsigned int count = 20; //迭代20次
	char series[] = "series version.txt", parallel1[] = "parallel version1.txt", parallel2[] = "parallel version2.txt";	// 输出文件名

	memcpy(copy_planet, array, max_size * sizeof(planet)); // 串行程序
	GET_TIME(start);
	series_stimulus(count, copy_planet);
	GET_TIME(end)
	cout << (double)(end - start) << endl;
	// 要写出如文件请去掉本行注释，下面同理
	// write(copy_planet, series);

	memcpy(copy_planet, array, max_size * sizeof(planet)); // openmp程序
	GET_TIME(start);
	parallel_stimulus1(count, copy_planet);
	GET_TIME(end)
	cout << (double)(end - start) << endl;
	// write(copy_planet, parallel1);

	memcpy(copy_planet, array, max_size * sizeof(planet)); // pthread程序
	GET_TIME(start);
	parallel_stimulus2(count, copy_planet);
	GET_TIME(end)
	cout << (double)(end - start) << endl;
	// write(copy_planet, parallel2);
}

二、总结

完整代码在这个地方连接。