高性能计算实验——基于Pthreads并行实现通用矩阵乘法、数组求和及二次方程组求解

1.实验目的

1.1.通过Pthreads实现通用矩阵乘法

熟练掌握Pthreads多线程编程方法，并将通用矩阵乘法转为多线程并行实现，进一步加深Pthreads的使用与理解。

1.2.基于Pthreads的数组求和

进一步熟悉Pthreads的写法，实现Pthreads的数组求和，加强对Pthreads的应用。

1.3.Pthreads求解二次方程组的根

进一步熟悉Pthreads的写法，实现Pthreads的二次方程组的根计算，加强对Pthreads的应用。

1.4.编写一个Pthreads多线程程序来实现基于monte-carlo方法的y=x^2阴影面积估算

通过复杂场景的应用进一步熟悉Pthreads的写法，实现Pthreads的面积计算，加强对Pthreads的应用。

2.实验过程和核心代码

2.1.通过Pthreads实现通用矩阵乘法

2.1.1 问题描述
通过 MPI 实现通用矩阵乘法（实验1）的并行版本，MPI并行进程（rank size）从 1 增加至 8，矩阵规模从 512 增加至 2048.
通用矩阵乘法（GEMM）通常定义为：
C = AB
输入：M , N, K 三个整数（512 ~2048）
问题描述：随机生成 MN 和NK 的两个矩阵 A,B,对这两个矩阵做乘法得到矩阵 C.
输出：A,B,C 三个矩阵以及矩阵计算的时间

2.1.2 实现过程
这里比MPI的实现更为简单，根据多线程的思想，只需要写好工作函数，就可以直接生成子线程来运行它们，如果将A、B函数设为全局，那么无需任何通信即可实现。
根据前一个实验的经验，这里我们可以对B函数进行转置，这样即可行访问（空间局部性优化）而不是列访问B了。

2.1.3.核心代码
①矩阵转置

这里将A、B设为全局函数，且我们的A、B、C矩阵采用一维数组的方式实现，这样可以方便的进行转置以及求C的相关操作。

②矩阵初始化

③工作函数（乘法函数）

运算时需要格外注意处理过的B矩阵相应的元素位置已经经过了转置；结果直接存储在全局C矩阵正确的位置，后续无需再处理。

④得到线程数并分配相应线程

⑤生成子线程

这里可以直接简单使用create、join两个函数完成任务的分发等待。

2.2.基于Pthreads的数组求和

2.2.1.问题描述
编写使用多个进程/线程对数组a[1000]求和的简单程序演示Pthreads 的用法。创建 n 个线程，每个线程通过共享单元
global_index 获取 a 数组的下一个未加元素，注意不能在临界段外访问全局下标 global_index
重写上面的例子，使得各进程可以一次最多提取 10 个连续的数， 以组为单位进行求和，从而减少对下标的访问

2.2.2.实现过程
这里题目的意思我的理解是编写的代码可以满足每次线程进入静态区后可以获得N个元素的求和权，这里的N是可变的。我也是这样实现的，通过bash输入N的方式控制线程每次进入临界区所能够获得的元素数。这里同时也可以控制线程的数目（同样可以通过bash）。
为了确保求和能够完成，工作线程需要设为while（1）的形式，当处理完毕后进行break退出。
已经完成统计的数组元素个数需要单独存储，并且它显然是一个临界区，每次只能有一个线程获取改变他，否则会出现重复统计的问题。

2.2.3.核心代码
①初始化锁


这里进入临界区需要锁，并且当前已经统计的元素个数、求总和sum都是需要锁的，因此初始化两个。

②获取线程数、每次线程进入临界区可以获得的元素数

还需要分配相应的线程。

③工作函数

这里为了提高灵活性，加快速度，在获得了相应的可以统计的元素坐标（begi）后将已统计的元素数（index）加上该进程将统计的数据后即可释放锁了，这样该进程在进行统计时不会影响到其他线程的工作。使用两个锁同样是这样考虑的。

④线程工作分发与释放

主线程额外输出时间的工作

2.3.Pthreads求解二次方程组的根

2.3.1.问题描述
编写一个多线程程序来求解二次方程组𝑎𝑥2 + 𝑏𝑥 + 𝑐 = 0的根，使用下面的公式

中间值被不同的线程计算，使用条件变量来识别何时所有的线程都完成了计算

2.3.2.实验过程
这个实验不是太明白，根据老师和助教的解释直接写最简单的版本，将其中的计算分给几个线程做，最后由主线程得到结果。
这里我将2a、b*b、4ac、 的计算分给了四个线程，其中最后一个计算需要等待前面的线程计算出结果才可以运算，我使用了最简单的停等来实现。

2.3.3.核心代码
①输入a、b、c

②工作函数

这里两个工作函数分别计算2a、b*b、4ac,为了使计算开方的线程可以判断是否这两步已经完成，特别使用一个标志flag，当flag=2时两个函数执行结束。

Kf使用轮询的方式看flag是否满足条件，满足则继续求开方。

③分发任务

这里将四个任务分发给四个线程。

④主线程得到结果

首先判断是否有解（whiq由kf函数进行赋值），有则计算出来，否则返回没有结果。

2.4.编写一个Pthreads多线程程序来实现基于monte-carlo方法的y=x^2阴影面积估算

核心思想为每一个进程都进行概率计算，即生成随机数在y=x^2曲线下方的次数 / 总次数

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <math.h>
 
#define intervals 10000     //区间数
int tnum=10             //进程数
int sum = 0;
pthread_mutex_t mutex=PTHREAD_MUTEX_INITIALIZER;
//计算区间
void* cast (void* args)
{
	int i;
	float x,y;
	int temp=0;
	for(i=0; i< (intervals/tnum); i++){
		x= (float)(rand() % 1001) * 0.001f;
		y= (float)(rand() % 1001) * 0.001f;
		if(y<x*x) temp++;
	}
    //加上
	pthread_mutex_lock(&mutex);
	sum += temp;
	pthread_mutex_unlock(&mutex);
}
 
int main(int argc, char *argv[]){
	printf("Please input the number of thread:\n");
    scanf("%d",&tnum);
	pthread_t thread[tnum];
	int i;
	for (i=0;i<tnum;i++)
	{
		if(pthread_create(&thread[i], NULL, cast, NULL))
			perror("Pthread Create Fails");
	}
	for (i=0;i<tnum;i++)
	{
		if(pthread_join(thread[i], NULL))
			perror("Pthread Join Fails");
	}
	float shadow=(float) sum/(float) intervals;
	printf("sum:%d\n",sum);
	printf("intervals:%d\n",intervals);
	printf("the area of the shadow is:%f\n",shadow);
	return 0;
}

3.实验结果

3.1.通过Pthreads实现通用矩阵乘法

执行得到N=2、M=3，K=4大小的矩阵乘法使用2个线程结果正确！
After B为B转置后的矩阵形式，正确！

执行得到2048大小的矩阵使用2个线程时间消耗为12s
这里比MPI快的原因就是通信量的差异以及B矩阵行访问与列访问之间的区别了，比MPI快很多。

执行得到2048大小的矩阵使用4个线程时间消耗为7.8s
结果可以看出在处理器足够的情况下，增加线程数目可以大大提高运算速度，符合预期。

3.2.基于Pthreads的数组求和

这里使用两个线程，每个线程每次仅能统计1个元素，最后得到的sum正确！

这里使用四个线程，每个线程每次仅能统计1个元素，最后得到的sum正确！

这里使用四个线程，每个线程每次能统计10个元素，符合第二个要求，最后得到的sum正确！

这里使用八个线程，每个线程每次能统计10个元素，符合第二个要求，最后得到的sum正确！

3.3.Pthreads求解二次方程式的根

首先输入1 4 4作为系数，得到的结果正确！

输入1 6 9作为系数，得到的结果正确！

输入1 7 6作为系数，得到的结果正确！

输入1 3 3作为系数，得到的结果正确！

3.4.编写一个Pthreads多线程程序来实现基于monte-carlo方法的y=x^2阴影面积估算

使用四个线程得到的结果接近正确结果，有效的统计数3359，结果0.3359

使用四个线程得到的结果接近正确结果，有效的统计数3286，结果0.3286

实验感想

本次实验是针对Pthreads编程，主要是使用Pthreads进行矩阵乘法的实现，总体来说并不难。
通过MPI与Pthreads两种方式实现矩阵乘法可以发现，MPI中通信量是占时间最大部分的，同时相应的内存分配与回收都是开销很大的操作，此外行访问较列访问有着空间局限性上的巨大优势。