diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw1.md" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw1.md" new file mode 100644 index 000000000..e2c31fffe --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw1.md" @@ -0,0 +1,63 @@ +## 1.1 + +设总共$n$个数待求和,总共有$P$个核 + +对于编号为$p(0 \le p \lt P)$的核 + +$my\_first\_i = p * \lfloor \frac{n}{P} \rfloor$ + +$my\_last\_i=\begin{cases} + n & p = P-1 + \\ + (p+1) * \lfloor \frac{n}{P} \rfloor & else + \end{cases}$ + +## 1.6 + +a. 需执行接收和加法各$P-1$次 + +b. 需执行接收和加法各$\lceil \log_2 P \rceil$次 + +(这是因为b的次数满足$f(P)=\begin{cases} + 0 & P=1 + \\ + f(\lceil \frac{P}{2} \rceil) + 1 & else + \end{cases}$) + +表格如下: +| P | a | b | +| ------ | ------ | ------ | +| 2 | 1 | 1 | +| 4 | 3 | 2 | +| 8 | 7 | 3 | +| 16 | 15 | 4 | +| 32 | 31 | 5 | +| 64 | 63 | 6 | +| 128 | 127 | 7 | +| 256 | 255 | 8 | +| 512 | 511 | 9 | +| 1024 | 1023 | 10 | + +## 1.7 + +我认为这二者的成分都有。可以认为总共有$\lceil \log_2 P\rceil + 1$种任务(分别对应参与$0, 1, ..., \lceil \log_2 P \rceil$次接受/求和),分配给$P$个核心来处理,所以存在执行相同任务而数据不同的核心,也存在执行不同任务的核心。 + + +## 1.9 + +我曾经在计算机图形学代码中使用并行计算,具体来说是利用多线程加速光线追踪。大致伪代码可以表述如下 +``` +func trace(ray) { + iterate over scene, find hit point, compute radiance, do reflect/refract and continue +} +``` +每个线程负责确定一定范围内的像素的颜色,它向这些像素的位置发射光线并调用trace;除了中途输出进度之外,不需要任何同步机制;所有线程的计算任务完成后即可输出图像,不需要额外的数据整合过程。 + +我使用了线程库(Rust,C++)和OpenMp(C++)两种方式,二者的效率没有明显差距。 + +除此之外,我还尝试了使用GPU来加速运算。GPU的核心数目远远多于CPU,而且此问题本身特别适合GPU并行执行,所以可以很容易的得到更好的效果。 + +最后,在我的硬件配置(CPU:R7-2700,GPU:GTX 1060)上渲染smallpt的场景(只是场景相同,代码基本上没有任何关系,我自己还进行了很多优化;图片参数为1024px * 1024 px,10000spp),CPU单核大概耗时25min,CPU并行大概耗时2min40s,GPU并行大概耗时8s。 + +我认为这是典型的数据并行,因为它是在一列数据(像素点位置)上并行执行相同的算法(trace)。 + diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw2.md" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw2.md" new file mode 100644 index 000000000..dd34935ed --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw2.md" @@ -0,0 +1,76 @@ +## Q1 当讨论浮点数加法时,我们曾简单地假设每个功能单元都花费相同的时间。现在,如果每个取命令与存命令都耗费2纳秒,其余的每个操作耗费1纳秒。 +### a.在上述假设下,每个浮点数加法要耗费多少时间? +有两个耗时2ns的指令和5个耗时1ns的指令,总共耗时9ns + +### b.非流水线1000对浮点数加法要耗费多少时间? +9ns * 1000 = 9000ns + +### c.流水线1000对浮点数加法要耗费多少时间? +2ns * 1000 + 7ns = 2007ns + +### d. 如果操作数/结果储存在不同级的内存层级上,那么取命令与存命令所要耗费的时间可能会差别非常大。假设从一级缓存上取数据/指令要耗费2纳秒,从二级缓存上去数据/指令要耗费5纳秒,从主存取数据/指令要耗费50纳秒,当执行某条指令,取其中一个操作数时,发生了一次一级缓存失效,那么流水线会发生什么情况?如果又发生了二级缓存失效,又会怎样? +这里假定读取耗时已经包括了在上层的存储结构中的查找过程的时间 + +一级缓存失效时:5ns后下一次读取才能开始,而7ns后上一次加法才结束,因此任一时刻还是有>=2次加法在执行,可以认为流水线没有中断 + +二级缓存失效时:50ns后下一次读取才开始,因此存在43ns的时间CPU只在等待这一次读取,可以认为流水线在这里中断了 + +## Q2 回顾之前一个从缓存读取二维数组的例子(见第2.2.3小节)。请问一个更大的矩阵和一个更大的缓存是如何影响上面两段嵌套循环程序的性能的?如果MAX=8,缓存可以储存4个缓存行,情况又会怎样?在第一段嵌套循环程序中对A的读操作,会导致发生多少次读缺失?第二段嵌套循环程序的缺失次数又是多少? +## 注:一个缓存行可以放四个数组元素。 +假定缓存的行数增多,缓存行的大小不变,而矩阵的行列等比例增大,那么版本1缓存缺失的次数会增多,因为对于矩阵的每一行也会发生多次缓存缺失(大约$\frac{矩阵元素个数}{缓存行大小}$次);如果缓存的行数与矩阵的大小差不多,那么版本2的缓存缺失的次数可以减小至与版本1相当,否则如果缓存行数明显小于矩阵行数,则大约会发生$矩阵元素个数$次缓存缺失。 + +MAX = 8,缓存行数 = 4时,版本1会有16次缓存缺失,版本2会有64次缓存缺失。 + +## Q3 在表2.2中,虚拟地址由12位字节偏移量和20位的虚拟页号组成。如果一个程序运行的系统上拥有这样的页大小和虚拟地址空间,这个程序最多可以有多少页? +2 ^ 20 = 1048576个 + +## Q4 假设一个程序需要运行$10^{12}$条指令来解决一个特定问题,一个单处理器系统可以在$10^6$秒(大约11.6天)内完成。所以一个单处理器系统平均每秒运行$10^6$条指令。现在假设程序已经实现并行化,可以在分布式内存系统上运行。该并行程序使用p个处理器,每个处理器执行$10^{12}/p$条指令,且需要发送$10^9(p-1)$条消息。执行该程序时,不会有额外的开销,即每个处理器执行完所有指令并发送完所有消息之后,程序就结束了,而不会有诸如等待消息等时间造成的延迟。那么, +### a.假设发送一条消息需要耗时$10^{-9}$秒。如果并行程序使用1000个处理器来运行,每个处理器速度和单个处理器运行串行程序速度一样,那么该并行程序的运行需要多少时间? +$10^{12} / 10^6 / 1000 + 10^{-9}*10^9*999 = 1999s$ + +### b.假设发送一条消息需要耗费$10^{-3}$秒,这个使用1000个处理器的并行程序运行需要多少时间? +$10^{12} / 10^6 / 1000 + 10^{-3}*10^9*999 = 999001000s = 11562.5d$ + +## Q5 对各种分布式内存互连网络,写出其包含的链路总数的表达式。 +## 注:只计算网络内部的链路数,即交换器之间的链路。 +设节点个数为p + +| 结构 | 链路总数 | bisection width | +| --- | ------- | --------------- | +| 环 | $p$ | $2$ | +| torodial mesh($p=q^2,q为偶数$) | $2p$ | $2\sqrt{p}(即2q)$ | +| 完全图 | $\frac{p(p-1)}{2}$ | $\frac{p}{2}$ | +| 超立方体($p=2^d$) | $\frac{p}{2}\log_2p$ | $\frac{p}{2}$ | + +| 结构 | 交换机总数 | 链路总数 | +| --- | ------- | ------- | +| 交叉开关矩阵 | $p^2$ | $2p^2(若只考虑矩阵内部,则为2p(p-1))$ | +| omega网络($p=2^d$) | $\frac{p}{2}\log_2p$ | $p(\log_2(p)+1)(若只考虑omega网络内部,则为p(\log_2(p)-1))$ | + +## Q6 为了定义间接网络的等分宽度,我们将处理器分为两组,每组拥有一半数量的处理器。然后,在网络的任意处移除链接,使两组之间不在连接。移除的最小连接数就是该网络的等分宽度。请说明一个8x8的交叉开关矩阵的等分宽度小于或等于8,并说明一个拥有8个处理器的omega网络的等分宽度小于等于4 +1. 8x8的交叉开关矩阵 +- 无论怎样划分,只需要切断每个处理器进入交换网络的总计8条链路,即可让任何两个处理器都不能通信 +2. p=8的omega网络 +- 从上到下,依次轮流把处理器分到第一组和第二组,只需切断第一组进入omega网络的总计4条链路,即可让第一组和第二组不能通信 + +## Q7 +### a.假定一个共享内存系统使用监听缓存一致性协议和写回缓存。并且假设0号核缓存里有变量x,并执行赋值语句x=5。1号核的缓存里没有变量x。当0号核更新x后,1号核开始尝试执行y=x。y被赋的值是多少?为什么? +是x的原始值,因为0号核执行了`x = 5`后没有把它写会内存 + +### b.假定上面的共享内存系统使用的是基于目录的协议,则y的值将是多少?为什么? +结论和原因都同上 + +### c.你能否为前面两部分中发现的问题提出解决方案? +在硬件设计的层面:对于write-back的缓存,需要在后来的核发出读请求时在核间广播,监听到这个广播的核如果发现自己的cache内含有这个读的地址的已经修改过的cache,需要把它写回内存 + +如果硬件这样实现的开销太大,也可以由程序员来手动控制这一过程:通过设置某个标志来表示x已经成功更新,从而允许其它的核来读取它 + +## Q8 假定$T_{串行} = n$,并且$T_{并行} = n / p + \log_2p$。时间单位为毫秒。如果以倍率k增加p,那么为了保持效率值的恒定,需要如何增加n?请给出公式。如果我们将进程数从8加倍到16,则n的增加又是多少?该并行程序是可拓展的吗?该并行程序是弱可扩展的吗? +$$E = \frac{n / p}{n / p + \log_2p} = \frac{n' / kp}{n' / kp + \log_2kp}$$ +$$solved\ to\ n' = n\frac{k\log_2kp}{\log_2p}$$ + +对于从p=8增加到p=16的例子,$\frac{n'}{n} = \frac{2\log_216}{\log_28} = \frac{8}{3}$ + +它是可拓展的,因为$n'$有解 + +它不是弱可拓展的,因为$\frac{n'}{n} > k$ \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/1.c" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/1.c" new file mode 100644 index 000000000..22af9f7ca --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/1.c" @@ -0,0 +1,45 @@ +// Exercise 3.11(d) +#include +#include +#include +#include + +#define LOCAL_N 10 + +int main(int argc, char **argv) { + int pid, nproc; + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &pid); + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + srand(pid + time(NULL)); + int a[LOCAL_N]; + for (int i = 0; i < LOCAL_N; ++i) { + a[i] = rand() % 10; + } + for (int i = 1; i < LOCAL_N; ++i) { + a[i] += a[i - 1]; + } + // Note: MPI_Scan(range, out range) won't work here + // MPI_Scan actually do the op on every scalar in the range and result length won't change + // so MPI_Scan(a, res, LOCAL_N, MPI_INT, MPI_SUM, MPI_COMM_WORLD) + // will result in res = [prefix sum of a[0], ..., prefix sum of a[LOCAL_N - 1], 0 repeat LOCAL_N * (nproc - 1)] + // so for proc 0, res = [self a[0], ... self a[LOCAL_N - 1], 0 repeat LOCAL_N * (nproc - 1)] + int pre; + MPI_Scan(&a[LOCAL_N - 1], &pre, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + pre -= a[LOCAL_N - 1]; + for (int i = 0; i < LOCAL_N; ++i) { + a[i] += pre; + } + if (pid == 0) { + int *res = malloc(nproc * LOCAL_N * sizeof(int)); + MPI_Gather(a, LOCAL_N, MPI_INT, res, LOCAL_N, MPI_INT, 0, MPI_COMM_WORLD); + for (int i = 0; i < nproc * LOCAL_N; ++i) { + printf("%d ", res[i]); + } + puts(""); + free(res); + } else { + MPI_Gather(a, LOCAL_N, MPI_INT, NULL, LOCAL_N, MPI_INT, 0, MPI_COMM_WORLD); + } + MPI_Finalize(); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/2.c" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/2.c" new file mode 100644 index 000000000..23a32fbbf --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/2.c" @@ -0,0 +1,281 @@ +/* + * Purpose: Use MPI to implement a program that creates a histogram + * + * Input: Number of bins + * Minimum measurement + * Maximum measurement + * Number of data items + * + * Output: Histogram of the data. If DEBUG flag is set, the data + * generated. + * + * Notes: + * 1. The number of data must be evenly divisible by the number of processes + * 2. The program generates random floats x in the range + * + * min measurement <= x < max measurement + * + * 3. If i >= 1, the ith bin contains measurements x in the + * range + * + * bin_maxes[i-1] <= x < bin_maxes[i] + * + * Bin 0 will contain measurements x in the range + * + * min measurement <= x < bin_maxes[0] + * + * + * IPP: Programming Assignment 3.1 + */ + +#include +#include +#include + +void Get_input(int* bin_count_p, float* min_meas_p, float* max_meas_p, + int* data_count_p, int* local_data_count_p, int my_rank, + int comm_sz, MPI_Comm comm); +void Check_for_error(int local_ok, char fname[], char message[], + MPI_Comm comm); +void Gen_data(float local_data[], int local_data_count, int data_count, + float min_meas, float max_meas, int my_rank, MPI_Comm comm); +void Set_bins(float bin_maxes[], int loc_bin_cts[], float min_meas, + float max_meas, int bin_count, int my_rank, MPI_Comm comm); +void Find_bins(int bin_counts[], float local_data[], int loc_bin_cts[], + int local_data_count, float bin_maxes[], + int bin_count, float min_meas, MPI_Comm comm); +int Which_bin(float data, float bin_maxes[], int bin_count, + float min_meas); +void Print_histo(float bin_maxes[], int bin_counts[], int bin_count, + float min_meas); + +/*---------------------------------------------------------------------*/ +int main(void) { + int bin_count; + + float min_meas, max_meas; + float* bin_maxes; + int* bin_counts; + int* loc_bin_cts; + + int data_count; + int local_data_count; + + float* data; + float* local_data; + + int my_rank, comm_sz; + MPI_Comm comm; + + MPI_Init(NULL, NULL); + comm = MPI_COMM_WORLD; + MPI_Comm_size(comm, &comm_sz); + MPI_Comm_rank(comm, &my_rank); + + // get user inputs for bin_count, max_meas, min_meas, and data_count + Get_input(&bin_count, &min_meas, &max_meas, &data_count, + &local_data_count, my_rank, comm_sz, comm); + + // allocate arrays + bin_maxes = malloc(bin_count*sizeof(float)); + bin_counts = malloc(bin_count*sizeof(int)); + loc_bin_cts = malloc(bin_count*sizeof(int)); + data = malloc(data_count*sizeof(float)); + local_data = malloc(local_data_count*sizeof(float)); + + Set_bins(bin_maxes, loc_bin_cts, min_meas, max_meas, bin_count, + my_rank, comm); + Gen_data(local_data, local_data_count, data_count, min_meas, + max_meas, my_rank, comm); + Find_bins(bin_counts, local_data, loc_bin_cts, local_data_count, + bin_maxes, bin_count, min_meas, comm); + + if (my_rank == 0) + Print_histo(bin_maxes, bin_counts, bin_count, min_meas); + + free(bin_maxes); + free(bin_counts); + free(loc_bin_cts); + free(data); + free(local_data); + MPI_Finalize(); + return 0; + +} /* main */ + + +/*---------------------------------------------------------------------*/ +void Print_histo( + float bin_maxes[] /* in */, + int bin_counts[] /* in */, + int bin_count /* in */, + float min_meas /* in */) { + int i, j; + float bin_max, bin_min; + + for (i = 0; i < bin_count; i++) { + bin_max = bin_maxes[i]; + bin_min = (i == 0) ? min_meas: bin_maxes[i-1]; + printf("%.3f-%.3f:\t", bin_min, bin_max); + for (j = 0; j < bin_counts[i]; j++) + printf("X"); + printf("\n"); + } +} /* Print_histo */ + + +/*---------------------------------------------------------------------*/ +void Find_bins( + int bin_counts[] /* out */, + float local_data[] /* in */, + int loc_bin_cts[] /* out */, + int local_data_count /* in */, + float bin_maxes[] /* in */, + int bin_count /* in */, + float min_meas /* in */, + MPI_Comm comm){ + /* Use a for loop to find bins, the statement in the loop can be: + bin = Which_bin(local_data[i], bin_maxes, bin_count, min_meas); + Then, calculate the global sum using collective communication. + */ + for (int i = 0; i < local_data_count; ++i) { + ++loc_bin_cts[Which_bin(local_data[i], bin_maxes, bin_count, min_meas)]; + } + MPI_Reduce(loc_bin_cts, bin_counts, bin_count, MPI_INT, MPI_SUM, 0, comm); +} /* Find_bins */ + + +/*---------------------------------------------------------------------*/ +int Which_bin(float data, float bin_maxes[], int bin_count, + float min_meas) { + if (min_meas <= data && data < bin_maxes[0]) { + return 0; + } else { + for (int i = 1; i < bin_count; ++i) { + if (bin_maxes[i - 1] <= data && data < bin_maxes[i]) { + return i; + } + } + return bin_count - 1; + } +} /* Which_bin */ + + +/*---------------------------------------------------------------------*/ +void Set_bins( + float bin_maxes[] /* out */, + int loc_bin_cts[] /* out */, + float min_meas /* in */, + float max_meas /* in */, + int bin_count /* in */, + int my_rank /* in */, + MPI_Comm comm /* in */) { + + if (my_rank == 0) { + int i; + float bin_width; + bin_width = (max_meas - min_meas) / bin_count; + + for (i = 0; i < bin_count; i++) { + loc_bin_cts[i] = 0; + bin_maxes[i] = min_meas + (i+1)*bin_width; + } + + } + + // set bin_maxes for each proc + MPI_Bcast(bin_maxes, bin_count, MPI_FLOAT, 0, comm); + + // reset loc_bin_cts of each proc + MPI_Bcast(loc_bin_cts, bin_count, MPI_INT, 0, comm); +} /* Set_bins */ + + +/*---------------------------------------------------------------------*/ +void Gen_data( + float local_data[] /* out */, + int local_data_count /* in */, + int data_count /* in */, + float min_meas /* in */, + float max_meas /* in */, + int my_rank /* in */, + MPI_Comm comm /* in */) { + int i; + float* data = NULL; + + if (my_rank ==0) { + data = malloc(data_count*sizeof(float)); + srandom(1); + for (i = 0; i < data_count; i++) + data[i] = + min_meas + (max_meas - min_meas)*random()/((double) RAND_MAX); +# ifdef DEBUG + printf("Generated data:\n "); + for (i = 0; i < data_count; i++) + printf("%.3f ", data[i]); + printf("\n\n"); +# endif + MPI_Scatter(data, local_data_count, MPI_FLOAT, local_data, + local_data_count, MPI_FLOAT, 0, comm); + free(data); + } else { + MPI_Scatter(data, local_data_count, MPI_FLOAT, local_data, + local_data_count, MPI_FLOAT, 0, comm); + } + +} /* Gen_data */ + +/*---------------------------------------------------------------------*/ +void Get_input(int* bin_count_p, float* min_meas_p, float* max_meas_p, + int* data_count_p, int* local_data_count_p, int my_rank, + int comm_sz, MPI_Comm comm) { + + int local_ok = 1; + + if (my_rank == 0) { + printf("Enter the number of bins\n"); + scanf("%d", bin_count_p); + printf("Enter the minimum measurement\n"); + scanf("%f", min_meas_p); + printf("Enter the maximum measurement\n"); + scanf("%f", max_meas_p); + printf("Enter the number of data\n"); + scanf("%d", data_count_p); + } + + MPI_Bcast(bin_count_p, 1, MPI_INT, 0, comm); + MPI_Bcast(min_meas_p, 1, MPI_INT, 0, comm); + MPI_Bcast(max_meas_p, 1, MPI_INT, 0, comm); + MPI_Bcast(data_count_p, 1, MPI_INT, 0, comm); + + if(*data_count_p % comm_sz != 0) local_ok = 0; + + Check_for_error(local_ok, "Get_input", + "data_count must be evenly divisible by comm_sz", + comm); + *local_data_count_p = *data_count_p / comm_sz; + +} /* Get_input */ + + +/*---------------------------------------------------------------------*/ +void Check_for_error( + int local_ok /* in */, + char fname[] /* in */, + char message[] /* in */, + MPI_Comm comm /* in */) { + int ok; + + MPI_Allreduce(&local_ok, &ok, 1, MPI_INT, MPI_MIN, comm); + if (ok == 0) { + int my_rank; + MPI_Comm_rank(comm, &my_rank); + if (my_rank == 0) { + fprintf(stderr, "Proc %d > In %s, %s\n", my_rank, fname, + message); + fflush(stderr); + } + MPI_Finalize(); + exit(-1); + } +} /* Check_for_error */ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/3.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/3.cpp" new file mode 100644 index 000000000..380f4e23c --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/3.cpp" @@ -0,0 +1,124 @@ +#include +#include +#include +#include +#include + +std::pair, std::unique_ptr> gen(int n) { + std::unique_ptr mat(new double[n * n]); + std::unique_ptr vec(new double[n]); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + mat[i * n + j] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + } + for (int i = 0; i < n; ++i) { + vec[i] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + return {std::move(mat), std::move(vec)}; +} + +std::unique_ptr serial(const std::unique_ptr &mat, const std::unique_ptr &vec, int n) { + std::unique_ptr ret(new double[n]); + for (int i = 0; i < n; ++i) { + double tmp = 0.0; + for (int j = 0; j < n; ++j) { + tmp += mat[i * n + j] * vec[j]; + } + ret[i] = tmp; + } + return ret; +} + +void run(int n) { + int pid, nproc; + double beg[2], elapsed[2]; + + MPI_Init(nullptr, nullptr); + MPI_Comm_rank(MPI_COMM_WORLD, &pid); + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + + if (n % nproc != 0) { + printf("n must be divisible by nproc, but n = %d, nproc= %d\n", n, nproc); + exit(-1); + } + + int each = n / nproc; + std::unique_ptr loc_mat(new double[n * each]), loc_vec(new double[each]), ser_res; + + MPI_Datatype vec_t, col_t; + // vec_t is a type like struct { double[each]@0, double[each]@(n * sizeof(double)), ..., double[each]@((n * n - 1) * sizeof(double))} + MPI_Type_vector(n, each, n, MPI_DOUBLE, &vec_t); + // col_t is a type like vec_t, but it behaves like sizeof(col_t) == each * sizeof(double) in MPI_Scatter + MPI_Type_create_resized(vec_t, 0, each * sizeof(double), &col_t); + MPI_Type_commit(&col_t); + std::unique_ptr size(new int[nproc]), off(new int[nproc]); + for (int i = 0; i < nproc; ++i) { + size[i] = 1; + off[i] = i; + } + + decltype(gen(n)) mat_vec(nullptr, nullptr); + if (pid == 0) { + mat_vec = gen(n); + ser_res = serial(mat_vec.first, mat_vec.second, n); + } + + MPI_Barrier(MPI_COMM_WORLD); + beg[0] = MPI_Wtime(); + + MPI_Scatterv(mat_vec.first.get(), size.get(), off.get(), col_t, loc_mat.get(), n * each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Scatter(mat_vec.second.get(), each, MPI_DOUBLE, loc_vec.get(), each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + MPI_Barrier(MPI_COMM_WORLD); + beg[1] = MPI_Wtime(); + + std::unique_ptr y(new double[n]); + for (int i = 0; i < n; ++i) { + double tmp = 0.0; + for (int j = 0; j < each; ++j) { + tmp += loc_mat[i * each + j] * loc_vec[j]; + } + y[i] = tmp; + } + + if (pid == 0) { + std::unique_ptr res(new double[n]); + MPI_Reduce(y.get(), res.get(), n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + + double now = MPI_Wtime(); + beg[0] = now - beg[0], beg[1] = now - beg[1]; + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + printf("elapsed(+MPI_Scatter): %.4lfs\n", elapsed[0]); + printf("elapsed(-MPI_Scatter): %.4lfs\n", elapsed[1]); + + double l1 = 0.0, l2 = 0.0, linf = 0.0; + for (int i = 0; i < n; ++i) { + double diff = fabs(ser_res[i] - res[i]); + l1 += diff; + l2 += diff * diff; + linf = fmax(linf, diff); + } + l2 = sqrtf(l2); + printf("error(1 norm): %.20lf\n", l1); + printf("error(2 norm): %.20lf\n", l2); + printf("error(infinite norm): %.20lf\n", linf); + } else { + MPI_Reduce(y.get(), nullptr, n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + + double now = MPI_Wtime(); + beg[0] -= now, beg[1] -= now; + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + } + MPI_Type_free(&col_t); + MPI_Finalize(); +} + +int main(int argc, char *argv[]) { + if (argc < 2) { + printf("Usage: %s n\n", argv[0]); + return -1; + } + int n = atoi(argv[1]); + run(n); +} diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/4.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/4.cpp" new file mode 100644 index 000000000..7fd64a691 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/4.cpp" @@ -0,0 +1,135 @@ +#include +#include +#include +#include +#include + +std::pair, std::unique_ptr> gen(int n) { + std::unique_ptr mat(new double[n * n]); + std::unique_ptr vec(new double[n]); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + mat[i * n + j] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + } + for (int i = 0; i < n; ++i) { + vec[i] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + return {std::move(mat), std::move(vec)}; +} + +std::unique_ptr serial(const std::unique_ptr &mat, const std::unique_ptr &vec, int n) { + std::unique_ptr ret(new double[n]); + for (int i = 0; i < n; ++i) { + double tmp = 0.0; + for (int j = 0; j < n; ++j) { + tmp += mat[i * n + j] * vec[j]; + } + ret[i] = tmp; + } + return ret; +} + +void run(int n) { + int pid, nproc; + double beg[2], elapsed[2]; + + MPI_Init(nullptr, nullptr); + MPI_Comm_rank(MPI_COMM_WORLD, &pid); + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + + int sq_nproc = sqrtf(nproc); + if (sq_nproc * sq_nproc != nproc) { + printf("nproc must be complete square number, but nproc= %d\n", nproc); + exit(-1); + } + if (n % sq_nproc != 0) { + printf("n must be divisible by sqrt(nproc), but n = %d, nproc= %d\n", n, nproc); + exit(-1); + } + + int each = n / sq_nproc; + std::unique_ptr loc_mat(new double[each * each]), vec(nullptr), ser_res; + + MPI_Datatype vec_t, block_t; + MPI_Type_vector(each, each, n, MPI_DOUBLE, &vec_t); + MPI_Type_create_resized(vec_t, 0, each * sizeof(double), &block_t); + MPI_Type_commit(&block_t); + std::unique_ptr size(new int[nproc]), off(new int[nproc]); + for (int i = 0; i < sq_nproc; ++i) { + for (int j = 0; j < sq_nproc; ++j) { + size[i * sq_nproc + j] = 1; + off[i * sq_nproc + j] = i * n + j; + } + } + + std::unique_ptr mat(nullptr); + if (pid == 0) { + auto mat_vec = gen(n); + ser_res = serial(mat_vec.first, mat_vec.second, n); + mat = std::move(mat_vec.first); + vec = std::move(mat_vec.second); + } else { + vec.reset(new double[n]); + } + + MPI_Barrier(MPI_COMM_WORLD); + beg[0] = MPI_Wtime(); + + MPI_Scatterv(mat.get(), size.get(), off.get(), block_t, loc_mat.get(), each * each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(vec.get(), n, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + MPI_Barrier(MPI_COMM_WORLD); + beg[1] = MPI_Wtime(); + + std::unique_ptr y(new double[n]{}); + int vec_off = pid % sq_nproc * each, y_off = pid / sq_nproc * each; + for (int i = 0; i < each; ++i) { + double tmp = 0.0; + for (int j = 0; j < each; ++j) { + tmp += loc_mat[i * each + j] * vec[vec_off + j]; + } + y[y_off + i] = tmp; + } + + if (pid == 0) { + std::unique_ptr res(new double[n]); + MPI_Reduce(y.get(), res.get(), n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + + double now = MPI_Wtime(); + beg[0] = now - beg[0], beg[1] = now - beg[1]; + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + printf("elapsed(+MPI_Scatter): %.4lfs\n", elapsed[0]); + printf("elapsed(-MPI_Scatter): %.4lfs\n", elapsed[1]); + + double l1 = 0.0, l2 = 0.0, linf = 0.0; + for (int i = 0; i < n; ++i) { + double diff = fabs(ser_res[i] - res[i]); + l1 += diff; + l2 += diff * diff; + linf = fmax(linf, diff); + } + l2 = sqrtf(l2); + printf("error(1 norm): %.20lf\n", l1); + printf("error(2 norm): %.20lf\n", l2); + printf("error(infinite norm): %.20lf\n", linf); + } else { + MPI_Reduce(y.get(), nullptr, n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + + double now = MPI_Wtime(); + beg[0] -= now, beg[1] -= now; + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + } + + MPI_Type_free(&block_t); + MPI_Finalize(); +} + +int main(int argc, char *argv[]) { + if (argc < 2) { + printf("Usage: %s n\n", argv[0]); + return -1; + } + int n = atoi(argv[1]); + run(n); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/makefile" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/makefile" new file mode 100644 index 000000000..16a9f1db5 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/makefile" @@ -0,0 +1,18 @@ +.PHONY: all clean + +all: 1 2 3 4 + +1: 1.c + mpicc -O3 -std=c11 1.c -o 1 + +2: 2.c + mpicc -O3 -std=c11 2.c -o 2 + +3: 3.cpp + mpiicc -O3 -std=c++11 3.cpp -o 3 + +4: 4.cpp + mpiicc -O3 -std=c++11 4.cpp -o 4 + +clean: + rm 1 2 3 4 \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/pic1.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/pic1.png" new file mode 100644 index 000000000..25d4058b4 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/pic1.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/pic2.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/pic2.png" new file mode 100644 index 000000000..0ee588c1e Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/pic2.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/report.md" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/report.md" new file mode 100644 index 000000000..7f78bc908 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/report.md" @@ -0,0 +1,300 @@ +--- +title: "第三次作业报告" +author: "李晨昊 2017011466" +date: "2019-4-9" +output: + pdf_document: + latex_engine: xelatex + number_sections: yes + toc: yes + word_document: + toc: yes + html_document: + toc: yes +header-includes: \usepackage{ctex} \usepackage{graphicx} +--- + +\newpage + +# 运行方式 +直接执行`make`即可分别编译`1.c`,`2.c`,`3.cpp`,`4.cpp`得到四个可执行文件。其中1与2均可直接用`srun -n [n proc] -l [1|2]`运行,3和4需要提供为程序提供一个命令行参数表示矩阵大小`srun -n [n proc] -l [3|4] [matrix size]`,程序内部会对这个参数是否合法做出检查。 + +# Q1 +## 解题思路 +本题要求每个进程生成有10个元素的随机数组,按顺序拼接成10p个元素的全局数组,计算这个全局数组的10p个前缀和, 存在p个进程中。最后,将前缀和传输给0号进程输出。 + +使用`MPI_Scan`来从各个进程收集局部的前缀和信息。先在每个进程中计算出局部的前缀和,再调用`MPI_Scan`计算从开头到本进程的数字之和。此值减去本进程的所有数字之和,即是本进程的计算出的前缀和距离全局前缀和的差距,将本进程的计算结果都加上这个差距后,汇总到0号进程输出即可。 + +## 关键程序 +```c + // 计算局部前缀和 + for (int i = 1; i < LOCAL_N; ++i) { + a[i] += a[i - 1]; + } + int pre; + // 计算从开头到本进程的数字之和 + MPI_Scan(&a[LOCAL_N - 1], &pre, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + // pre - a[LOCAL_N - 1]即是本进程之前的所有数字之和 + pre -= a[LOCAL_N - 1]; + // 加上这个值后,本进程计算出的值变成全局前缀和 + for (int i = 0; i < LOCAL_N; ++i) { + a[i] += pre; + } +``` + +## 测试结果 +由于本题的设定,数据的规模本身是正比于进程数目的,因此无法对于一个固定的数据规模来计算不同进程数下的加速比。由于加速比(一般而言)总是小于核数,因此本题在较多进程上的运行速度一定慢于较少进程,只用考察速度之比与1差距多少即可。 + +为了方便时间测定,去掉了程序中的输出,并且让每个进程生成$10^7$个随机数而非10个。注意表中的`与p=1的速度比`并不是平时讨论的加速比。 + +时间消耗表格如下 + +| p | 耗时 | 与p=1的时间比 | 与p=1的速度比 | +| -- | --- | ------------ | ----------- | +| 1 | 0.212s | 1.000 | 1.000 | +| 2 | 0.337s | 1.590 | 1.258 | +| 3 | 0.396s | 1.868 | 1.606 | +| 4 | 0.424s | 2.000 | 2.000 | +| 5 | 0.395s | 1.863 | 2.684 | +| 6 | 0.416s | 1.962 | 3.058 | +| 7 | 0.438s | 2.066 | 3.388 | +| 8 | 0.452s | 2.132 | 3.752 | +| 9 | 0.472s | 2.226 | 4.042 | +| 10 | 0.490s | 2.311 | 4.327 | +| 11 | 0.516s | 2.434 | 4.519 | +| 12 | 0.536s | 2.528 | 4.746 | +| 13 | 0.563s | 2.656 | 4.895 | +| 14 | 0.591s | 2.788 | 5.022 | +| 15 | 0.604s | 2.849 | 5.265 | +| 16 | 0.627s | 2.958 | 5.410 | +| 17 | 0.646s | 3.047 | 5.579 | +| 18 | 0.668s | 3.151 | 5.713 | +| 19 | 0.694s | 3.274 | 5.804 | +| 20 | 0.711s | 3.354 | 5.963 | +| 21 | 0.740s | 3.491 | 6.016 | +| 22 | 0.752s | 3.547 | 6.202 | +| 23 | 0.772s | 3.642 | 6.316 | +| 24 | 0.812s | 3.830 | 6.266 | + +# Q2 +## 解题思路 +本题要求完成`Find_bins`函数和`Which_bin`函数。 + +`Which_bin`的作用为在一个已排序的数组中寻找某个位置,使得某值位于两个数组元素之间。其实本来在本题的设定中,可以直接使用一步除法就得到这个结果的,但是既然要求这样做,我也没太大必要来优化,直接采用朴素的线性查找即可。 + +`Find_bins`的作用为为每个进程局部的桶统计计数,这一步只需要调用`Which_bin`即可。完成之后,使用`MPI_Reduce`将所有进程的统计结果求和。 + +## 关键程序 +```c +int Which_bin(float data, float bin_maxes[], int bin_count, + float min_meas) { + if (min_meas <= data && data < bin_maxes[0]) { + return 0; // 特判data位于第一个桶内 + } else { + // 线性查找 + for (int i = 1; i < bin_count; ++i) { + if (bin_maxes[i - 1] <= data && data < bin_maxes[i]) { + return i; + } + } + return bin_count - 1; // 特判data比所有元素都大(其实这不会发生) + } +} + +void Find_bins( + int bin_counts[] /* out */, + float local_data[] /* in */, + int loc_bin_cts[] /* out */, + int local_data_count /* in */, + float bin_maxes[] /* in */, + int bin_count /* in */, + float min_meas /* in */, + MPI_Comm comm){ + // 统计每个桶的计数 + for (int i = 0; i < local_data_count; ++i) { + ++loc_bin_cts[Which_bin(local_data[i], bin_maxes, bin_count, min_meas)]; + } + // 全局求和 + MPI_Reduce(loc_bin_cts, bin_counts, bin_count, MPI_INT, MPI_SUM, 0, comm); +} +``` + +## 测试结果 +运行`srun -n 10 ./2`,输入`10 0 100 400`,得到 +``` +0.000-10.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +10.000-20.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +20.000-30.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +30.000-40.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +40.000-50.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +50.000-60.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +60.000-70.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +70.000-80.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +80.000-90.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +90.000-100.000: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +``` +可见结果比较均匀,是合理的。输出的结果与指定的进程数目无关,基本证明实现是正确的。 + +# Q3 +## 解题思路 +本题要求将矩阵按列分配到不同进程,这就导致每个进程分配到的内存并不是连续的。经测试,循环分配每行元素虽然可行,但是效率较低,因此考虑采用MPI的类型定义接口,实现在一次`MPI_Scatter`中将矩阵的所有元素分配出去。 + +这里使用`MPI_Type_vector`和`MPI_Type_create_resized`来得到所需的类型。达到的效果是,定义这样一种类型:它由`n`块连续的`n / nproc`个`double`构成,相应`double`块间的间隔为`n * sizeof(double)`,利用这种类型的指针遍历内存的步长为`n / nproc * sizeof(double)`。这样以来,在原始的矩阵上(row-major连续存储),利用这种类型步进`nproc`次就遍历了所有的元素,这样的分配是符合要求的。 + +我认为到此为止,应该是可以使用`MPI_Scatter`来分配矩阵元素的。然而因为一些目前无法解释的原因,这会在某些数据下发生segment fault,因此只能改用`MPI_Scatterv`来实现。 + +按列分配了矩阵元素后,再均分向量元素,这样每个进程就得到了自己计算需要的所有数据。它的计算结果是包括最终结果的每个分量,但都只是部分和,因此最后对中间结果进行`MPI_Reduce`即可。 + +## 关键程序 +类型定义部分 +```c + MPI_Datatype vec_t, col_t; + MPI_Type_vector(n, each, n, MPI_DOUBLE, &vec_t); + MPI_Type_create_resized(vec_t, 0, each * sizeof(double), &col_t); + MPI_Type_commit(&col_t); + std::unique_ptr size(new int[nproc]), off(new int[nproc]); + for (int i = 0; i < nproc; ++i) { + size[i] = 1; + off[i] = i; + } +``` +数据分发部分 +```c + // 按列分配数组元素 + MPI_Scatterv(mat_vec.first.get(), size.get(), off.get(), + col_t, loc_mat.get(), n * each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + // 分块分配向量元素 + MPI_Scatter(mat_vec.second.get(), each, MPI_DOUBLE, + loc_vec.get(), each, MPI_DOUBLE, 0, MPI_COMM_WORLD); +``` +计算部分 +```c + std::unique_ptr y(new double[n]); + for (int i = 0; i < n; ++i) { + double tmp = 0.0; + for (int j = 0; j < each; ++j) { + tmp += loc_mat[i * each + j] * loc_vec[j]; + } + y[i] = tmp; + } +``` + +## 测试结果 +考虑数据分发,时间/加速比表格如下 + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +| -- | ------ | ------ | ------- | ------- | +| 1 | 0.0395s/1.000 | 0.1563s/1.000 | 0.6401s/1.000 | 2.6686s/1.000 | +| 2 | 0.0423s/1.071 | 0.1506s/1.038 | 0.5847s/1.095 | 2.4319s/1.097 | +| 3 | 0.0306s/1.291 | 0.1061s/1.473 | 0.4078s/1.570 | 1.5989s/1.669 | +| 4 | 0.0254s/1.555 | 0.0930s/1.681 | 0.3158s/2.027 | 1.2397s/2.153 | +| 5 | 0.0241s/1.639 | 0.0864s/1.809 | 0.2642s/2.422 | 1.0355s/2.577 | +| 6 | 0.0241s/1.639 | 0.0862s/1.813 | 0.2304s/2.778 | 0.8889s/3.002 | +| 8 | 0.0234s/1.688 | 0.0833s/1.876 | 0.3149s/2.032 | 0.7288s/3.662 | +| 9 | 0.0230s/1.717 | 0.0825s/1.894 | 0.3196s/2.002 | 0.6896s/3.870 | +| 10 | 0.0229s/1.724 | 0.0830s/1.883 | 0.3156s/2.028 | 0.6728s/3.966 | +| 12 | 0.0230s/1.717 | 0.0831s/1.880 | 0.3142s/2.037 | 0.6780s/3.936 | +| 15 | 0.0242s/1.632 | 0.0837s/1.867 | 0.3298s/1.940 | 1.2153s/2.196 | +| 16 | 0.0246s/1.605 | 0.0890s/1.756 | 0.3295s/1.942 | 1.2026s/2.219 | +| 18 | 0.0257s/1.536 | 0.0867s/1.802 | 0.3291s/1.945 | 1.2537s/2.129 | +| 20 | 0.0260s/1.519 | 0.0886s/1.764 | 0.3264s/1.961 | 1.2059s/2.213 | +| 24 | 0.0262s/1.507 | 0.0884s/1.768 | 0.3174s/2.016 | 1.2251s/2.178 | + +不考虑数据分发,时间/加速比表格如下 + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +| -- | ------ | ------ | ------- | ------- | +| 1 | 0.0084s/1.000 | 0.0328s/1.000 | 0.1368s/1.000 | 0.6466s/1.000 | +| 2 | 0.0043s/1.953 | 0.0171s/1.918 | 0.0675s/2.027 | 0.2852s/2.267 | +| 3 | 0.0031s/2.710 | 0.0119s/2.756 | 0.0472s/2.898 | 0.1905s/3.394 | +| 4 | 0.0026s/3.231 | 0.0092s/3.565 | 0.0369s/3.707 | 0.1479s/4.372 | +| 5 | 0.0023s/3.652 | 0.0081s/4.049 | 0.0300s/4.560 | 0.1189s/5.438 | +| 6 | 0.0023s/3.652 | 0.0079s/4.152 | 0.0264s/5.182 | 0.1141s/5.667 | +| 8 | 0.0022s/3.818 | 0.0072s/4.556 | 0.0219s/6.247 | 0.1072s/6.032 | +| 9 | 0.0022s/3.818 | 0.0072s/4.556 | 0.0246s/5.561 | 0.1067s/6.060 | +| 10 | 0.0025s/3.360 | 0.0075s/4.373 | 0.0249s/5.494 | 0.1066s/6.066 | +| 12 | 0.0023s/3.652 | 0.0073s/4.493 | 0.0209s/6.545 | 0.0980s/6.598 | +| 15 | 0.0022s/3.818 | 0.0061s/5.377 | 0.0207s/6.609 | 0.0854s/7.571 | +| 16 | 0.0019s/4.421 | 0.0058s/5.655 | 0.0207s/6.609 | 0.0738s/8.762 | +| 18 | 0.0022s/3.818 | 0.0056s/5.857 | 0.0186s/7.355 | 0.0659s/9.812 | +| 20 | 0.0019s/4.421 | 0.0052s/6.308 | 0.0170s/8.047 | 0.0597s/10.831 | +| 24 | 0.0018s/4.667 | 0.0044s/7.455 | 0.0142s/9.634 | 0.0544s/11.886 | + +加速比折线图如下。可以发现并行效率在n较大而p较小的时候可以略微超过1,我认为这与cache有关。举例来说,假设某个任务需要对一个数组迭代多次,而这个数组无法完整装入per core cache(L1/L2),那么每次迭代的时候都会发生cache缺失(也不一定需要访存,可能在共享的L3 cache处解决);但是如果将任务分配到两个核上,这个数组的两半可能可以分别装入per core cache,这样除了第一次迭代,后面的迭代都不会发生cache缺失,这样每个核的耗时都可能会小于原本耗时的一半,并行效率就超过了1。 + +\begin{center} +\includegraphics[width=0.9\linewidth]{pic1.png} +\end{center} + +# Q4 +## 解题思路 +与Q3类似,本题需要自定义数据类型才能一次把矩阵分配到各个子进程中。除此之外,观察得每个子矩阵在原始矩阵中的偏移量的差值并不是常数,因此需要使用`MPI_Scatterv`来指定子块的大小和偏移量。 + +对于向量的分配,这里简单的把它广播到每个进程中,让它们自己选择自己需要的部分。 + +## 关键程序 +类型定义部分 +```c + MPI_Datatype vec_t, block_t; + MPI_Type_vector(each, each, n, MPI_DOUBLE, &vec_t); + MPI_Type_create_resized(vec_t, 0, each * sizeof(double), &block_t); + MPI_Type_commit(&block_t); + std::unique_ptr size(new int[nproc]), off(new int[nproc]); + for (int i = 0; i < sq_nproc; ++i) { + for (int j = 0; j < sq_nproc; ++j) { + // 每个进程接收到一个块 + size[i * sq_nproc + j] = 1; + // 这里的偏移量是和MPI_Type_create_resized中制定的大小一起考虑的 + // 换算成字节的偏移量,是(i * n + j) * each * sizeof(double) + off[i * sq_nproc + j] = i * n + j; + } + } +``` +数据分发部分 +```c + // 分块分配矩阵 + MPI_Scatterv(mat.get(), size.get(), off.get(), block_t, loc_mat.get(), + each * each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + // 广播向量 + MPI_Bcast(vec.get(), n, MPI_DOUBLE, 0, MPI_COMM_WORLD); +``` +计算部分 +```c + // 注意后面的花括号,这里就把y清零了 + std::unique_ptr y(new double[n]{}); + // 自身只使用vec的一部分,只填y的一部分,根据进程序号计算具体是哪一段 + int vec_off = pid % sq_nproc * each, y_off = pid / sq_nproc * each; + for (int i = 0; i < each; ++i) { + double tmp = 0.0; + for (int j = 0; j < each; ++j) { + tmp += loc_mat[i * each + j] * vec[vec_off + j]; + } + y[y_off + i] = tmp; + } +``` + +## 测试结果 +考虑数据分发,时间/加速比表格如下 + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +| -- | ------ | ------ | ------- | ------- | +| 1 | 0.0394s/1.000 | 0.1591s/1.000 | 0.6341s/1.000 | 2.6252s/1.000 | +| 4 | 0.0262s/1.504 | 0.0804s/1.979 | 0.3091s/2.051 | 1.2312s/2.132 | +| 9 | 0.0223s/1.767 | 0.0442s/3.600 | 0.1631s/3.888 | 0.6327s/4.149 | +| 16 | 0.0234s/1.684 | 0.0842s/1.890 | 0.1400s/4.529 | 0.5583s/4.702 | + + +不考虑数据分发,时间/加速比表格如下 + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +| -- | ------ | ------ | ------- | ------- | +| 1 | 0.0084s/1.000 | 0.0333s/1.000 | 0.1404s/1.000 | 0.6461s/1.000 | +| 4 | 0.0026s/3.231 | 0.0094s/3.543 | 0.0371s/3.784 | 0.1528s/4.228 | +| 9 | 0.0026s/3.231 | 0.0069s/4.826 | 0.0272s/5.162 | 0.0968s/6.675 | +| 16 | 0.0019s/4.421 | 0.0056s/5.946 | 0.0207s/6.783 | 0.0748s/8.638 | + +加速比折线图如下 + +\begin{center} +\includegraphics[width=0.9\linewidth]{pic2.png} +\end{center} diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/report.pdf" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/report.pdf" new file mode 100644 index 000000000..4d7d8d0c9 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/report.pdf" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/run.py" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/run.py" new file mode 100644 index 000000000..cbf43527e --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw3/run.py" @@ -0,0 +1,16 @@ +import os + +for n in (3600, 7200, 14400, 28800): + # for i in (1, 2, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24): + for i in (3, 4): + cmd = 'srun -n %d ./3 %d' % (i, n) + print cmd + for _ in range(5): + os.system(cmd) + +# for n in (3600, 7200, 14400, 28800): +# for i in (1, 4, 9, 16): +# cmd = 'srun -n %d ./4 %d' % (i, n) +# print cmd +# for _ in range(5): +# os.system(cmd) diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_1.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_1.cpp" new file mode 100644 index 000000000..04c2971bc --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_1.cpp" @@ -0,0 +1,33 @@ +#include "util.hpp" + +const u32 MAX_MSG = 100; +Mutex mu; +bool ok; +i8 msg[MAX_MSG]; + +void* cons(void*) { + while (true) { + Locker _lk(&mu); + if (ok) { + printf("get msg %s\n", msg); + break; + } + } + return nullptr; +} + +void* prod(void*) { + Locker _lk(&mu); + strcpy(msg, "hello world!"); + ok = true; + return nullptr; +} + +i32 main() { + pthread_t cons, prod; + pthread_create(&cons, nullptr, ::cons, nullptr); + pthread_create(&prod, nullptr, ::prod, nullptr); + + pthread_join(cons, nullptr); + pthread_join(prod, nullptr); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_2.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_2.cpp" new file mode 100644 index 000000000..75466bf43 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_2.cpp" @@ -0,0 +1,50 @@ +#include "util.hpp" + +const u32 MAX_MSG = 100; +Mutex mu; +bool ok; +i8 msg[MAX_MSG]; + +void *thread_fn(void *_tid) { + usize tid = (usize)_tid; + if (tid % 2 == 0) { // prod + while (true) { + Locker _lk(&mu); + if (!ok) { + sprintf(msg, "hello world from tid = %lu", tid); + ok = true; + break; + } + } + } else { // cons + while (true) { + Locker _lk(&mu); + if (ok) { + printf("%lu: get msg: %s\n", tid, msg); + ok = false; + break; + } + } + } + return nullptr; +} + +i32 main(i32 argc, i8 **argv) { + u32 n; + if (argc < 2 || (n = atoi(argv[1])) == 0 || n % 2 == 1) { + printf("usage %s [thread count](positive even)\n", argv[0]); + exit(-1); + } + + Mutex mu; + std::unique_ptr msg; + std::unique_ptr ths(new pthread_t[n]); + + for (usize i = 0; i < n; ++i) { + pthread_create(&ths[i], nullptr, thread_fn, (void *)i); + } + + for (u32 i = 0; i < n; ++i) { + pthread_join(ths[i], nullptr); + } +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_3.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_3.cpp" new file mode 100644 index 000000000..29498dd0b --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/4_7_3.cpp" @@ -0,0 +1,42 @@ +#include "util.hpp" + +const u32 MAX_MSG = 100; +Mutex mu; +bool ok; +i8 msg[MAX_MSG]; +u32 recv, n; + +void *thread_fn(void *_tid) { + usize tid = (usize)_tid; + bool sended = false, recved = false; + while (!sended || !recved) { + Locker _lk(&mu); + if (ok) { + if (!recved && recv == tid) { + printf("%lu: get msg: %s\n", tid, msg); + ok = false; + recved = true; + } + } else if (!sended) { + sprintf(msg, "hello world from tid = %lu", tid); + ok = true; + sended = true; + recv = (tid + 1) % n; + } + } + return nullptr; +} + +i32 main(i32 argc, i8 **argv) { + if (argc < 2 || (n = atoi(argv[1])) == 0) { + printf("usage %s [thread count](positive)\n", argv[0]); + exit(-1); + } + std::unique_ptr ths(new pthread_t[n]); + for (usize i = 0; i < n; ++i) { + pthread_create(&ths[i], nullptr, thread_fn, (void *)i); + } + for (u32 i = 0; i < n; ++i) { + pthread_join(ths[i], nullptr); + } +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bad_volatile.c" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bad_volatile.c" new file mode 100644 index 000000000..05645b618 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bad_volatile.c" @@ -0,0 +1,111 @@ +#include +#include +#include +#include + +#define ALL (10000000) +#define THREAD_NUM (4) +#define EACH (ALL / THREAD_NUM) + +volatile intptr_t sum; +volatile intptr_t flag; + +void *inc(void *_tid) { + intptr_t tid = (intptr_t)_tid; + for (int i = 0; i < EACH; ++i) { + while (flag != tid) + ; + ++sum; + flag = (flag + 1) % THREAD_NUM; + } + return NULL; +} + +int main() { + pthread_t th[THREAD_NUM]; + for (intptr_t i = 0; i < THREAD_NUM; ++i) { + pthread_create(&th[i], NULL, inc, (void *)i); + } + for (intptr_t i = 0; i < THREAD_NUM; ++i) { + pthread_join(th[i], NULL); + } + printf("%ld %d\n", sum, ALL); +} + +// void *hello(void *_) { +// printf("hello world from %ld\n", (long)(_)); +// return NULL; +// } + +// int main() { +// pthread_t ths[10]; +// for (long i = 0; i < 10; ++i) { +// pthread_create(&ths[i], NULL, hello, (void *)i); +// } +// for (long i = 0; i < 10; ++i) { +// pthread_join(ths[i], NULL); +// } +// } + +// intptr_t val1, val2; + +// void *f1(void *_) { +// val1 = 1; +// return (void *)val2; +// } + +// void *f2(void *_) { +// val2 = 1; +// return (void *)val1; +// } + +// int main() { +// int fail = 0; +// for (int i = 0; i < 1000000; ++i) { +// val1 = val2 = 0; +// pthread_t th1, th2; +// intptr_t ret1, ret2; +// pthread_create(&th1, NULL, f1, NULL); +// pthread_create(&th2, NULL, f2, NULL); +// pthread_join(th1, (void **)&ret1); +// pthread_join(th2, (void **)&ret2); +// fail += (ret1 + ret2 == 0); +// } +// printf("%d/%d\n", fail, 1000000); +// } + +// volatile bool val; + +// void f1() { +// sleep(1); +// val = true; +// while (true) +// ; +// } + +// void f2() { +// while (!val) +// ; +// puts("thread 2 exit!"); +// } + +// int main() { +// bool val = false; + +// std::thread th1(f1); +// std::thread th2(f2); +// // std::thread th1([&] { +// // sleep(1); +// // val = true; +// // while (true) +// // ; +// // }); +// // std::thread th2([&] { +// // while (!val) +// // ; +// // puts("thread 2 exit!"); +// // }); + +// th1.join(); +// th2.join(); +// } \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bin_omp.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bin_omp.cpp" new file mode 100644 index 000000000..59ffe9a7d --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bin_omp.cpp" @@ -0,0 +1,57 @@ +#include "util.hpp" + +std::unique_ptr gen_data(u32 n, f32 min, f32 max) { + std::unique_ptr ret(new f32[n]); + for (u32 i = 0; i < n; ++i) { + ret[i] = min + (max - min) * (f32(rand()) / RAND_MAX); + } + return ret; +} + +i32 main() { + u32 bin, n, th; + f32 min, max; + + puts("Enter the number of bins"); + scanf("%d", &bin); + puts("Enter the minimum measurement"); + scanf("%f", &min); + puts("Enter the maximum measurement"); + scanf("%f", &max); + puts("Enter the number of data"); + scanf("%d", &n); + puts("Enter the number of threads"); + scanf("%d", &th); + + std::unique_ptr loc_cnt(new u32[bin * th]); + std::unique_ptr data = gen_data(n, min, max); + +#pragma omp parallel num_threads(th) + { +#pragma omp for + for (u32 tid = 0; tid < th; ++tid) { + u32 beg = n / th * tid, end = (tid == th - 1) ? n : n / th * (tid + 1); + u32 *x = loc_cnt.get() + bin * tid; + memset(x, 0, bin * sizeof(u32)); + for (u32 i = beg; i < end; ++i) { + ++x[u32((data[i] - min) / ((max - min) / bin))]; + } + } +#pragma omp for + for (u32 i = 0; i < bin; ++i) { + u32 sum = 0; + for (u32 j = 0; j < th; ++j) { + sum += loc_cnt[j * bin + i]; + } + loc_cnt[i] = sum; + } + } + + for (u32 i = 0; i < bin; ++i) { + f32 loc_min = min + (max - min) / bin * i; + f32 loc_max = min + (max - min) / bin * (i + 1); + printf("%.3f-%.3f:\t", loc_min, loc_max); + for (u32 j = 0; j < loc_cnt[i]; ++j) putchar('X'); + putchar('\n'); + } +} diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bin_pth.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bin_pth.cpp" new file mode 100644 index 000000000..0fb2bf37e --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/bin_pth.cpp" @@ -0,0 +1,79 @@ +#include "util.hpp" + +u32 bin, n, th; +f32 min, max; +std::unique_ptr loc_cnt; +std::unique_ptr data; +pthread_barrier_t calc, reduce[32]; + +std::unique_ptr gen_data(u32 n, f32 min, f32 max) { + std::unique_ptr ret(new f32[n]); + for (u32 i = 0; i < n; ++i) { + ret[i] = min + (max - min) * (f32(rand()) / RAND_MAX); + } + return ret; +} + +void *thread_fn(void *_tid) { + usize tid = (usize)_tid; + u32 beg = n / th * tid, end = (tid == th - 1) ? n : n / th * (tid + 1); + u32 *x = loc_cnt.get() + bin * tid; + memset(x, 0, bin * sizeof(u32)); + for (u32 i = beg; i < end; ++i) { + ++x[u32((data[i] - min) / ((max - min) / bin))]; + } + pthread_barrier_wait(&calc); + + for (u32 i = 0; tid + (1 << i) < th && !(tid >> i & 1); ++i) { + u32 *y = x + (bin << i); + for (u32 j = 0; j < bin; ++j) { + x[j] += y[j]; + } + pthread_barrier_wait(&reduce[i]); + } + return nullptr; +} + +i32 main() { + puts("Enter the number of bins"); + scanf("%d", &bin); + puts("Enter the minimum measurement"); + scanf("%f", &min); + puts("Enter the maximum measurement"); + scanf("%f", &max); + puts("Enter the number of data"); + scanf("%d", &n); + puts("Enter the number of threads"); + scanf("%d", &th); + + data = gen_data(n, min, max); + loc_cnt.reset(new u32[bin * th]); + + std::unique_ptr ths(new pthread_t[th]); + + pthread_barrier_init(&calc, nullptr, th); + for (u32 i = th, idx = 0; i != 1; ++idx, i = (i + 1) >> 1) { + pthread_barrier_init(&reduce[idx], nullptr, i >> 1); + } + for (usize i = 0; i < th; ++i) { + pthread_create(&ths[i], nullptr, thread_fn, (void *)i); + } + for (u32 i = 0; i < th; ++i) { + pthread_join(ths[i], nullptr); + } + + for (u32 i = 0; i < bin; ++i) { + f32 loc_min = min + (max - min) / bin * i; + f32 loc_max = min + (max - min) / bin * (i + 1); + printf("%.3f-%.3f:\t", loc_min, loc_max); + for (u32 j = 0; j < loc_cnt[i]; ++j) putchar('X'); + putchar('\n'); + } +#ifdef CHECK_SUM + u32 sum = 0; + for (u32 i = 0; i < bin; ++i) { + sum += loc_cnt[i]; + } + assert(sum == n); +#endif +} diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_omp_section.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_omp_section.cpp" new file mode 100644 index 000000000..7d365a01b --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_omp_section.cpp" @@ -0,0 +1,67 @@ +#include "util.hpp" + +usize seq_fib(usize n) { + if (n <= 2) { + return 1; + } + return seq_fib(n - 1) + seq_fib(n - 2); +} + +OmpMutex mu; +u32 remain; + +usize par_fib(usize n) { + if (n <= 20) { + return seq_fib(n); + } + bool spawn = false; + { + Locker _lk(&mu); + if (remain) { + --remain; + spawn = true; + } + } + if (spawn) { + usize lres, rres; +#pragma omp parallel sections + { +#pragma omp section + lres = par_fib(n - 1); +#pragma omp section + rres = par_fib(n - 2); + } + { + Locker _lk(&mu); + ++remain; + } + return lres + rres; + } else { + return par_fib(n - 1) + par_fib(n - 2); + } +} + +i32 main(i32 argc, i8 **argv) { + if (argc < 3) { + printf("usage %s [n] [th]\n", argv[0]); + exit(-1); + } + usize n = atoi(argv[1]); + remain = atoi(argv[2]); + + using namespace std::chrono; + auto now = high_resolution_clock::now; + + auto beg = now(); + + omp_set_dynamic(1); + omp_set_nested(1); + usize res = par_fib(n); + + printf("par res: %lu\n", res); + printf("par time: %.3f\n", duration(now() - beg).count()); + + beg = now(); + printf("seq res: %lu\n", seq_fib(n)); + printf("seq time: %.3f\n", duration(now() - beg).count()); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_omp_task.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_omp_task.cpp" new file mode 100644 index 000000000..af16971a2 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_omp_task.cpp" @@ -0,0 +1,46 @@ +#include "util.hpp" + +usize seq_fib(usize n) { + if (n <= 2) { + return 1; + } + return seq_fib(n - 1) + seq_fib(n - 2); +} + +usize par_fib(usize n) { + if (n <= 20) { + return seq_fib(n); + } + usize lres, rres; +#pragma omp task shared(lres) + lres = par_fib(n - 1); +#pragma omp task shared(rres) + rres = par_fib(n - 2); +#pragma omp taskwait + return lres + rres; +} + +i32 main(i32 argc, i8 **argv) { + if (argc < 3) { + printf("usage %s [n] [th]\n", argv[0]); + exit(-1); + } + usize n = atoi(argv[1]); + + using namespace std::chrono; + auto now = high_resolution_clock::now; + + auto beg = now(); + + usize res; +#pragma omp parallel num_threads(atoi(argv[2]) + 1) +#pragma omp single + res = par_fib(n); + + printf("par res: %lu\n", res); + printf("par time: %.3f\n", duration(now() - beg).count()); + + beg = now(); + printf("seq res: %lu\n", seq_fib(n)); + printf("seq time: %.3f\n", duration(now() - beg).count()); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_pth.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_pth.cpp" new file mode 100644 index 000000000..e3c57ab47 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/fib_pth.cpp" @@ -0,0 +1,61 @@ +#include "util.hpp" + +Mutex mu; +u32 remain; + +usize seq_fib(usize n) { + if (n <= 2) { + return 1; + } + return seq_fib(n - 1) + seq_fib(n - 2); +} + +void *par_fib(void *_n) { + usize n = (usize)_n; + if (n <= 20) { + return (void *)seq_fib(n); + } + bool spawn = false; + { + Locker _lk(&mu); + if (remain) { + --remain; + spawn = true; + } + } + + if (spawn) { + pthread_t r; + usize rres; + pthread_create(&r, nullptr, par_fib, (void *)(n - 2)); + usize lres = (usize)par_fib((void *)(n - 1)); + pthread_join(r, (void **)&rres); + { + Locker _lk(&mu); + ++remain; + } + return (void *)(lres + rres); + } else { + return (void *)((usize)par_fib((void *)(n - 1)) + (usize)par_fib((void *)(n - 2))); + } +} + +i32 main(i32 argc, i8 **argv) { + if (argc < 3) { + printf("usage %s [n] [th]\n", argv[0]); + exit(-1); + } + usize n = atoi(argv[1]); + remain = atoi(argv[2]); + + using namespace std::chrono; + auto now = high_resolution_clock::now; + + auto beg = now(); + printf("par res: %lu\n", (usize)par_fib((void *)n)); + printf("par time: %.3f\n", duration(now() - beg).count()); + + beg = now(); + printf("seq res: %lu\n", seq_fib(n)); + printf("seq time: %.3f\n", duration(now() - beg).count()); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/good_atomic.c" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/good_atomic.c" new file mode 100644 index 000000000..279125343 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/good_atomic.c" @@ -0,0 +1,36 @@ +#include +#include +#include +#include +#include + +#define ALL (10000000) +#define THREAD_NUM (4) +#define EACH (ALL / THREAD_NUM) + +intptr_t sum; +_Atomic intptr_t flag; + +void *inc(void *_tid) { + intptr_t tid = (intptr_t)_tid; + for (int i = 0; i < EACH; ++i) { + while (atomic_load_explicit(&flag, memory_order_acquire) != tid) + ; + ++sum; + intptr_t new_flag = atomic_load_explicit(&flag, memory_order_relaxed); + new_flag = (new_flag + 1) % THREAD_NUM; + atomic_store_explicit(&flag, new_flag, memory_order_release); + } + return NULL; +} + +int main() { + pthread_t th[THREAD_NUM]; + for (intptr_t i = 0; i < THREAD_NUM; ++i) { + pthread_create(&th[i], NULL, inc, (void *)i); + } + for (intptr_t i = 0; i < THREAD_NUM; ++i) { + pthread_join(th[i], NULL); + } + printf("%ld %d\n", sum, ALL); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/makefile" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/makefile" new file mode 100644 index 000000000..59250bd9b --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/makefile" @@ -0,0 +1,32 @@ +PTH_FLAGS = -Ofast -pthread -std=c++11 -Wall -Wextra +OMP_FLAGS = -Ofast -fopenmp -std=c++11 -Wall -Wextra + +all: 4_7_1 4_7_2 4_7_3 bin_pth bin_omp fib_pth fib_omp_task fib_omp_section + +4_7_1: 4_7_1.cpp util.hpp + g++ $(PTH_FLAGS) 4_7_1.cpp -o 4_7_1 + +4_7_2: 4_7_2.cpp util.hpp + g++ $(PTH_FLAGS) 4_7_2.cpp -o 4_7_2 + +4_7_3: 4_7_3.cpp util.hpp + g++ $(PTH_FLAGS) 4_7_3.cpp -o 4_7_3 + +bin_pth: bin_pth.cpp util.hpp + g++ $(PTH_FLAGS) bin_pth.cpp -o bin_pth + +bin_omp: bin_omp.cpp util.hpp + g++ $(OMP_FLAGS) bin_omp.cpp -o bin_omp + +fib_pth: fib_pth.cpp util.hpp + g++ $(PTH_FLAGS) fib_pth.cpp -o fib_pth + +fib_omp_task: fib_omp_task.cpp util.hpp + g++ $(OMP_FLAGS) fib_omp_task.cpp -o fib_omp_task + +fib_omp_section: fib_omp_section.cpp util.hpp + g++ $(OMP_FLAGS) fib_omp_section.cpp -o fib_omp_section + +clean: + rm 4_7_1 4_7_2 4_7_3 bin_pth bin_omp fib_pth fib_omp_task fib_omp_section + \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/report.md" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/report.md" new file mode 100644 index 000000000..37b1888f5 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/report.md" @@ -0,0 +1,820 @@ +--- +title: "高性能计算导论第四次作业" +author: "李晨昊 2017011466" +date: "2019-5-5" +output: + pdf_document: + latex_engine: xelatex + number_sections: yes + toc: yes + word_document: + toc: yes + html_document: + toc: yes +header-includes: \usepackage{ctex} \usepackage{graphicx} \usepackage{hyperref} \hypersetup { colorlinks = true, linkcolor=blue } \usepackage{multirow} +--- + +\newpage + +# Exercise 4.7 + +1. 一个生产者+一个消费者 + +生产者和消费者的逻辑完全无关,也不用放在循环里,所以没必要写成一个函数,就写成生产者+消费者两个函数即可。 + +```cpp +void* cons(void*) { + while (true) { // 循环直到成功消费一次为止 + // 本质上来讲这个锁的意义不是特别大 + // 只要能保证ok和msg的内存序的话,其实是不用上锁的 + Locker _lk(&mu); + if (ok) { + printf("get msg %s\n", msg); + break; + } + } + return nullptr; +} + +void* prod(void*) { // 生产一次就退出 + // 上锁可以保证消费者不会看到ok为true而msg未设置的状态 + // 但是其实只要正确地设置内存序,也是可以的 + Locker _lk(&mu); + strcpy(msg, "hello world!"); + ok = true; + return nullptr; +} + + // caller + pthread_t cons, prod; + pthread_create(&cons, nullptr, ::cons, nullptr); + pthread_create(&prod, nullptr, ::prod, nullptr); + + pthread_join(cons, nullptr); + pthread_join(prod, nullptr); +``` + +关于这里提到的内存序的问题,看起来是`volatile`试图解决的问题,实际上不是的。具体可以参见\hyperref[附]{附}。 + +2. 偶数生产者,奇数消费者 + +与1的逻辑完全一致,只是需要循环创建线程,而偶数号线程承担生产者工作,奇数号线程承担消费者工作。 + +```cpp +void *thread_fn(void *_tid) { + usize tid = (usize)_tid; + if (tid % 2 == 0) { // prod + while (true) { + Locker _lk(&mu); + if (!ok) { + sprintf(msg, "hello world from tid = %lu", tid); + ok = true; + break; + } + } + } else { // cons + while (true) { + Locker _lk(&mu); + if (ok) { + printf("%lu: get msg: %s\n", tid, msg); + ok = false; + break; + } + } + } + return nullptr; +} +``` + +3. 消费前一个,生产给后一个 + +```cpp +u32 recv, n; + +void *thread_fn(void *_tid) { + usize tid = (usize)_tid; + bool sended = false, recved = false; + while (!sended || !recved) { // 收到一次且发出一次时退出 + Locker _lk(&mu); + if (ok) { + if (!recved && recv == tid) { // 没收到过,且允许的接收方是自身 + printf("%lu: get msg: %s\n", tid, msg); + ok = false; + recved = true; + } + } else if (!sended) { + sprintf(msg, "hello world from tid = %lu", tid); + ok = true; + sended = true; + recv = (tid + 1) % n; // 设置允许的接收方为下一个线程 + } + } + return nullptr; +} +``` + +# Exercise 4.11 +1. 两个Delete操作同时执行 + +线程$T_0$删除头结点,线程$T_1$删除头结点下一个节点,执行顺序如下: + +```cpp +0: if (pred_p == NULL) // true +0: *head_p = curr_p->next; +0: free(curr_p); +1: if (pred_p == NULL) // false +1: pred_p->next = curr_p->next; // access freed memory +``` + +2. 一个Insert和一个Delete操作同时执行。 + +线程$T_0$删除头结点,线程$T_1$插入在原先链表中能排次小的节点,执行顺序如下: + +```cpp +0: if (pred_p == NULL) // true +0: *head_p = curr_p->next; +0: free(curr_p); +1: if (pred_p == NULL) // false +1: pred_p->next = temp_p; // access freed memory +``` + +3. 一个Member和一个Delete操作同时执行。 + +线程$T_0$删除头结点,线程$T_1$查找头结点,执行顺序如下: + +```cpp +0: if (pred_p == NULL) // true +0: *head_p = curr_p->next; +0: free(curr_p); +1: while (curr_p != NULL && curr_p->date < value) // access freed memory +``` + +4. 两个Insert同时执行。 + +线程$T_0$和$T_1$都插入在原先链表中能排最小的节点(这两个值的大小关系随意),执行顺序如下: + +```cpp +0: if (pred_p == NULL) // true +0: *head_p = temp_p; +1: if (pred_p == NULL) // true +1: *head_p = temp_p; // override T0's change +``` + +5. 一个Insert和一个Member同时执行。 + +线程$T_0$插入的节点正是线程$T_1$正在寻找的。严格来说,这种情况下的"正确"与否并不是很好界定,如果定义为:输出必须与操作序列按顺序到来时的输出完全一致,则如果插入先开始,但是查找结束时仍未成功插入,则这次"不存在"的返回结果是错误的;如果定义为:输出应该符合操作到来时刻的链表状态,则如果插入先开始,但是查找结束时已经成功插入,则这次"存在"的返回结果是错误的。 + +但是,还是有一个确定的标准可以判定,这两个操作同时执行是会出错的,即内存安全的标准。表面上来看,插入和查找同时执行是不会产生内存不安全的,但实际情况并非如此。 + +首先编译器可以进行代码重排,对于插入的这个代码片段: + +```cpp + ... +1: temp_p->next = curr_p; +2: if (pred_p == NULL) +3: *head_p = temp_p; +4: else +5: pred_p->next = temp_p; +``` + +2345行可能被调整到1行前执行,此时查找过程可能读取到一个`next`尚未赋值的节点,从而发生内存错误。 + +即使在语法层面指示了编译器不要这样重排,处理器执行过程中仍然不能保证内存读写顺序就是汇编代码的顺序。如果发生store-store重排,上面的情况仍然可能发生。幸运(?)的是,平时常用的`x86_64`架构下这种重排是不会发生的。但是这对别的架构就不一定成立了。具体可以参见\hyperref[附]{附}。 + +# Exercise 4.12 +不安全。 + +同时删除:如果在删除的查找阶段只使用读锁,可能会发生两个线程同时试图删除同一个节点,它们也都能顺利地找到这个节点,虽然后续的删除过程会分为一先一后,但是后执行的线程会对一个已经删除过的节点再删除一次,可能发生包括double free之类的各种错误。 + +同时插入:如果两个试图插入的节点都位于原链表的同两个节点中间,则它们会找到同一个插入位置,后续插入时不能保证这两个节点中,大的被放在小的后面。 + +一插入一删除:如果插入节点的位置在待删除节点之后,且删除操作先执行,则删除完后插入时,会访问被free了的内存。 + +# Exercise 5.4 +初值是对应运算的单位元。 + +| 操作 | 初值 | +|-----|------| +| && | true(1) | +| \|\|| false(0) | +| & | $\underbrace{11111...11111_2}_{\text{sizeof(T) * CHAR\_BIT}}$ | +| \| | 0 | +| ^ | 0 | + +# Exercise 5.5 +1. 计算完四次加法后得到$1008$,存入内存后四舍五入到$101*10^1$,所以输出为`1010.0` + +2. 线程0计算出$4.000$,线程1计算出$1003$,存入内存后,参与下一次加法的是$4.00$和$100*10^1$。加法结果为$1004$,保存后结果为$100*10^1$,输出结果为`1000.0`。 + +# Exercise 5.8 +只需利用等差数列求和公式即可: + +```cpp +#pragma omp parallel for + for (int i = 0; i < n; ++i) { + a[i] = i * (i + 1) / 2; + } +``` + +# Exercise 5.14 +1. 为了储存向量y,最少需要多少个缓存行? + +$y \in R^8$,最少需要$\frac{8 * 8}{64} = 1$个缓存行。 + +2. 为了储存向量y,最多需要多少个缓存行? + +$y$至多可以存在于两个连续的缓存行,故最多需要2个缓存行。 + +3. 如果缓存行的边界和双精度数的边界始终一致,那么一共有多少种方式将y中的元素分配给缓存行? + +按照$y$的第一个元素在其所在的缓存行中的8个位置分类,得一共有8种方式将y中的元素分配给缓存行。 + +4. 如果我们只考虑两个线程共享一个处理器,那么在我们的计算机上一共有多少种方式将四个线程分配给处理器?我们假定在同一个处理器上的内核共享缓存。 + +$\frac{C_4^2}{2} = 3$种。其中$C_4^2$是让一个处理器选择两个线程,除以二是因为选择某两个线程和选择除他们外的两个线程被视作相同的选择。 + +5. 在我们的例子中,有没有哪一种对向量的分配和对线程的分配方式不会产生伪共享?换句话说,对于处于不同处理器上的线程,他们各自要处理的向量的元素不会出现在同一个缓存行内? + +有。例如:$y$的前四个元素在一个缓存行,后四个元素在另一个;线程01位于核0上,它们负责计算$y$的前四个元素;线程23位于核1上,它们负责计算$y$的后四个元素。 + +6. 有多少种方法将向量元素分配给缓存行和将线程分配给处理器? + +$3 * 8 =24$ + +7. 在这些分配中,有多少会使得程序没有伪共享? + + 1. 若不是$y$的前四个元素在一个缓存行,后四个元素在另一个,则必有一个缓存行拥有超过4个$y$中元素,不能被一个核上的两个线程完全处理,从而一定会需要两个核的缓存。 + 2. (假定线程0,1,2,3分别负责$y$的01,23,45,67元素)若不是线程01在一个核,线程23在另一个核,则必定导致一个缓存行被两个核处理。 + + 由12知,有且仅有一种分配方法,就是5中提到的那种。 + +# Programming Assignments 5.1 + +我编写了`pthread`和`openmp`两个版本的~~(其实是因为没看清要求,先写了个`pthread`版的)~~。 + +`openmp`版非常简单,开启一个`omp parallel`块后,执行两次`omp for`,一次是统计每个线程内的桶的计数,一次是把每个桶的计数合并起来,两次`omp for`间用`omp barrier`作为路障,保证前一阶段执行完成。关键代码如下 + +```cpp +#pragma omp parallel num_threads(th) + { +#pragma omp for + for (u32 tid = 0; tid < th; ++tid) { + // 分配任务,这样的分配方法并不是最均匀的,但是比较容易写 + // 而且n一般远大于p,所以影响不大 + u32 beg = n / th * tid, end = (tid == th - 1) ? n : n / th * (tid + 1); + u32 *x = loc_cnt.get() + bin * tid; + memset(x, 0, bin * sizeof(u32)); + // 直接根据除法来计算位于哪个桶 + for (u32 i = beg; i < end; ++i) { + ++x[u32((data[i] - min) / ((max - min) / bin))]; + } + } +#pragma omp barrier +#pragma omp for + for (u32 i = 0; i < bin; ++i) { + // 使用临时变量保存一个桶在所有线程上的和 + // 确保false sharing不会发生(虽然我也相信编译器可以自动做出这个转换) + u32 sum = 0; + for (u32 j = 0; j < th; ++j) { + sum += loc_cnt[j * bin + i]; + } + loc_cnt[i] = sum; + } + } +``` + +`pthread`版本我写的稍复杂一些,我只开启了一次线程,在这一个函数中执行统计和树形规约的操作。在主线程中可以计算出树形规约的每一层有多少线程需要参与规约,然后每一层都使用`pthread_barrier`来同步这些线程。 + +```cpp + // 主线程中初始化路障,这个计算方法是我自己推导出来的 + // 每一层,剩下的线程数是上一层的线程数/2上取整 + // 参与计算的线程数是本层线程数/2下取整 + for (u32 i = th, idx = 0; i != 1; ++idx, i = (i + 1) >> 1) { + pthread_barrier_init(&reduce[idx], nullptr, i >> 1); + } + + ... + + u32 beg = n / th * tid, end = (tid == th - 1) ? n : n / th * (tid + 1); + u32 *x = loc_cnt.get() + bin * tid; + memset(x, 0, bin * sizeof(u32)); + for (u32 i = beg; i < end; ++i) { + ++x[u32((data[i] - min) / ((max - min) / bin))]; + } + pthread_barrier_wait(&calc); // 所有线程等待统计完成 + + // 树形求和 + // 退出条件为本层的加法的右手项超出了线程总数,或者本层的编号不能整除2^(高度) + for (u32 i = 0; tid + (1 << i) < th && !(tid >> i & 1); ++i) { + u32 *y = x + (bin << i); + for (u32 j = 0; j < bin; ++j) { + x[j] += y[j]; + } + pthread_barrier_wait(&reduce[i]); + } +``` + +简单测试一下,两边输入数据都为: + +``` +Enter the number of bins +10 +Enter the minimum measurement +0 +Enter the maximum measurement +100 +Enter the number of data +1000000000 +Enter the number of threads +10 +``` + +统计实际计算耗时,`openmp`版本为1.574s,`pthread`版本为1.382s。可见后者虽然稍复杂些,但是性能也好一些。 + +其实分析一下,二者的前半部分是一样的。对于后半部分,`openmp`版本理论耗时为$O(\frac{bin * th}{th} = bin)$,`pthread`版本理论耗时为$O(bin \log th)$,后者应该慢一些才对。实际结果并非如此,我猜测可能与访存的局部性优劣有关。`openmp`版本每个线程内部以步长$bin$来访问内存,虽然所有线程总体上访问的内存还是连续的,但是这也只能在L3 cache的层面体现出来,而L3 cache的速度已经不算很快了。而`pthread`版本每个线程内部都是以步长1访问内存,那么L1 cache和L2 cache的作用则能体现出来。 + +# Pthread Programming + +我尝试了一下直接将文字说明翻译成代码,发现这样的效率实在太低了,完全不能接受。因此我直接实现了一个比较优化的版本,主要改动了两点: + +1. $n$较小($n \le 20$)时调用串行版本 +2. 可以开启线程时,不开启两个,而是只开启一个,本线程充当另一个线程继续计算 + +关键代码如下: + +```cpp +Mutex mu; +u32 remain; + +usize seq_fib(usize n) { + if (n <= 2) { + return 1; + } + return seq_fib(n - 1) + seq_fib(n - 2); +} + +void *par_fib(void *_n) { + usize n = (usize)_n; + if (n <= 20) { + return (void *)seq_fib(n); + } + bool spawn = false; + { + Locker _lk(&mu); + if (remain) { + --remain; + spawn = true; + } + } + + if (spawn) { + pthread_t r; + usize rres; + pthread_create(&r, nullptr, par_fib, (void *)(n - 2)); + usize lres = (usize)par_fib((void *)(n - 1)); + pthread_join(r, (void **)&rres); + { + Locker _lk(&mu); + ++remain; + } + return (void *)(lres + rres); + } else { + return (void *)((usize)par_fib((void *)(n - 1)) + (usize)par_fib((void *)(n - 2))); + } +} +``` + +测试结果如下,注意表格中的`p`实际上是`th + 1`,其中`th`是输入参数,表示允许**新生成**的线程数: + +\newpage + +\begin{table}[!htbp] +\centering +\begin{tabular}{|c|c|c|c|c|c|c|} +\hline +\multicolumn{2}{|c|}{ } & \multicolumn{5}{c|}{n} \\ +\cline{3-7} +\multicolumn{2}{|c|}{ } & 30 & 34 & 38 & 42 & 46 \\ +\hline +\multicolumn{2}{|c|}{ seq\_fib } & 0.002 & 0.011 & 0.049 & 0.34 & 2.321 \\ +\hline +\multirow{6}*{par\_fib} & p = 1& 0.002 & 0.011 & 0.047 & 0.321 & 2.199 \\ +\cline{2-7} +& p = 2 & 0.002 & 0.008 & 0.043 & 0.201 & 1.366 \\ +\cline{2-7} +& p = 4 & 0.001 & 0.003 & 0.019 & 0.15 & 0.689 \\ +\cline{2-7} +& p = 8 & 0.001 & 0.003 & 0.015 & 0.096 & 0.443 \\ +\cline{2-7} +& p = 16 & 0.005 & 0.007 & 0.012 & 0.066 & 0.26 \\ +\cline{2-7} +& p = 24 & 0.009 & 0.009 & 0.012 & 0.055 & 0.196 \\ +\hline +\end{tabular} +\end{table} + +数据基本符合预期,有一点比较奇怪的地方是`p=1`时的`par_fib`略快于`seq_fib`。我猜测这可能是因为`par_fib`中调用了`seq_fib`,在这个地方编译器可以多展开几层函数调用,从而减少了递归的开销。 + +使用`openmp`来实现也是可行的,基本的逻辑都没变,调用的东西稍有改变而已: + +```cpp +OmpMutex mu; +u32 remain; + +usize par_fib(usize n) { + if (n <= 20) { + return seq_fib(n); + } + bool spawn = false; + { + Locker _lk(&mu); + if (remain) { + --remain; + spawn = true; + } + } + if (spawn) { + usize lres, rres; + // 开启两个并行的section,正常的实现应该只会启动一个新线程 +#pragma omp parallel sections + { +#pragma omp section + lres = par_fib(n - 1); +#pragma omp section + rres = par_fib(n - 2); + } + { + Locker _lk(&mu); + ++remain; + } + return lres + rres; + } else { + return par_fib(n - 1) + par_fib(n - 2); + } +} + + // caller + // 设置调度策略,经测试这样能显著变快,可能是因为计算两个fib的耗时差距很大 + // 如果动态分配线程的话等待时间可能更少 + omp_set_dynamic(1); + // 设置允许并行区域嵌套(默认不允许),否则最多只能在两个线程上运行 + omp_set_nested(1); + usize res = par_fib(n); +``` + +测试结果如下: + +\begin{table}[!htbp] +\centering +\begin{tabular}{|c|c|c|c|c|c|c|} +\hline +\multicolumn{2}{|c|}{ } & \multicolumn{5}{c|}{n} \\ +\cline{3-7} +\multicolumn{2}{|c|}{ } & 30 & 34 & 38 & 42 & 46 \\ +\hline +\multirow{6}*{par\_fib} & p = 1 & 0.002 & 0.012 & 0.048 & 0.327 & 2.239 \\ +\cline{2-7} +& p = 2 & 0.001 & 0.005 & 0.031 & 0.203 & 1.389 \\ +\cline{2-7} +& p = 4 & 0.001 & 0.004 & 0.021 & 0.137 & 0.766 \\ +\cline{2-7} +& p = 8 & 0.001 & 0.005 & 0.014 & 0.097 & 0.465 \\ +\cline{2-7} +& p = 16 & 0.006 & 0.007 & 0.012 & 0.066 & 0.272 \\ +\cline{2-7} +& p = 24 & 0.008 & 0.013 & 0.025 & 0.068 & 0.228 \\ +\hline +\end{tabular} +\end{table} + +可见使用`omp section`的性能与直接使用`pthread`差不多。 + +除此之外,`openmp`还有一套更高层的api:`omp task`,它可以动态的分发任务,只需要在一开始设定好最大线程数,则`openmp`的运行时库会帮我们处理好所需的逻辑。通过使用它,基本可以直接完成要求: + +```cpp +usize par_fib(usize n) { + if (n <= 20) { + return seq_fib(n); + } + usize lres, rres; +#pragma omp task shared(lres) + lres = par_fib(n - 1); +#pragma omp task shared(rres) + rres = par_fib(n - 2); +#pragma omp taskwait + return lres + rres; +} + + // caller + usize res; + // +1是因为包括了主线程在内 +#pragma omp parallel num_threads(atoi(argv[2]) + 1) +#pragma omp single + res = par_fib(n); +``` + +测试结果如下: + +\begin{table}[!htbp] +\centering +\begin{tabular}{|c|c|c|c|c|c|c|} +\hline +\multicolumn{2}{|c|}{ } & \multicolumn{5}{c|}{n} \\ +\cline{3-7} +\multicolumn{2}{|c|}{ } & 30 & 34 & 38 & 42 & 46 \\ +\hline +\multirow{6}*{par\_fib} & p = 1 & 0.002 & 0.012 & 0.049 & 0.336 & 2.302 \\ +\cline{2-7} +& p = 2 & 0.002 & 0.007 & 0.041 & 0.232 & 1.558 \\ +\cline{2-7} +& p = 4 & 0.001 & 0.006 & 0.041 & 0.241 & 1.077 \\ +\cline{2-7} +& p = 8 & 0.001 & 0.005 & 0.03 & 0.171 & 0.918 \\ +\cline{2-7} +& p = 16 & 0.002 & 0.005 & 0.031 & 0.154 & 1.215 \\ +\cline{2-7} +& p = 24 & 0.003 & 0.008 & 0.044 & 0.305 & 1.926 \\ +\hline +\end{tabular} +\end{table} + +性能测试结果不是很理想,并且观察到进程的CPU占用率并未达到$(th + 1) * 100\%$,这可能是`openmp`的调度策略不够理想导致的。 + +\newpage + +# 附:对课本/PPT上关于volatile的说法的质疑 + +\label{附} + +## volatile简介 + +选自\href{https://en.cppreference.com/w/cpp/language/cv}{cppreference: cv type qualifiers} + +> ... Every access (read or write operation, member function call, etc.) made through a glvalue expression of volatile-qualified type is treated as a **visible side-effect** for the purposes of optimization (that is, within a single thread of execution, volatile accesses cannot be optimized out or reordered with another visible side effect that is sequenced-before or sequenced-after the volatile access. This makes volatile objects **suitable for communication with a signal handler**, but **not with another thread of execution**, see std::memory_order). ... + +C/C++中`volatile`关键字的唯一语义是:每次对改变量的读写都会产生副作用,并且以此为基础指示(限制)编译器优化。`volatile`的直接推论如下: + +1. `volatile`变量的每次内存读写都需经过存储设备(包括内存和高速缓存)。 +2. `volatile`变量不能编译器被完全优化掉。 +3. `volatile`变量间读写顺序不能**在汇编层面**被改变。 + +## 课本/PPT上的用法的错误之处 + +看起来一切都很正常,对于课本上计算$\pi$的程序,似乎只需要把`sum`和`flag`均标记成`volatile`,即可保证生成的汇编代码不会随意调整临界区附近的代码顺序。 + +问题在于,这种顺序仅限于汇编代码层面,并不直接对应于CPU实际执行时的访存顺序。现代CPU都具备乱序执行的能力,即使汇编代码顺序看起来是合乎逻辑要求的,实际的执行结果也不完全保证符合汇编代码的顺序。 + +这其中,关于内存读写的顺序,不同平台的处理器的标准不尽相同。一些平台的内存顺序可参考\href{https://en.wikipedia.org/wiki/Memory_ordering}{wikipedia: memory ordering}。这里简单总结几点: + +1. 对于桌面电脑最常用的`x86_64`和`amd64`,唯一的内存重排序是load前的store可能被重排到load后 + - 所以对于计算$\pi$的程序,临界区的语义是可以保证的 +2. 但是很多其它平台则没有这么严格的内存序,例如手机最常用的`arm64`平台。 + - 由于store-store乱序的存在,`flag`可能在`sum`前被更新,这样一来对于`sum`的更新就位于临界区外了 + +我编写了一个简单的程序来证明这个理论结果,简单起见,也为了避免浮点误差对结果的干扰,这里就测试一系列整数的自增: + +```c +#include +#include +#include + +#define ALL (10000000) +#define THREAD_NUM (4) +#define EACH (ALL / THREAD_NUM) + +volatile intptr_t sum; +volatile intptr_t flag; + +void *inc(void *_tid) { + intptr_t tid = (intptr_t)_tid; + for (int i = 0; i < EACH; ++i) { + while (flag != tid) + ; + ++sum; + flag = (flag + 1) % THREAD_NUM; + } + return NULL; +} + +int main() { + pthread_t th[THREAD_NUM]; + for (intptr_t i = 0; i < THREAD_NUM; ++i) { + pthread_create(&th[i], NULL, inc, (void *)i); + } + for (intptr_t i = 0; i < THREAD_NUM; ++i) { + pthread_join(th[i], NULL); + } + printf("%ld %d\n", sum, ALL); +} +``` + +其中主线程生成四个线程,每个线程对全局的`sum`自增`EACH`次。`sum`和`flag`均用`volatile`标记。 + +在桌面平台(`AMD R7-2700, gcc 8.3.0, -O3`)上测试10次,结果都是正确的。在手机平台(`OnePlus3 MSM8996, aarch64-linux-gnu-gcc 8.3.0, -O3`)上测试10次,结果如下: + +\begin{center} +9999952 10000000 \\ +9998286 10000000 \\ +9999807 10000000 \\ +9999920 10000000 \\ +9997961 10000000 \\ +9997514 10000000 \\ +9997219 10000000 \\ +9999809 10000000 \\ +9999899 10000000 \\ +9999654 10000000 \\ +\end{center} + +最高的出错率达到了$\frac{10000000 - 9997219}{10000000} = 0.02781\%$。虽然这数字并不算大,但它毕竟不是0。 + +这里分别截取一段汇编来分析一下: + +`amd64`平台: + +```x86asm + ... + # 把全局变量flag读到rax +1 mov 0x2e39(%rip),%rax # 4058 + # 比较flag与tid +2 cmp %rdi,%rax + # 不等,则跳回1 +3 jne 1218 + # 全局变量sum自增1 +4 mov 0x2e25(%rip),%rax # 4050 +5 add $0x1,%rax +6 mov %rax,0x2e1a(%rip) # 4050 + # 更新flag,由于编译器的除/模常数优化,生成了一大堆代码 +7 mov 0x2e1b(%rip),%rax # 4058 +8 add $0x1,%rax +9 cqto +10 shr $0x3e,%rdx +11 add %rdx,%rax +12 and $0x3,%eax +13 sub %rdx,%rax + # 把更新好的flag写回全局变量 +14 mov %rax,0x2e01(%rip) # 4058 + # 外层循环更新和跳转 +15 sub $0x1,%ecx +16 jne 1218 + ... +``` + +`arm64`平台: + +```armasm + ... + # 辅助计算sum和flag的地址的一堆代码 + # 计算好后sum地址存在x5中,flag地址存在x2中 +1 adrp x2, 496000 +2 adrp x5, 496000 +3 mov w4, #0x25a0 // #9632 +4 add x2, x2, #0x300 +5 add x5, x5, #0x2f8 +6 movk w4, #0x26, lsl #16 +7 ldr x1, [x2] + # 比较flag与tid +8 cmp x1, x0 + # 不等,则跳回7 +9 b.ne 400678 // b.any + # 全局变量sum自增1 + # 12可以重排到19之后,从而失去保护 +10 ldr x1, [x5] +11 add x1, x1, #0x1 +12 str x1, [x5] + # 更新flag,由于编译器的除/模常数优化,生成了一大堆代码 +13 ldr x1, [x2] +14 add x1, x1, #0x1 +15 negs x3, x1 +16 and x1, x1, #0x3 +17 and x3, x3, #0x3 +18 csneg x1, x1, x3, mi // mi = first + # 把更新好的flag写回全局变量 +19 str x1, [x2] + # 外层循环更新和跳转 +20 subs w4, w4, #0x1 +21 b.ne 400678 // b.any + ... +``` + +基本上可以说是`C`语句和汇编语句间一一对应,看不出`volatile`的功能,**而这也恰好就是它的功能**。但是很可惜,对于`arm64`这种内存序较弱的架构,`volatile`的功能用来实现临界区是不够的,就算已经是这么清晰的汇编顺序,最终结果仍然可能是错的。 + +## 可行的方案 + +`cppreference`中已经写的比较明确了:`volatile`并不适用于多线程间的通信。它当然不会让正确的程序变错,但是也不能让错误的程序变对;它只能让正确的程序变慢,或者让错误的程序变得看起来正确。然而我们一般并不认为"看起来正确"比"错误"要更优越一些。 + +回到最初的需求,本质上是希望能够保证"变量间读写顺序不能被改变",而非仅仅"在汇编层面"不改变。为了告知处理器这一需求,需要使用特殊的机器指令,原有的指令无法表达这个语义。 + +`C11`和`C++11`中引入了原子类型和内存序,它们可以表达这种语义。当然,一个可行的方案是直接把`sum`定义成原子类型,然后调用原子自增即可,但是这样并不能用来实现所有的临界区。这里还是使用`flag`来表示临界区的语义,虽然复杂一些,但也更通用一些: + +```c +#include + +intptr_t sum; +_Atomic intptr_t flag; + +void *inc(void *_tid) { + intptr_t tid = (intptr_t)_tid; + for (int i = 0; i < EACH; ++i) { +1 while (atomic_load_explicit(&flag, memory_order_acquire) != tid) + ; +2, 3 ++sum; +4 intptr_t new_flag = atomic_load_explicit(&flag, memory_order_relaxed); + new_flag = (new_flag + 1) % THREAD_NUM; +5 atomic_store_explicit(&flag, new_flag, memory_order_release); + } + return NULL; +} +``` + +完整代码见`good_atomic.c`。 + +总共有5个读写内存的地点,其中1处的`memory_order_acquire`和5处的`memory_order_release`保证了1一定最先执行,5一定最后执行,中间的234还有一些自由调整的余地,但是这都不会影响到临界区的语义了。 + +在我的电脑和手机上又都分别进行了10次测试,结果都是正确的。这的确解决了问题,那么时间开销是否会比原来的版本高呢?依然可以分析汇编: + +`amd64`平台: + +```x86asm + mov 0x2e39(%rip),%rax # 4058 + cmp %rdi,%rax + jne 1218 + # 用一条addq指令代替了原先的读,加,存 + addq $0x1,0x2e24(%rip) # 4050 + mov 0x2e25(%rip),%rax # 4058 + add $0x1,%rax + cqto + shr $0x3e,%rdx + add %rdx,%rax + and $0x3,%eax + sub %rdx,%rax + mov %rax,0x2e0b(%rip) # 4058 + sub $0x1,%ecx + jne 1218 +``` +因为在`amd64`平台上这样的内存序要求是天然成立的,所以并不需要任何特别的指令。不仅如此,由于去掉了`sum`变量的`volatile`,编译器不再**错认为**对于`sum`读写有副作用,从而原先的读,加,存被用一条`addq`指令代替了,性能反而更好。值得注意的是,虽然有`addq`这种看起来很"原子"的指令,但实际上这条指令本身不是原子操作,多线程同时执行它仍然是存在竞争条件的,具体可参见\href{https://stackoverflow.com/questions/39393850/can-num-be-atomic-for-int-num/39394630#39394630}{StackOverFlow: Can num++ be atomic for 'int num'?}。 + +`arm64`平台: + +```armasm + adrp x2, 496000 + adrp x5, 496000 + mov w4, #0x25a0 // #9632 + add x2, x2, #0x300 + add x5, x5, #0x2f8 + movk w4, #0x26, lsl #16 + # 从ldr指令变成了ldar,表示acquire语义 + # LDAR Wt, [base{,#0}] + # Load-Acquire Register: loads a word from memory addressed by base to Wt. + ldar x1, [x2] + cmp x1, x0 + b.ne 400678 // b.any + # 对sum的操作没有任何变化 + ldr x1, [x5] + add x1, x1, #0x1 + str x1, [x5] + # 对flag的relaxed load,还是普通的ldr + ldr x1, [x2] + add x1, x1, #0x1 + negs x3, x1 + and x1, x1, #0x3 + and x3, x3, #0x3 + csneg x1, x1, x3, mi // mi = first + # 从str指令变成了stlr,表示release语义 + # STLR Wt, [base{,#0}] + # Store-Release Register: stores a word from Wt to memory addressed by base. + stlr x1, [x2] + subs w4, w4, #0x1 + b.ne 400678 // b.any +``` + +在`arm64`上,情况反过来,没有因为去掉`volatile`而得到了额外的优化机会,但是为一次`acquire load`和一次`release store`专门生成了对应的指令,预计速度可能会有所下降。 + +为了验证我的预测,测试了运行10次的总时间如下: + +\begin{center} +\begin{tabular}{|c|c|c|} +\hline + & amd64 & arm64 \\ +\hline +volatile & 0m15.561s & 0m27.67s \\ +\hline +atomic & 0m14.158s & 0m21.77s \\ +\hline +\end{tabular} +\end{center} + +结果发现`arm64 + atomic`的组合性能并没有下降,反而显著提升了。我猜测这可能是由于语义的改变导致它在忙等待中浪费的时间更少。 + +## 结论 + +多线程编程时,可以根据自身的需求来选择合适的同步手段,如锁,信号量,条件变量等;如果临界区是简单的自增等操作,可以直接用原子变量实现;如果临界区没这么简单,但是也相当简单,可以考虑自己用原子变量+忙等待来实现临界区,开销可能会比使用线程库提供的锁小一些。 + +如果不在乎可移植性,或者对程序的运行环境非常确定,又由于某种原因用不了`C11`和`C++11`中的新feature,可以考虑使用内嵌汇编来指示一个内存屏障: + +```cpp + asm volatile("" ::: "memory"); +``` + +这并不会生成任何指令,只是会阻止汇编层面的重排。相比于`volatile`,它不会让编译器错认为对于变量的读写有副作用(就如同之前的`sum`的例子),从而会有更多的优化机会。 + +无论如何,不要(为了同步的目的,在`C`/`C++`中)用`volatile`。 diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/report.pdf" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/report.pdf" new file mode 100644 index 000000000..6a8e5c4e4 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/report.pdf" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/run.py" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/run.py" new file mode 100644 index 000000000..191cc60cb --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/run.py" @@ -0,0 +1,24 @@ +from __future__ import print_function +import os +import subprocess +import re +import sys + +min_seqs = [] + + +for p in (1, 2, 4, 8, 16, 24): + for n in (30, 34, 38, 42, 46): + min_par = 1e9 + min_seq = 1e9 + for r in range(3): + o = subprocess.check_output(['./fib_omp_task', str(n), str(p-1)]) + o = o.split() + min_par = min(min_par, float(o[1][10:])) + min_seq = min(min_seq, float(o[3][10:])) + print(min_par, end=' ') + if p == 1: + min_seqs.append(min_seq) + print() + +print(min_seqs) diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/util.hpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/util.hpp" new file mode 100644 index 000000000..f004b181d --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw4/util.hpp" @@ -0,0 +1,44 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +using i8 = char; +using u8 = unsigned char; +using i32 = int; +using u32 = unsigned; +using usize = std::size_t; +using f32 = float; + +struct Mutex { + pthread_mutex_t inner; + Mutex() { pthread_mutex_init(&inner, nullptr); } + Mutex(const Mutex &) = delete; + Mutex &operator=(const Mutex &) = delete; + void lock() { pthread_mutex_lock(&inner); } + void unlock() { pthread_mutex_unlock(&inner); } + ~Mutex() { pthread_mutex_destroy(&inner); } +}; + +struct OmpMutex { + omp_lock_t inner; + OmpMutex() { omp_init_lock(&inner); } + OmpMutex(const Mutex &) = delete; + OmpMutex &operator=(const OmpMutex &) = delete; + void lock() { omp_set_lock(&inner); } + void unlock() { omp_unset_lock(&inner); } + ~OmpMutex() { omp_destroy_lock(&inner); } +}; + +template +struct Locker { + M *inner; + Locker(M *mu) : inner(mu) { inner->lock(); } + Locker(const Locker &) = delete; + Locker &operator=(const Locker &) = delete; + ~Locker() { inner->unlock(); } +}; \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/1_1.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/1_1.cpp" new file mode 100644 index 000000000..0664f0ff4 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/1_1.cpp" @@ -0,0 +1,120 @@ +#include +#include "util.hpp" + +std::pair, std::unique_ptr> gen(u32 n) { + std::unique_ptr mat(new f64[n * n]); + std::unique_ptr vec(new f64[n]); + for (u32 i = 0; i < n; ++i) { + for (u32 j = 0; j < n; ++j) { + mat[i * n + j] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + } + for (u32 i = 0; i < n; ++i) { + vec[i] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + return {std::move(mat), std::move(vec)}; +} + +std::unique_ptr serial(const std::unique_ptr &mat, const std::unique_ptr &vec, u32 n) { + std::unique_ptr ret(new f64[n]); + for (u32 i = 0; i < n; ++i) { + f64 tmp = 0.0; + for (u32 j = 0; j < n; ++j) { + tmp += mat[i * n + j] * vec[j]; + } + ret[i] = tmp; + } + return ret; +} + +void run(u32 n) { + i32 pid, nproc; + f64 beg[2], elapsed[2]; + + MPI_Init(nullptr, nullptr); + MPI_Comm_rank(MPI_COMM_WORLD, &pid); + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + + if (n % nproc != 0) { + printf("n must be divisible by nproc, but n = %d, nproc= %d\n", n, nproc); + exit(-1); + } + + u32 each = n / nproc; + std::unique_ptr loc_mat(new f64[n * each]), loc_vec(new f64[each]), ser_res; + + MPI_Datatype vec_t, col_t; + // vec_t is a type like struct { f64[each]@0, f64[each]@(n * sizeof(f64)), ..., f64[each]@((n * n - 1) * sizeof(f64))} + MPI_Type_vector(n, each, n, MPI_DOUBLE, &vec_t); + // col_t is a type like vec_t, but it behaves like sizeof(col_t) == each * sizeof(f64) in MPI_Scatter + MPI_Type_create_resized(vec_t, 0, each * sizeof(f64), &col_t); + MPI_Type_commit(&col_t); + std::unique_ptr size(new int[nproc]), off(new int[nproc]); + for (u32 i = 0; i < nproc; ++i) { + size[i] = 1; + off[i] = i; + } + + decltype(gen(n)) mat_vec(nullptr, nullptr); + if (pid == 0) { + mat_vec = gen(n); + ser_res = serial(mat_vec.first, mat_vec.second, n); + } + + MPI_Barrier(MPI_COMM_WORLD); + beg[0] = MPI_Wtime(); + + MPI_Scatterv(mat_vec.first.get(), size.get(), off.get(), col_t, loc_mat.get(), n * each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Scatter(mat_vec.second.get(), each, MPI_DOUBLE, loc_vec.get(), each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + MPI_Barrier(MPI_COMM_WORLD); + beg[1] = MPI_Wtime(); + + std::unique_ptr y(new f64[n]); + for (u32 i = 0; i < n; ++i) { + f64 tmp = 0.0; + for (u32 j = 0; j < each; ++j) { + tmp += loc_mat[i * each + j] * loc_vec[j]; + } + y[i] = tmp; + } + + beg[1] = MPI_Wtime() - beg[1]; + + if (pid == 0) { + std::unique_ptr res(new f64[n]); + MPI_Reduce(y.get(), res.get(), n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + beg[0] = MPI_Wtime() - beg[0]; + + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + printf("prod(+MPI_Scatter): %.5lfs\n", elapsed[0]); + printf("prod(-MPI_Scatter): %.5lfs\n", elapsed[1]); + + beg[0] = MPI_Wtime(); + + f64 l2 = 0.0; +#pragma omp parallel for reduction(+ : l2) + for (u32 i = 0; i < n; ++i) { + l2 += (ser_res[i] - res[i]) * (ser_res[i] - res[i]); + } + l2 = sqrt(l2); + + beg[0] = MPI_Wtime() - beg[0]; + printf("calc 2 norm: %.5lfs\n", beg[0]); + printf("2 norm: %.20lf\n", l2); + } else { + MPI_Reduce(y.get(), nullptr, n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + beg[0] = MPI_Wtime() - beg[0]; + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + } + MPI_Type_free(&col_t); + MPI_Finalize(); +} + +i32 main(i32 argc, i8 **argv) { + if (argc < 2) { + printf("usage: %s [n]\n", argv[0]); + exit(1); + } + run(atoi(argv[1])); +} diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/1_2.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/1_2.cpp" new file mode 100644 index 000000000..cdf88e214 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/1_2.cpp" @@ -0,0 +1,131 @@ +#include +#include "util.hpp" + +std::pair, std::unique_ptr> gen(u32 n) { + std::unique_ptr mat(new f64[n * n]); + std::unique_ptr vec(new f64[n]); + for (u32 i = 0; i < n; ++i) { + for (u32 j = 0; j < n; ++j) { + mat[i * n + j] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + } + for (u32 i = 0; i < n; ++i) { + vec[i] = -1 + rand() * 1.0 / RAND_MAX * 2; + } + return {std::move(mat), std::move(vec)}; +} + +std::unique_ptr serial(const std::unique_ptr &mat, const std::unique_ptr &vec, u32 n) { + std::unique_ptr ret(new f64[n]); + for (u32 i = 0; i < n; ++i) { + f64 tmp = 0.0; + for (u32 j = 0; j < n; ++j) { + tmp += mat[i * n + j] * vec[j]; + } + ret[i] = tmp; + } + return ret; +} + +void run(int n) { + i32 pid, nproc; + f64 beg[2], elapsed[2]; + + MPI_Init(nullptr, nullptr); + MPI_Comm_rank(MPI_COMM_WORLD, &pid); + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + + u32 sq_nproc = sqrtf(nproc); + if (sq_nproc * sq_nproc != nproc) { + printf("nproc must be complete square number, but nproc= %d\n", nproc); + exit(-1); + } + if (n % sq_nproc != 0) { + printf("n must be divisible by sqrt(nproc), but n = %d, nproc= %d\n", n, nproc); + exit(-1); + } + + u32 each = n / sq_nproc; + std::unique_ptr loc_mat(new f64[each * each]), vec(nullptr), ser_res; + + MPI_Datatype vec_t, block_t; + MPI_Type_vector(each, each, n, MPI_DOUBLE, &vec_t); + MPI_Type_create_resized(vec_t, 0, each * sizeof(f64), &block_t); + MPI_Type_commit(&block_t); + std::unique_ptr size(new i32[nproc]), off(new i32[nproc]); + for (u32 i = 0; i < sq_nproc; ++i) { + for (u32 j = 0; j < sq_nproc; ++j) { + size[i * sq_nproc + j] = 1; + off[i * sq_nproc + j] = i * n + j; + } + } + + std::unique_ptr mat(nullptr); + if (pid == 0) { + auto mat_vec = gen(n); + ser_res = serial(mat_vec.first, mat_vec.second, n); + mat = std::move(mat_vec.first); + vec = std::move(mat_vec.second); + } else { + vec.reset(new f64[n]); + } + + MPI_Barrier(MPI_COMM_WORLD); + beg[0] = MPI_Wtime(); + + MPI_Scatterv(mat.get(), size.get(), off.get(), block_t, loc_mat.get(), each * each, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(vec.get(), n, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + MPI_Barrier(MPI_COMM_WORLD); + beg[1] = MPI_Wtime(); + + std::unique_ptr y(new f64[n]{}); + int vec_off = pid % sq_nproc * each, y_off = pid / sq_nproc * each; + for (int i = 0; i < each; ++i) { + f64 tmp = 0.0; + for (int j = 0; j < each; ++j) { + tmp += loc_mat[i * each + j] * vec[vec_off + j]; + } + y[y_off + i] = tmp; + } + + beg[1] = MPI_Wtime() - beg[1]; + + if (pid == 0) { + std::unique_ptr res(new f64[n]); + MPI_Reduce(y.get(), res.get(), n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + beg[0] = MPI_Wtime() - beg[0]; + + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + printf("prod(+MPI_Scatter): %.5lfs\n", elapsed[0]); + printf("prod(-MPI_Scatter): %.5lfs\n", elapsed[1]); + + beg[0] = MPI_Wtime(); + + f64 l2 = 0.0; +#pragma omp parallel for reduction(+ : l2) + for (u32 i = 0; i < n; ++i) { + l2 += (ser_res[i] - res[i]) * (ser_res[i] - res[i]); + } + l2 = sqrt(l2); + + beg[0] = MPI_Wtime() - beg[0]; + printf("calc 2 norm: %.5lfs\n", beg[0]); + printf("2 norm: %.20lf\n", l2); + } else { + MPI_Reduce(y.get(), nullptr, n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + beg[0] = MPI_Wtime() - beg[0]; + MPI_Reduce(beg, elapsed, 2, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + } + + MPI_Type_free(&block_t); + MPI_Finalize(); +} + +i32 main(i32 argc, i8 **argv) { + if (argc < 2) { + printf("Usage: %s n\n", argv[0]); + return -1; + } + run(atoi(argv[1])); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/2.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/2.cpp" new file mode 100644 index 000000000..f002e80f4 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/2.cpp" @@ -0,0 +1,113 @@ +#include +#include +#include +#include "util.hpp" + +template +struct ThreadPool { + std::unique_ptr ths; + u32 n; + std::atomic no_more; + Mutex mu; + CondVar cv; + std::queue q; + + static void *task(void *arg) { + ThreadPool *self = (ThreadPool *)arg; + while (true) { + self->mu.lock(); + // Mesa monitor + while (self->q.empty() && !self->no_more) { + self->cv.wait(&self->mu); + } + if (self->q.empty()) { // no more + self->mu.unlock(); + break; + } + T t = self->q.front(); + self->q.pop(); + self->mu.unlock(); + t(); + } + return nullptr; + } + + ThreadPool(u32 n) : ths(new pthread_t[n]), n(n), no_more(false) { + for (u32 i = 0; i < n; ++i) { + pthread_create(&ths[i], nullptr, task, this); + } + } + + // this will wait until all tasks are finished + ~ThreadPool() { + no_more = true; + cv.notify_all(); + for (u32 i = 0; i < n; ++i) { + pthread_join(ths[i], nullptr); + } + } + + void push(T t) { + mu.lock(); + q.push(t); + cv.notify_one(); + mu.unlock(); + } +}; + +Mutex mu; +std::list lst; + +struct ListOp { + u32 ty, val; + + void operator()() { + Locker _lk(&mu); + switch (ty) { + case 0: + case 1: + printf("insert %d\n", val); + lst.push_back(val); + break; + case 2: { + auto it = std::find(lst.begin(), lst.end(), val); + if (it != lst.end()) { + lst.erase(it); + printf("remove %d succeed\n", val); + } else { + printf("remove %d fail\n", val); + } + } break; + case 3: { + auto it = std::find(lst.begin(), lst.end(), val); + if (it != lst.end()) { + printf("find %d succeed\n", val); + } else { + printf("find %d fail\n", val); + } + } break; + case 4: + printf("display ["); + for (u32 x : lst) { + printf("%d, ", x); + } + printf("]\n"); + break; + default: + __builtin_unreachable(); + } + } +}; + +i32 main(i32 argc, i8 **argv) { + if (argc != 3) { + printf("usage: %s [th] [task]\n", argv[0]); + exit(1); + } + u32 th = atoi(argv[1]); + u32 task = atoi(argv[2]); + ThreadPool pool(th); + for (u32 i = 0; i < task; ++i) { + pool.push({u32(rand()) % 5, u32(rand()) % 1000}); + } +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/3.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/3.cpp" new file mode 100644 index 000000000..15c18f377 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/3.cpp" @@ -0,0 +1,174 @@ +#include +#include "util.hpp" + +void count(i32 *a, u32 n) { + i32 *aux = (i32 *)malloc(n * sizeof(i32)); + for (u32 i = 0; i < n; ++i) { + u32 cnt = 0; + for (u32 j = 0; j < i; ++j) { + cnt += a[j] <= a[i]; + } + for (u32 j = i; j < n; ++j) { + cnt += a[j] < a[i]; + } + aux[cnt] = a[i]; + } + memcpy(a, aux, n * sizeof(i32)); + free(aux); +} + +void count_par(i32 *a, u32 n, u32 th) { + i32 *aux = (i32 *)malloc(n * sizeof(i32)); +#pragma omp parallel for num_threads(th) + for (u32 i = 0; i < n; ++i) { + u32 cnt = 0; + for (u32 j = 0; j < i; ++j) { + cnt += a[j] <= a[i]; + } + for (u32 j = i; j < n; ++j) { + cnt += a[j] < a[i]; + } + aux[cnt] = a[i]; + } + memcpy(a, aux, n * sizeof(i32)); + free(aux); +} + +namespace detail { +i32 *a, *aux; +u32 n, th, each; +pthread_t *ths; +pthread_barrier_t calc, reduce[32]; + +void merge(i32 *first, i32 *mid, i32 *last, i32 *aux) { + i32 *pos1 = first, *pos2 = mid, *pos = aux; + while (pos1 != mid && pos2 != last) { + i32 a = *pos1, b = *pos2; + i32 cmp = b < a; + *pos++ = cmp ? b : a; + pos1 += cmp ^ 1, pos2 += cmp; + } + memcpy(pos, pos1, sizeof(i32) * (mid - pos1)), pos += mid - pos1; + memcpy(first, aux, sizeof(i32) * (pos - aux)); +} + +void *thread_fn(void *_tid) { + u32 tid = (u32)(usize)_tid; + u32 beg = each * tid, end = (tid == th - 1) ? n : each * (tid + 1); + i32 *a_off = a + beg, *aux_off = aux + beg; + for (u32 i = beg; i < end; ++i) { + u32 cnt = 0; + for (u32 j = beg; j < i; ++j) { + cnt += a[j] <= a[i]; + } + for (u32 j = i; j < end; ++j) { + cnt += a[j] < a[i]; + } + aux_off[cnt] = a[i]; + } + memcpy(a_off, aux_off, sizeof(i32) * (end - beg)); + pthread_barrier_wait(&calc); + for (u32 i = 0; tid + (1 << i) < th && !(tid >> i & 1); ++i) { + u32 right_tid = tid + (1 << (i + 1)); + u32 right = (right_tid >= th) ? n : each * right_tid; + merge(a_off, a_off + (each << i), a + right, aux_off); + pthread_barrier_wait(&reduce[i]); + } + return nullptr; +} +} // namespace detail + +void count_llx(i32 *a, u32 n, u32 th) { + detail::a = a; + detail::aux = (i32 *)malloc(n * sizeof(i32)); + detail::n = n; + detail::th = th; + detail::each = n / th; + pthread_barrier_init(&detail::calc, nullptr, th); + u32 reduce_idx = 0; + for (u32 i = th; i != 1; ++reduce_idx, i = (i + 1) >> 1) { + pthread_barrier_init(&detail::reduce[reduce_idx], nullptr, i >> 1); + } + detail::ths = (pthread_t *)malloc(th * sizeof(pthread_t)); + for (usize i = 0; i < th; ++i) { + pthread_create(&detail::ths[i], nullptr, detail::thread_fn, (void *)i); + } + for (u32 i = 0; i < th; ++i) { + pthread_join(detail::ths[i], nullptr); + } + free(detail::ths); + pthread_barrier_destroy(&detail::calc); + for (u32 i = 0; i < reduce_idx; ++i) { + pthread_barrier_destroy(&detail::reduce[i]); + } + free(detail::aux); +} + +void radix(i32 *_a, u32 n) { + const u32 U = 256; + u32 cnt[U], *a = (u32 *)_a; + u32 *aux = (u32 *)malloc(n * sizeof(u32)); + memset(cnt, 0, sizeof cnt); + for (u32 i = 0; i < n; ++i) ++cnt[a[i] & U - 1]; + for (u32 i = 1; i < U; ++i) cnt[i] += cnt[i - 1]; + for (u32 i = n - 1; ~i; --i) aux[--cnt[a[i] & U - 1]] = a[i]; + memset(cnt, 0, sizeof cnt); + for (u32 i = 0; i < n; ++i) ++cnt[aux[i] >> 8 & U - 1]; + for (u32 i = 1; i < U; ++i) cnt[i] += cnt[i - 1]; + for (u32 i = n - 1; ~i; --i) a[--cnt[aux[i] >> 8 & U - 1]] = aux[i]; + memset(cnt, 0, sizeof cnt); + for (u32 i = 0; i < n; ++i) ++cnt[a[i] >> 16 & U - 1]; + for (u32 i = 1; i < U; ++i) cnt[i] += cnt[i - 1]; + for (u32 i = n - 1; ~i; --i) aux[--cnt[a[i] >> 16 & U - 1]] = a[i]; + memset(cnt, 0, sizeof cnt); + for (u32 i = 0; i < n; ++i) ++cnt[aux[i] >> 24 & U - 1]; + // notice this loop + for (u32 i = U / 2 + 1; i < U + U / 2; ++i) cnt[i & U - 1] += cnt[(i - 1) & U - 1]; + for (u32 i = n - 1; ~i; --i) a[--cnt[aux[i] >> 24 & U - 1]] = aux[i]; + free(aux); +} + +i32 main(i32 argc, i8 **argv) { + if (argc < 3) { + printf("usage: %s [algo] [n] \n", argv[0]); + puts("algo: 0 => count, 1 => count_par, 2 => count_llx, 3 => qsort, 4 => std::sort, 5 => radix"); + exit(1); + } + u32 algo = atoi(argv[1]); + u32 n = atoi(argv[2]); + std::unique_ptr a(new i32[n]); + for (u32 i = 0; i < n; ++i) { + a[i] = rand(); + } + + auto now = std::chrono::high_resolution_clock::now; + auto beg = now(); + switch (algo) { + case 0: + count(a.get(), n); + break; + case 1: + count_par(a.get(), n, argc == 4 ? atoi(argv[3]) : omp_get_max_threads()); + break; + case 2: + count_llx(a.get(), n, argc == 4 ? atoi(argv[3]) : omp_get_max_threads()); + break; + case 3: + // this may overflow in a general case + // but rand() always return a integer in [0, RAND_MAX], so no problem here + qsort(a.get(), n, sizeof(i32), [](const void *l, const void *r) { return *(i32 *)l - *(i32 *)r; }); + break; + case 4: + std::sort(a.get(), a.get() + n); + break; + case 5: + radix(a.get(), n); + break; + default: + printf("invalid algo type:%d", algo); + exit(1); + break; + } + printf("%.5f\n", std::chrono::duration(now() - beg).count()); + assert(std::is_sorted(a.get(), a.get() + n)); +} diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/4_1.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/4_1.cpp" new file mode 100644 index 000000000..d87af461c --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/4_1.cpp" @@ -0,0 +1,57 @@ +#include "util.hpp" + +struct XorShiftRNG { + u32 seed; + + f64 gen_f64() { + seed ^= seed << 13; + seed ^= seed >> 17; + seed ^= seed << 5; + return seed * (1.0 / -1u); + } + + u32 gen_u32() { + seed ^= seed << 13; + seed ^= seed >> 17; + seed ^= seed << 5; + return seed; + } +}; + +std::pair, std::unique_ptr> gen(u32 n) { + std::unique_ptr a(new f64[n * n]); + std::unique_ptr b(new u32[n]); + XorShiftRNG rng{u32(time(nullptr))}; + for (u32 i = 0; i < n; ++i) { + for (u32 j = 0; j < n; ++j) { + a[i * n + j] = rng.gen_f64(); + } + } + for (u32 i = 0; i < n; ++i) { + b[i] = rng.gen_u32() % n + 1; + } + return {std::move(a), std::move(b)}; +} + +i32 main(i32 argc, i8** argv) { + if (argc != 3) { + printf("usage: %s [n] [th]\n", argv[0]); + exit(1); + } + u32 n = atoi(argv[1]); + u32 th = atoi(argv[2]); + auto pr = gen(n); + std::unique_ptr s(new f64[n]); + auto now = std::chrono::high_resolution_clock::now; + auto beg = now(); +#pragma omp parallel for schedule(runtime) num_threads(th) + for (u32 i = 0; i < n; ++i) { + f64 x = 0; + u32 b = pr.second[i]; + for (u32 j = 0; j < b; ++j) { + x += pr.first[i * n + j]; + } + s[i] = x / b; + } + printf("%.5f\n", std::chrono::duration(now() - beg).count()); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/4_2.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/4_2.cpp" new file mode 100644 index 000000000..2805e21d8 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/4_2.cpp" @@ -0,0 +1,130 @@ +#include +#include "util.hpp" + +template +struct ThreadPool { + std::unique_ptr ths; + u32 n; + std::atomic no_more; + Mutex mu; + CondVar cv; + std::queue q; + + static void *task(void *arg) { + ThreadPool *self = (ThreadPool *)arg; + while (true) { + self->mu.lock(); + // Mesa monitor + while (self->q.empty() && !self->no_more) { + self->cv.wait(&self->mu); + } + if (self->q.empty()) { // no more + self->mu.unlock(); + break; + } + T t = self->q.front(); + self->q.pop(); + self->mu.unlock(); + t(); + } + return nullptr; + } + + ThreadPool(u32 n) : ths(new pthread_t[n]), n(n), no_more(false) { + for (u32 i = 0; i < n; ++i) { + pthread_create(&ths[i], nullptr, task, this); + } + } + + // this will wait until all tasks are finished + ~ThreadPool() { + no_more = true; + cv.notify_all(); + for (u32 i = 0; i < n; ++i) { + pthread_join(ths[i], nullptr); + } + } + + void push(T t) { + mu.lock(); + q.push(t); + cv.notify_one(); + mu.unlock(); + } +}; + +struct XorShiftRNG { + u32 seed; + + f64 gen_f64() { + seed ^= seed << 13; + seed ^= seed >> 17; + seed ^= seed << 5; + return seed * (1.0 / -1u); + } + + u32 gen_u32() { + seed ^= seed << 13; + seed ^= seed >> 17; + seed ^= seed << 5; + return seed; + } +}; + +std::pair, std::unique_ptr> gen(u32 n) { + std::unique_ptr a(new f64[n * n]); + std::unique_ptr b(new u32[n]); + XorShiftRNG rng{u32(time(nullptr))}; + for (u32 i = 0; i < n; ++i) { + for (u32 j = 0; j < n; ++j) { + a[i * n + j] = rng.gen_f64(); + } + } + for (u32 i = 0; i < n; ++i) { + b[i] = rng.gen_u32() % n + 1; + } + return {std::move(a), std::move(b)}; +} + +u32 n; +std::pair, std::unique_ptr> pr; +std::unique_ptr s; + +struct Op { + u32 beg, end; + void operator()() { + for (u32 i = beg; i < end; ++i) { + f64 x = 0; + u32 b = pr.second[i]; + for (u32 j = 0; j < b; ++j) { + x += pr.first[i * n + j]; + } + s[i] = x / b; + } + } +}; + +i32 main(i32 argc, i8 **argv) { + if (argc != 4) { + printf("usage: %s [n] [th] [blk]\n", argv[0]); + exit(1); + } + n = atoi(argv[1]); + u32 th = atoi(argv[2]); + u32 blk = atoi(argv[3]); + pr = gen(n); + s.reset(new f64[n]); + auto now = std::chrono::high_resolution_clock::now; + auto beg = now(); + { + ThreadPool pool(th); + u32 i = 0; + for (; i + blk <= n; i += blk) { + pool.push(Op{i, i + blk}); + } + if (i != n) { + pool.push(Op{i, n}); + } + } + printf("%.5f\n", std::chrono::duration(now() - beg).count()); +} \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/makefile" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/makefile" new file mode 100644 index 000000000..3f9156f03 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/makefile" @@ -0,0 +1,22 @@ +all: 1_1 1_2 2 3 4_1 4_2 + +1_1: 1_1.cpp util.hpp + mpiicc -Ofast -march=native -fopenmp -std=c++11 1_1.cpp -o 1_1 + +1_2: 1_2.cpp util.hpp + mpiicc -Ofast -march=native -fopenmp -std=c++11 1_2.cpp -o 1_2 + +2: 2.cpp util.hpp + g++ -O2 -pthread -std=c++11 2.cpp -o 2 + +3: 3.cpp util.hpp + icc -O3 -march=native -fopenmp -std=c++11 3.cpp -o 3 + +4_1: 4_1.cpp util.hpp + g++ -Ofast -march=native -fopenmp -std=c++11 4_1.cpp -o 4_1 + +4_2: 4_2.cpp util.hpp + g++ -Ofast -march=native -pthread -fopenmp -std=c++11 4_2.cpp -o 4_2 + +clean: + rm 1_1 1_2 2 3 4_1 4_2 diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_1.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_1.png" new file mode 100644 index 000000000..c009771ba Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_1.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_2.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_2.png" new file mode 100644 index 000000000..cd0005ab8 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_2.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_3.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_3.png" new file mode 100644 index 000000000..be2d63ddd Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_3.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_4.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_4.png" new file mode 100644 index 000000000..92eaa2245 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/1_4.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_1.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_1.png" new file mode 100644 index 000000000..b942f683c Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_1.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_2.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_2.png" new file mode 100644 index 000000000..4b4d6f3ed Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_2.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_3.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_3.png" new file mode 100644 index 000000000..a47cf0588 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_3.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_4.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_4.png" new file mode 100644 index 000000000..4bf53c0cf Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/3_4.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/4_1.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/4_1.png" new file mode 100644 index 000000000..905c1d195 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/4_1.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/4_2.png" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/4_2.png" new file mode 100644 index 000000000..a7b003f05 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/4_2.png" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot1.py" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot1.py" new file mode 100644 index 000000000..9b088de33 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot1.py" @@ -0,0 +1,72 @@ +import matplotlib.pyplot as plt +import numpy as np + +# data = '''0.00854s | 0.03399s | 0.13942s | 0.65621s | +# 0.00422s | 0.01728s | 0.06795s | 0.29894s | +# 0.00302s | 0.01205s | 0.04818s | 0.19448s | +# 0.00247s | 0.00908s | 0.03776s | 0.15100s | +# 0.00223s | 0.00788s | 0.03150s | 0.12579s | +# 0.00211s | 0.00728s | 0.02834s | 0.11284s | +# 0.00206s | 0.00700s | 0.02673s | 0.10631s | +# 0.00190s | 0.00699s | 0.02403s | 0.10599s | +# 0.00204s | 0.00657s | 0.02682s | 0.08540s | +# 0.00205s | 0.00668s | 0.02459s | 0.09724s | +# 0.00170s | 0.00590s | 0.02163s | 0.07809s | +# 0.00162s | 0.00520s | 0.02050s | 0.07308s | +# 0.00170s | 0.00484s | 0.01808s | 0.06524s | +# 0.00125s | 0.00459s | 0.01531s | 0.05890s | +# 0.00114s | 0.00394s | 0.01380s | 0.05368s |''' +# data = '''0.04101s | 0.16158s | 0.65047s | 2.68258s | +# 0.04319s | 0.15372s | 0.59899s | 2.36608s | +# 0.03149s | 0.10811s | 0.44896s | 1.68811s | +# 0.02759s | 0.09536s | 0.32899s | 1.31387s | +# 0.02540s | 0.09263s | 0.30892s | 1.04187s | +# 0.02378s | 0.08665s | 0.23575s | 0.89775s | +# 0.02327s | 0.08849s | 0.31522s | 0.72554s | +# 0.02297s | 0.08393s | 0.31131s | 0.68641s | +# 0.02305s | 0.08256s | 0.33165s | 0.67264s | +# 0.02351s | 0.08366s | 0.31419s | 0.67175s | +# 0.02503s | 0.08926s | 0.31526s | 1.21293s | +# 0.02523s | 0.08426s | 0.32665s | 1.21559s | +# 0.02560s | 0.08943s | 0.31148s | 1.25053s | +# 0.02554s | 0.08869s | 0.32243s | 1.23589s | +# 0.02619s | 0.08890s | 0.31714s | 1.20727s |''' +# data = '''0.00998s | 0.03824s | 0.16265s | 0.71438s | +# 0.00298s | 0.01066s | 0.04209s | 0.17797s | +# 0.00221s | 0.00725s | 0.02752s | 0.09978s | +# 0.00192s | 0.00556s | 0.02050s | 0.07388s |''' +data = '''0.04069s | 0.16151s | 0.64825s | 2.66666s | +0.02576s | 0.08296s | 0.31809s | 1.24863s | +0.02205s | 0.04532s | 0.16570s | 0.63102s | +0.02364s | 0.08014s | 0.14077s | 0.52902s |''' + +data = list(map(float, data.split('s |')[:-1])) +# th = (1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24) +th = (1, 4, 9, 16) +n3600 = data[0] / np.array(data[0::4]) +n7200 = data[1] / np.array(data[1::4]) +n14400 = data[2] / np.array(data[2::4]) +n28800 = data[3] / np.array(data[3::4]) + +plt.xlabel('thread') +plt.ylabel('accelerate factor') + +plt.plot(th, n3600, label='n=3600', color='red', marker='x') +plt.plot(th, n7200, label='n=7200', color='blue', marker='x') +plt.plot(th, n14400, label='n=14400', color='green', marker='x') +plt.plot(th, n28800, label='n=28800', color='yellow', marker='x') +plt.plot(th, th, label='y=x', color='purple', marker='x') + +for a, b in zip(th, n3600): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, n7200): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, n14400): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, n28800): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, th): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') + +plt.legend() +plt.show() diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot3_1.py" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot3_1.py" new file mode 100644 index 000000000..10c933e0b --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot3_1.py" @@ -0,0 +1,54 @@ +import matplotlib.pyplot as plt +import numpy as np + +data = '''0.02705s | 0.10475s | 0.30914s | 0.65088s | 1.15128s | +0.01180s | 0.01668s | 0.02386s | 0.04726s | 0.06355s | +0.00185s | 0.00379s | 0.00592s | 0.00790s | 0.00996s | +0.00107s | 0.00222s | 0.00347s | 0.00468s | 0.00608s | +0.00017s | 0.00039s | 0.00065s | 0.00089s | 0.00111s |''' +# data = '''0.02705s | 0.10475s | 0.30914s | 0.65088s | 1.15128s | +# 0.01180s | 0.01668s | 0.02386s | 0.04726s | 0.06355s | +# 0.00218s | 0.00279s | 0.00348s | 0.00447s | 0.00714s | +# 0.00185s | 0.00379s | 0.00592s | 0.00790s | 0.00996s | +# 0.00107s | 0.00222s | 0.00347s | 0.00468s | 0.00608s | +# 0.00017s | 0.00039s | 0.00065s | 0.00089s | 0.00111s |''' + + +data = list(map(float, data.split('s |')[:-1])) + +ns = (1e4, 4e4, 6e4, 8e4, 1e5) +count = data[0:5] +count_par = data[5:10] +# count_llx = data[10:15] +# qsort = data[15:20] +# std_sort = data[20:25] +# radix = data[25:30] +qsort = data[10:15] +std_sort = data[15:20] +radix = data[20:25] + +plt.xlabel('n') +plt.ylabel('time') + +plt.semilogy(ns, count, label='count', color='red', marker='x') +plt.semilogy(ns, count_par, label='count_par', color='blue', marker='x') +# plt.semilogy(ns, count_llx, label='count_llx', color='purple', marker='x') +plt.semilogy(ns, qsort, label='qsort', color='green', marker='x') +plt.semilogy(ns, std_sort, label='std_sort', color='yellow', marker='x') +plt.semilogy(ns, radix, label='radix', color='cyan', marker='x') + +for a, b in zip(ns, count): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(ns, count_par): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +# for a, b in zip(ns, count_llx): +# plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(ns, qsort): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(ns, std_sort): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(ns, radix): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') + +plt.legend() +plt.show() diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot3_2.py" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot3_2.py" new file mode 100644 index 000000000..f738f1f4e --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot3_2.py" @@ -0,0 +1,64 @@ +import matplotlib.pyplot as plt +import numpy as np + +# data = ''' 0.02905s | 0.10707s | 0.27759s | 0.62348s | 1.14682s | +# 0.02637s | 0.05793s | 0.12984s | 0.33981s | 0.54269s | +# 0.01441s | 0.03087s | 0.08311s | 0.18087s | 0.29576s | +# 0.01119s | 0.02185s | 0.05794s | 0.12896s | 0.21299s | +# 0.01019s | 0.02480s | 0.04664s | 0.09645s | 0.15156s | +# 0.00961s | 0.01864s | 0.05524s | 0.07933s | 0.16566s | +# 0.00954s | 0.02151s | 0.05835s | 0.09079s | 0.11473s | +# 0.00915s | 0.01995s | 0.04034s | 0.08840s | 0.12367s | +# 0.00854s | 0.01309s | 0.03502s | 0.07951s | 0.10757s | +# 0.00843s | 0.01373s | 0.03304s | 0.07125s | 0.10845s | +# 0.00791s | 0.01489s | 0.03750s | 0.05872s | 0.10018s | +# 0.00915s | 0.01260s | 0.02596s | 0.05035s | 0.07533s | +# 0.00831s | 0.01325s | 0.03003s | 0.04433s | 0.07256s |''' +data = '''0.02956s | 0.11585s | 0.30878s | 0.66505s | 1.14823s | +0.01412s | 0.02988s | 0.07663s | 0.13395s | 0.23511s | +0.00409s | 0.01290s | 0.02252s | 0.03510s | 0.04931s | +0.00259s | 0.00696s | 0.01416s | 0.01626s | 0.03212s | +0.00191s | 0.00444s | 0.00784s | 0.01264s | 0.01446s | +0.00139s | 0.00346s | 0.00547s | 0.00776s | 0.01606s | +0.00123s | 0.00261s | 0.00444s | 0.00669s | 0.00908s | +0.00123s | 0.00253s | 0.00384s | 0.00541s | 0.00710s | +0.00098s | 0.00229s | 0.00341s | 0.00468s | 0.00553s | +0.00098s | 0.00218s | 0.00338s | 0.00467s | 0.00602s | +0.00099s | 0.00194s | 0.00284s | 0.00405s | 0.00563s | +0.00093s | 0.00173s | 0.00288s | 0.00388s | 0.00510s | +0.00104s | 0.00172s | 0.00258s | 0.00380s | 0.00516s | ''' + +data = list(map(float, data.split('s |')[:-1])) + +th = (1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24) +n2e4 = data[0] / np.array(data[0::5]) +n4e4 = data[1] / np.array(data[1::5]) +n6e4 = data[2] / np.array(data[2::5]) +n8e4 = data[3] / np.array(data[3::5]) +n1e5 = data[4] / np.array(data[4::5]) + +plt.xlabel('thread') +plt.ylabel('accelerate factor') + +plt.plot(th, n2e4, label='n=2e4', color='red', marker='x') +plt.plot(th, n4e4, label='n=4e4', color='blue', marker='x') +plt.plot(th, n6e4, label='n=6e4', color='green', marker='x') +plt.plot(th, n8e4, label='n=8e4', color='yellow', marker='x') +plt.plot(th, n1e5, label='n=1e5', color='cyan', marker='x') +plt.plot(th, th, label='y=x', color='purple', marker='x') + +for a, b in zip(th, n2e4): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, n4e4): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, n6e4): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, n8e4): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, n1e5): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(th, th): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') + +plt.legend() +plt.show() diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot4.py" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot4.py" new file mode 100644 index 000000000..c7cbfcfff --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/plot4.py" @@ -0,0 +1,39 @@ +import matplotlib.pyplot as plt + +blks = [1] + list(range(10, 90 + 10, 10)) + list(range(100, 2000 + 100, 100)) +static = '''0.274s | 0.287s | 0.284s | 0.286s | 0.285s | 0.286s | 0.288s | 0.286s | 0.283s | 0.287s | +0.286s | 0.285s | 0.287s | 0.283s | 0.286s | 0.290s | 0.281s | 0.290s | 0.294s | 0.286s | +0.281s | 0.290s | 0.287s | 0.294s | 0.300s | 0.301s | 0.304s | 0.308s | 0.304s | 0.285s |''' +dynamic = '''0.174s | 0.161s | 0.160s | 0.160s | 0.160s | 0.160s | 0.161s | 0.162s | 0.161s | 0.161s | +0.161s | 0.161s | 0.172s | 0.164s | 0.180s | 0.194s | 0.190s | 0.168s | 0.175s | 0.182s | +0.183s | 0.198s | 0.212s | 0.220s | 0.227s | 0.236s | 0.249s | 0.261s | 0.273s | 0.283s | ''' +guided = '''0.238s | 0.223s | 0.232s | 0.234s | 0.228s | 0.232s | 0.231s | 0.227s | 0.231s | 0.232s | +0.228s | 0.223s | 0.230s | 0.226s | 0.233s | 0.228s | 0.224s | 0.235s | 0.233s | 0.234s | +0.227s | 0.230s | 0.238s | 0.239s | 0.242s | 0.249s | 0.254s | 0.265s | 0.271s | 0.284s |''' +pthread = '''0.175s | 0.161s | 0.159s | 0.159s | 0.165s | 0.165s | 0.168s | 0.163s | 0.163s | 0.165s | +0.167s | 0.168s | 0.173s | 0.177s | 0.182s | 0.192s | 0.195s | 0.185s | 0.169s | 0.180s | +0.189s | 0.193s | 0.215s | 0.223s | 0.232s | 0.245s | 0.253s | 0.265s | 0.276s | 0.284s |''' + +static = list(map(float, static.split('s |')[:-1])) +dynamic = list(map(float, dynamic.split('s |')[:-1])) +guided = list(map(float, guided.split('s |')[:-1])) +pthread = list(map(float, pthread.split('s |')[:-1])) + +plt.xlabel('blk') +plt.ylabel('time') +# plt.plot(blks, static, label='static', color='red', marker='x') +plt.plot(blks, dynamic, label='dynamic', color='blue', marker='x') +# plt.plot(blks, guided, label='guided', color='green', marker='x') +plt.plot(blks, pthread, label='pthread', color='red', marker='x') + +# for a, b in zip(blks, static): +# plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(blks, dynamic): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +# for a, b in zip(blks, guided): +# plt.text(a, b, '%.3f' % b, ha='center', va='bottom') +for a, b in zip(blks, pthread): + plt.text(a, b, '%.3f' % b, ha='center', va='bottom') + +plt.legend() +plt.show() diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/run.py" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/run.py" new file mode 100644 index 000000000..fdd7d604d --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/pic/run.py" @@ -0,0 +1,101 @@ +from __future__ import print_function +import subprocess +import re + +# for blk in [1] + range(10, 100, 10): +# m = 1e9 +# for _ in range(5): +# o = subprocess.check_output(['./4_2', '48000', '24', str(blk)]) +# m = min(m, float(o)) +# print('%.3f' % m, end='s | ') +# print() +# for i in range(2): +# for j in range(1, 11): +# blk = i * 1000 + j * 100 +# m = 1e9 +# for _ in range(5): +# o = subprocess.check_output(['./4_2', '48000', '24', str(blk)]) +# m = min(m, float(o)) +# print('%.3f' % m, end='s | ') +# print() + +# for blk in range(100, 4000 + 100, 100): +# m = 1e9 +# for _ in range(5): +# o = subprocess.check_output( +# ['./4_1', '96000', '24'], env={'OMP_SCHEDULE': 'static,%d' % blk}) +# m = min(m, float(o)) +# pass +# print('%.5f' % m, end='s | ') +# print() +t1s = [] +t2s = [] + + +# for i in (1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24): +# t1s.append('') +# t2s.append('') +# for n in (3600, 7200, 14400, 28800): +# t1 = 1e9 +# t2 = 1e9 +# for _ in range(5): +# o = subprocess.check_output( +# ['srun', '-n', str(i), './1_1', str(n)]) +# t1 = min(t1, float(re.findall( +# r'prod\(\+MPI_Scatter\): ([^s]+)s', o)[0])) +# t2 = min(t1, float(re.findall( +# r'prod\(-MPI_Scatter\): ([^s]+)s', o)[0])) +# t1s[-1] += '%.5fs | ' % t1 +# t2s[-1] += '%.5fs | ' % t2 + +# for i in (1, 4, 9, 16): +# t1s.append('') +# t2s.append('') +# for n in (3600, 7200, 14400, 28800): +# t1 = 1e9 +# t2 = 1e9 +# for _ in range(5): +# o = subprocess.check_output( +# ['srun', '-n', str(i), './1_2', str(n)]) +# t1 = min(t1, float(re.findall( +# r'prod\(\+MPI_Scatter\): ([^s]+)s', o)[0])) +# t2 = min(t1, float(re.findall( +# r'prod\(-MPI_Scatter\): ([^s]+)s', o)[0])) +# t1s[-1] += '%.5fs | ' % t1 +# t2s[-1] += '%.5fs | ' % t2 + +# for t1 in t1s: +# print(t1) +# print() + +# for t2 in t2s: +# print(t2) +# print() + +# for algo in range(0, 6): +# for n in range(20000, 100000+20000, 20000): +# m = 1e9 +# for _ in range(0, 5): +# o = subprocess.check_output(['./3', str(algo), str(n)]) +# m = min(m, float(o)) +# print('%.5f' % m, end='s | ') +# print() + +# for th in [1] + list(range(2, 24 + 2, 2)): +# for n in range(20000, 100000+20000, 20000): +# m = 1e9 +# for _ in range(0, 5): +# o = subprocess.check_output(['./3', '1', str(n), str(th)]) +# m = min(m, float(o)) +# print('%.5f' % m, end='s | ') +# print() + +print() +for th in [1] + list(range(2, 24 + 2, 2)): + for n in range(20000, 100000+20000, 20000): + m = 1e9 + for _ in range(0, 5): + o = subprocess.check_output(['./3', '2', str(n), str(th)]) + m = min(m, float(o)) + print('%.5f' % m, end='s | ') + print() diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/prod.cpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/prod.cpp" new file mode 100644 index 000000000..df1337d1c --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/prod.cpp" @@ -0,0 +1,82 @@ +#include +#include +#include +#include + +using f64 = double; +using u32 = unsigned; +using f64x4 = __m256d; + +f64 dot(f64 *a, f64 *b, int n) { + f64x4 ymm0 = _mm256_setzero_pd(); + f64x4 ymm1 = _mm256_setzero_pd(); + f64x4 ymm2 = _mm256_setzero_pd(); + f64x4 ymm3 = _mm256_setzero_pd(); + f64x4 ymm4, ymm5, ymm6, ymm7; + for (int i = 0; i < n; i += 32) { + ymm4 = _mm256_loadu_pd(a + i + 0); + ymm5 = _mm256_loadu_pd(a + i + 4); + ymm6 = _mm256_loadu_pd(a + i + 8); + ymm7 = _mm256_loadu_pd(a + i + 12); + ymm4 = _mm256_fmadd_pd(ymm4, _mm256_loadu_pd(b + i + 0), ymm0); + ymm5 = _mm256_fmadd_pd(ymm5, _mm256_loadu_pd(b + i + 4), ymm1); + ymm6 = _mm256_fmadd_pd(ymm6, _mm256_loadu_pd(b + i + 8), ymm2); + ymm7 = _mm256_fmadd_pd(ymm7, _mm256_loadu_pd(b + i + 12), ymm3); + ymm0 = _mm256_loadu_pd(a + i + 16); + ymm1 = _mm256_loadu_pd(a + i + 20); + ymm2 = _mm256_loadu_pd(a + i + 24); + ymm3 = _mm256_loadu_pd(a + i + 28); + ymm0 = _mm256_fmadd_pd(ymm0, _mm256_loadu_pd(b + i + 16), ymm4); + ymm1 = _mm256_fmadd_pd(ymm1, _mm256_loadu_pd(b + i + 20), ymm5); + ymm2 = _mm256_fmadd_pd(ymm2, _mm256_loadu_pd(b + i + 24), ymm6); + ymm3 = _mm256_fmadd_pd(ymm3, _mm256_loadu_pd(b + i + 28), ymm7); + } + ymm0 = _mm256_add_pd(ymm0, ymm2); + ymm1 = _mm256_add_pd(ymm1, ymm3); + ymm0 = _mm256_add_pd(ymm0, ymm1); + f64 sum[4]; + _mm256_store_pd(sum, ymm0); + return sum[0] + sum[1] + sum[2] + sum[3]; +} + +std::unique_ptr prod(const std::unique_ptr &mat, const std::unique_ptr &vec, int n) { + std::unique_ptr ret(new double[n]); + for (int i = 0; i < n; ++i) { + ret[i] = dot(&mat[i * n], &vec[0], n); + } + return ret; +} + +struct XorShiftRNG { + u32 seed; + + XorShiftRNG(u32 seed) : seed(seed ? seed : 1) {} + + f64 gen() { + seed ^= seed << 13; + seed ^= seed >> 17; + seed ^= seed << 5; + return seed * (1.0 / -1u); + } +}; + +int main(int argc, char **argv) { + int N = atoi(argv[1]); + std::unique_ptr mat(new double[N * N]); + std::unique_ptr vec(new double[N]); + + XorShiftRNG rng{19260817}; + for (int i = 0; i < N; ++i) { + for (int j = 0; j < N; ++j) { + mat[i * N + j] = rng.gen(); + } + } + for (int i = 0; i < N; ++i) { + vec[i] = rng.gen(); + } + using namespace std::chrono; + auto now = high_resolution_clock::now; + auto beg = now(); + auto c = prod(mat, vec, N); + printf("%f\n", duration(now() - beg).count()); +} diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/report.md" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/report.md" new file mode 100644 index 000000000..a5bbdad3b --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/report.md" @@ -0,0 +1,487 @@ +--- +title: "高性能计算导论第五次作业" +author: "李晨昊 2017011466" +date: "2019-6-1" +output: + pdf_document: + latex_engine: xelatex + number_sections: yes + toc: yes + word_document: + toc: yes + html_document: + toc: yes +header-includes: \usepackage{ctex} \usepackage{graphicx} +--- + +\newpage + +# Programming Assignment +第三次的作业中我就已经测量过了性能数据,不过这次换了一下编译器和编译选项,所以还是全部重测了一遍。 + +既然本次对于性能做出了明确的要求,这里理应采用-Ofast -march=native优化,让编译器能够自由地用SIMD指令来优化浮点运算,不必担心浮点误差对结果的影响。不过经测试性能提升很小,于是我又尝试手写了一些循环展开和SIMD,发现提升也很小,我认为这可能有下面几方面的原因: + +1. 即使我没有写这些东西,编译器也已经优化的比较好了;测试表明,运行速度gcc -Ofast -march=native $\approx$ icc -Ofast -march=native $\approx$ icc -O3 $\approx$ 手写SIMD $>$ gcc -O3,不过我对此表示比较怀疑,因为这种优化往往能改变程序的运行结果,编译器应该不能擅自做这么激进的优化 + +2. SIMD优化本身对于双精度浮点数的效果就没有对于单精度浮点数大。由于二者数据宽度的差别,理想状态下单精度浮点数的运算效率二倍于双精度浮点数。但是既然题目要求使用双精度浮点数,我对此也没有什么办法 + +3. 任务本身瓶颈在访存,而非运算。矩阵乘向量过程运算简单但访存需求大,而这种简单浮点运算的速度其实是比访存快得多的。例如在intel手册上可以查到,在服务器的Haswell架构的CPU上,`_mm256_fmadd_pd`的latency为6,CPI为0.5,这比访存大概快了两个数量级 + +由于最后速度提升并不明显,我在最终版本里没有使用手写的SIMD,但是还是把它留在了prod.cpp中。这代码基本上是我的机器上的clang生成的汇编代码反向翻译回C++的,在我的机器上的效果也还算不错,明显地快于clang和gcc的-O2生成的代码。但是因为我没有icc(而且icc据说生成的代码在AMD平台上效率并不高),所以没办法在我的平台上和icc的-O2进行比较。 + +按列划分,不考虑数据分发和收集,耗时结果如下: + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +|---|--------|--------|---------|---------| +| 1 | 0.00854s | 0.03399s | 0.13942s | 0.65621s | +| 2 | 0.00422s | 0.01728s | 0.06795s | 0.29894s | +| 3 | 0.00302s | 0.01205s | 0.04818s | 0.19448s | +| 4 | 0.00247s | 0.00908s | 0.03776s | 0.15100s | +| 5 | 0.00223s | 0.00788s | 0.03150s | 0.12579s | +| 6 | 0.00211s | 0.00728s | 0.02834s | 0.11284s | +| 8 | 0.00206s | 0.00700s | 0.02673s | 0.10631s | +| 9 | 0.00190s | 0.00699s | 0.02403s | 0.10599s | +| 10 | 0.00204s | 0.00657s | 0.02682s | 0.08540s | +| 12 | 0.00205s | 0.00668s | 0.02459s | 0.09724s | +| 15 | 0.00170s | 0.00590s | 0.02163s | 0.07809s | +| 16 | 0.00162s | 0.00520s | 0.02050s | 0.07308s | +| 18 | 0.00170s | 0.00484s | 0.01808s | 0.06524s | +| 20 | 0.00125s | 0.00459s | 0.01531s | 0.05890s | +| 24 | 0.00114s | 0.00394s | 0.01380s | 0.05368s | + +加速比曲线如下: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/1_1.png} +\end{center} + +考虑数据分发和收集,耗时结果如下: + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +|---|--------|--------|---------|---------| +| 1 | 0.04101s | 0.16158s | 0.65047s | 2.68258s | +| 2 | 0.04319s | 0.15372s | 0.59899s | 2.36608s | +| 3 | 0.03149s | 0.10811s | 0.44896s | 1.68811s | +| 4 | 0.02759s | 0.09536s | 0.32899s | 1.31387s | +| 5 | 0.02540s | 0.09263s | 0.30892s | 1.04187s | +| 6 | 0.02378s | 0.08665s | 0.23575s | 0.89775s | +| 8 | 0.02327s | 0.08849s | 0.31522s | 0.72554s | +| 9 | 0.02297s | 0.08393s | 0.31131s | 0.68641s | +| 10 | 0.02305s | 0.08256s | 0.33165s | 0.67264s | +| 12 | 0.02351s | 0.08366s | 0.31419s | 0.67175s | +| 15 | 0.02503s | 0.08926s | 0.31526s | 1.21293s | +| 16 | 0.02523s | 0.08426s | 0.32665s | 1.21559s | +| 18 | 0.02560s | 0.08943s | 0.31148s | 1.25053s | +| 20 | 0.02554s | 0.08869s | 0.32243s | 1.23589s | +| 24 | 0.02619s | 0.08890s | 0.31714s | 1.20727s | + +加速比曲线如下: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/1_2.png} +\end{center} + +对于每个n,计算误差二范数的计算量都是一样的(而且都很小),比较不同线程数下的执行速度没什么价值,直接用默认的线程数即可: + +| n=3600 | n=7200 | n=14400 | n=28800 | +|--------|--------|---------|---------| +| 0.00744s | 0.00416s | 0.00393s | 0.00478s | + +这些数据都是多次测量取得的最小值,而测量过程中也发现时间的波动很大。我认为原因在于数据量太小,以至于计算耗时远小于开启线程和同步线程的额外开销。事实上,我把它们替换回了串行版本后,再测试一次,结果为: + +| n=3600 | n=7200 | n=14400 | n=28800 | +|--------|--------|---------|---------| +| 0.00000s | 0.00001s | 0.00003s | 0.00006s | + +这算是比较正常的结果,而上面的测试结果里的波动也就可以解释了:几乎所有的时间都属于$T_{overhead}$。对于数据规模如此小的计算,没有使用并行的意义。 + +按子矩阵划分,不考虑数据分发和收集,耗时结果如下: + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +| -- | ------ | ------ | ------- | ------- | +| 1 | 0.00998s | 0.03824s | 0.16265s | 0.71438s | +| 4 | 0.00298s | 0.01066s | 0.04209s | 0.17797s | +| 9 | 0.00221s | 0.00725s | 0.02752s | 0.09978s | +| 16 | 0.00192s | 0.00556s | 0.02050s | 0.07388s | + +加速比曲线如下: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/1_3.png} +\end{center} + +考虑数据分发和收集,耗时结果如下: + +| p | n=3600 | n=7200 | n=14400 | n=28800 | +| -- | ------ | ------ | ------- | ------- | +| 1 | 0.04069s | 0.16151s | 0.64825s | 2.66666s | +| 4 | 0.02576s | 0.08296s | 0.31809s | 1.24863s | +| 9 | 0.02205s | 0.04532s | 0.16570s | 0.63102s | +| 16 | 0.02364s | 0.08014s | 0.14077s | 0.52902s | + +加速比曲线如下: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/1_4.png} +\end{center} + +由于n没有变化,这里就不重复列出计算二范数的耗时了。 + +# Programming Assignment 4.5 +题目要求大概让我们实现一个线程池。至于那个给我们的主函数,我没有管它。 + +线程池的大致实现思路是,开启多个线程,它们一直都在等待任务而不退出,这样在创建任务时就避免了创建线程的开销。当然开销肯定还是存在的,这是线程间通信的开销,比创建线程要小的多。 + +"等待任务"的具体实现方法是使用条件变量,维护一个任务队列,当队列为空时需要等待,当加入新任务时可唤醒一个线程。此外,线程池还需要一个退出的逻辑,因此等待条件中还需要检查是否退出,如果已经退出,且任务队列为空,则本线程的执行结束。 + +关键代码如下: + +```cpp + static void *task(void *arg) { + ThreadPool *self = (ThreadPool *)arg; + while (true) { + self->mu.lock(); // 管程开始 + // Mesa 管程,用while + while (self->q.empty() && !self->no_more) { + self->cv.wait(&self->mu); + } + if (self->q.empty()) { // 证明self->no_more成立,可以退出 + self->mu.unlock(); + break; + } + T t = self->q.front(); + self->q.pop(); + self->mu.unlock(); // 管程结束 + t(); // 管程外执行任务 + } + return nullptr; + } +``` + +题目中要求的:"each time it generates a new block of tasks, it awakens a thread with a condition signal",对应于代码: + +```cpp + void push(T t) { + mu.lock(); + q.push(t); + cv.notify_one(); + mu.unlock(); + } +``` + +题目中要求的:"When the main thread completes generating tasks, it sets a global variable indicating that there will be no more tasks, and awakens all the threads with a condition broadcast",对应于代码: + +```cpp + ~ThreadPool() { + no_more = true; + cv.notify_all(); + ... + } +``` + +我对锁和条件变量等做了一点封装,这样自己用起来舒服一些。当然stl其实已经封装好了pthread,不过为了练习的目的,显然我们不应该使用使用stl的版本。不过这个no_more变量,我使用了std::atomic\,至于为什么不应该用普通bool或者volatile bool,我在上次报告中说的比较清楚了。也许它们生成的汇编是一模一样的,但是我们期望的语义是atomic,不是volatile。比较可惜pthread没有提供原子变量,因此这里是不得不用到stl的。 + +关于链表,直接用一个锁和一个std::list来实现。这里没有维护链表的有序性,因为题目中并没有说需要维护,它只是一个普通的链表而已。对于插入操作,直接使用了`push_back`将它插入到链表头。如果真的想要一个课件中的那种维护了有序性和元素唯一性的链表,毫无疑问应该使用平衡树,例如std::set来实现,使用链表是完全没有道理的。 + +# Programming Assignment 5.3 +## 普通版本 +这个小标题里回答了题目提出的几个问题,并且按照最直接的想法并行化了计数排序。 + +1. 如果我们试图并行化i的for循环(外层循环),哪些变量应该是私有的,哪些变量应该是共享的? + +私有:i,j,count + +共享:a,n,temp + +2. 如果我们使用前面定义的作用域来并行化i的for循环,是否存在循环依赖?请解释你的回答。 + +不存在。注意到循环体内没有对于共享可变变量的读操作,因此不存在循环依赖。 + +3. 我们是否能够并行化对memcpy的调用?我们是否能够修改代码,使得这部分代码可以并行化? + +可以,分段并行memcpy即可。但是这样很有可能并不能达到较好的效果,因为memcpy本身非常快,很有可能耗时无法弥补开启线程等操作的耗时。 + +4. 编写一个包含对计数排序程序并行化实现的C程序。 + +具体实现见代码。 + +这里对于内层循环做了一点优化,用两个循环代替了原来的一个: + +```cpp + u32 cnt = 0; + for (u32 j = 0; j < i; ++j) { + cnt += a[j] <= a[i]; + } + for (u32 j = i; j < n; ++j) { + cnt += a[j] < a[i]; + } +``` + +容易看出这和题目中给出的循环是等价的,它的好处在于大大减少了循环体内所需的判断。不过,其实这种优化应该算是比较平凡的,我很遗憾编译器没有自己看出来,还是需要人的智慧的帮助。 + +同时注意到这个循环可以很好的用SIMD指令来实现。但是,是否能够进行SIMD优化,以及优化能到什么程度,其实都完全取决于编译器的行为:我在我的机器(AMD R7-2700)上分别使用gcc 8.3.0和clang 8.0.0,在服务器上分别使用了gcc 4.8.5和icc 18.0.0(**注:icc编译出来的带openmp代码似乎只能在bootstraper上运行,在cn0xx上运行提示缺少动态链接库libiomp5.so**),对于$n=10^5$记录单线程版本耗时,记录如下: + +| | -O2 | -O3 -march=native | +|----------------------|--|--| +| AMD R7-2700, gcc 8.3.0 | 5.52603s | 0.87211s | +| AMD R7-2700, clang 8.0.0 | 0.78120s | 0.62353s | +| Intel Xeon E5-2680 v3, gcc 4.8.5 | 6.26367s | 3.89370s | +| Intel Xeon E5-2680 v3, icc 18.0.0 | 1.60733s | 1.12694s | + +不得不说clang的SIMD优化还是相当值得信赖的,符合我对它一贯的认知。$O(n^2)$能过$n=10^5$(过:指耗时小于1s),还是挺不错的。 + +然而,继续研究这几个编译器之间的差别意义不大,而循环展开/向量化这些本来应该是编译器自动完成的机械化的工作,如果让人来做意义也不太大(当然,的确在很多真实的应用场景中这种手工的优化仍然是有必要的,但是如果这里为了一个非常低效的算法去优化常数,实在没什么价值),因此我也没有继续优化下去了。 + +顺便说一下,希望贵课的服务器能早日更新编译器,一个现代的编译器对于代码速度的提升可能远远超过程序员的绞尽脑汁想出来的一些奇技淫巧。 + +5. 与串行化的计数排序程序相比,并行化计数排序的性能如何?与串行化的库函数qsort相比,并行化的计数排序性能如何? + +与串行化的计数排序程序相比,并行化计数排序的性能的提升很显著。 + +然而,这点常数的减小相比于它$O(n^2)$的复杂度来说,实在没什么价值,因此它显著地慢于$O(n\log n)$qsort和std::sort。而相比于真正高效的排序算法:$O(n\log U)$的基数排序来讲,它的速度更加显得可笑了。 + +测试结果如下: + +| | $n=2*10^4$ | $n=4*10^4$ | $n=6*10^4$ | $n=8*10^4$ | $n=10^5$ | +|--|------------|------------|------------|------------|----------| +| count | 0.02705s | 0.10475s | 0.30914s | 0.65088s | 1.15128s | +| count_par($th = 24$) | 0.01180s | 0.01668s | 0.02386s | 0.04726s | 0.06355s | +| qsort | 0.00185s | 0.00379s | 0.00592s | 0.00790s | 0.00996s | +| std::sort | 0.00107s | 0.00222s | 0.00347s | 0.00468s | 0.00608s | +| radix | 0.00017s | 0.00039s | 0.00065s | 0.00089s | 0.00111s | + +曲线图如下,注意时间是对数坐标的: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/3_1.png} +\end{center} + +单独测试计数排序,而选择不同的线程数,测试结果如下: + +| | $n=2*10^4$ | $n=4*10^4$ | $n=6*10^4$ | $n=8*10^4$ | $n=10^5$ | +|--|------------|------------|------------|------------|----------| +| $th = 1$ | 0.02905s | 0.10707s | 0.27759s | 0.62348s | 1.14682s | +| $th = 2$ | 0.02637s | 0.05793s | 0.12984s | 0.33981s | 0.54269s | +| $th = 4$ | 0.01441s | 0.03087s | 0.08311s | 0.18087s | 0.29576s | +| $th = 6$ | 0.01119s | 0.02185s | 0.05794s | 0.12896s | 0.21299s | +| $th = 8$ | 0.01019s | 0.02480s | 0.04664s | 0.09645s | 0.15156s | +| $th = 10$ | 0.00961s | 0.01864s | 0.05524s | 0.07933s | 0.16566s | +| $th = 12$ | 0.00954s | 0.02151s | 0.05835s | 0.09079s | 0.11473s | +| $th = 14$ | 0.00915s | 0.01995s | 0.04034s | 0.08840s | 0.12367s | +| $th = 16$ | 0.00854s | 0.01309s | 0.03502s | 0.07951s | 0.10757s | +| $th = 18$ | 0.00843s | 0.01373s | 0.03304s | 0.07125s | 0.10845s | +| $th = 20$ | 0.00791s | 0.01489s | 0.03750s | 0.05872s | 0.10018s | +| $th = 22$ | 0.00915s | 0.01260s | 0.02596s | 0.05035s | 0.07533s | +| $th = 24$ | 0.00831s | 0.01325s | 0.03003s | 0.04433s | 0.07256s | + +加速比曲线图如下: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/3_2.png} +\end{center} + +可见$n$较大时,并行效率还是比较优秀的。不过再次申明,本身这个算法是非常低效的,因此这个效果也没什么价值。如果能(高效地)并行化基数排序,那么意义将大得多,不过一眼看上去,很明显难度也大得多。 + +## 一个更好的版本 +llx同学在微信群中提出可以让每个线程负责排序一段,然后将所有排好序的序列归并起来。其实这个方法并不是很不平凡,我之前之所以没有这么做,只是因为我认为这样没有意义:既然是让每个线程负责排序一段,为什么每一段要用这个低效的算法?完全可以每个线程都使用基数排序或者快速排序,再把结果归并起来,效率一定会更高。这样做的唯一意义在于,它强行保留了题目提供的这个低效的算法,表面上看起来更加符合题目的要求,然而实际想一下,凭什么说这个算法还算是计数排序?很明显归并才是它的灵魂所在。 + +好的,上面抱怨完了。简单讲下我的实现思路:为灵活性我使用了pthread来实现这个算法。每个线程先是负责排序自己的部分,然后开始树形规约,规约代码如下: + +```cpp + pthread_barrier_wait(&calc); + for (u32 i = 0; tid + (1 << i) < th && !(tid >> i & 1); ++i) { + u32 right_tid = tid + (1 << (i + 1)); + u32 right = (right_tid >= th) ? n : each * right_tid; + merge(a_off, a_off + (each << i), a + right, aux_off); + pthread_barrier_wait(&reduce[i]); + } +``` + +设置路障个数的代码如下: + +```cpp + u32 reduce_idx = 0; + for (u32 i = th; i != 1; ++reduce_idx, i = (i + 1) >> 1) { + pthread_barrier_init(&detail::reduce[reduce_idx], nullptr, i >> 1); + } +``` + +如果助教还记得的话,会发现这代码和我上次提交的代码中的利用树形规约来加速直方图计算的代码基本是一样的,只是规约的地方需要额外处理一下边界情况。 + +顺便说一句,这个归并函数是我以前雕琢了很久得出的一个版本: + +```cpp +void merge(i32 *first, i32 *mid, i32 *last, i32 *aux) { + i32 *pos1 = first, *pos2 = mid, *pos = aux; + while (pos1 != mid && pos2 != last) { + i32 a = *pos1, b = *pos2; + i32 cmp = b < a; + *pos++ = cmp ? b : a; + pos1 += cmp ^ 1, pos2 += cmp; + } + memcpy(pos, pos1, sizeof(i32) * (mid - pos1)), pos += mid - pos1; + memcpy(first, aux, sizeof(i32) * (pos - aux)); +} +``` + +其核心在于循环内部的判断可以用条件传送指令来实现,有些人似乎觉得只要用了三元运算符就一定比if高效,这种想法是错误的,譬如说写: + +```cpp + *pos++ = *pos2 < *pos1 ? *pos2++ : *pos1++; +``` + +和写 + +```cpp + if (*pos2 < *pos1) { + *pos++ = *pos2++; + } else { + *pos++ = *pos1++; + } +``` + +经测试在最新版的gcc和clang上生成的汇编都是一模一样的,都使用了分支指令,而我的版本则真的会被用条件传送指令来实现。当然这三种写法实际上都是等价的,只是目前的编译器还都没有能力看出来而已。 + +测试结果如下,其它方法的数据都与上一节的相同: + +| | $n=2*10^4$ | $n=4*10^4$ | $n=6*10^4$ | $n=8*10^4$ | $n=10^5$ | +|--|------------|------------|------------|------------|----------| +| count | 0.02705s | 0.10475s | 0.30914s | 0.65088s | 1.15128s | +| count_par($th = 24$) | 0.01180s | 0.01668s | 0.02386s | 0.04726s | 0.06355s | +| count_llx($th = 24$) | 0.00218s | 0.00279s | 0.00348s | 0.00447s | 0.00714s | +| qsort | 0.00185s | 0.00379s | 0.00592s | 0.00790s | 0.00996s | +| std::sort | 0.00107s | 0.00222s | 0.00347s | 0.00468s | 0.00608s | +| radix | 0.00017s | 0.00039s | 0.00065s | 0.00089s | 0.00111s | + +曲线图如下,时间仍是对数坐标的: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/3_3.png} +\end{center} + +使用这个方法,选择不同的线程数,测试结果如下: + +| | $n=2*10^4$ | $n=4*10^4$ | $n=6*10^4$ | $n=8*10^4$ | $n=10^5$ | +|--|------------|------------|------------|------------|----------| +| $th = 1$ | 0.02956s | 0.11585s | 0.30878s | 0.66505s | 1.14823s | +| $th = 2$ | 0.01412s | 0.02988s | 0.07663s | 0.13395s | 0.23511s | +| $th = 4$ | 0.00409s | 0.01290s | 0.02252s | 0.03510s | 0.04931s | +| $th = 6$ | 0.00259s | 0.00696s | 0.01416s | 0.01626s | 0.03212s | +| $th = 8$ | 0.00191s | 0.00444s | 0.00784s | 0.01264s | 0.01446s | +| $th = 10$ | 0.00139s | 0.00346s | 0.00547s | 0.00776s | 0.01606s | +| $th = 12$ | 0.00123s | 0.00261s | 0.00444s | 0.00669s | 0.00908s | +| $th = 14$ | 0.00123s | 0.00253s | 0.00384s | 0.00541s | 0.00710s | +| $th = 16$ | 0.00098s | 0.00229s | 0.00341s | 0.00468s | 0.00553s | +| $th = 18$ | 0.00098s | 0.00218s | 0.00338s | 0.00467s | 0.00602s | +| $th = 20$ | 0.00099s | 0.00194s | 0.00284s | 0.00405s | 0.00563s | +| $th = 22$ | 0.00093s | 0.00173s | 0.00288s | 0.00388s | 0.00510s | +| $th = 24$ | 0.00104s | 0.00172s | 0.00258s | 0.00380s | 0.00516s | + +注:这里测出来的$th = 24$时似乎比上表中快了不少,这可能是因为每次测量时间的访存和分支规律都是一样的,所以运行次数变多时效率有一定提升。 + +加速比曲线图如下: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/3_4.png} +\end{center} + +并行效率很轻松地就超过了1。这很容易从算法复杂度角度分析出来。排序阶段的复杂度为$O(\frac{n^2}{p^2})$;假设并行的线程数不超过核心数,则归并阶段的复杂度为 + +$$\frac{1}{p}n + \frac{2}{p}n + ... + \frac{2^{\lceil \log_2 p \rceil}}{p}n = O(n)$$ + +因此总复杂度为$O(\frac{n^2}{p^2} + n)$。在$p$较小的时候可以达到接近平方级别的加速比。不过,这种加速是否真的是我们平常讨论的"加速"仍然存疑,本质上这样属于换了一种算法来实现。 + +如果$p$继续增大,超过核心数,可以预料一定范围内耗时仍会下降,因为$p$越大这个算法中计数排序的部分越少,归并排序的部分越多。 + +# Programming +我依然没有管那个main.c。本来gen.c中提供的代码是生成单精度浮点数的,但是根据助教的要求,还是使用了双精度浮点数。这挺可惜的,因为单精度浮点数的效率确实明显高于双精度浮点数。 + +运行了几次发现,主要的时间瓶颈在数据生成上,rand()函数的效率非常低。于是我使用了\href{https://en.wikipedia.org/wiki/Xorshift}{Xorshift}算法生成随机数。 + +## openmp版本 + +在代码中使用了runtime调度,这样比较方便测试不同的调度方式的速度差别。输入数据为$n = 48000 = 24 * 2000$,$th = 24$,选择这个$n$主要是为了分块的数据更好看一些。 + +### static调度 + +| 块大小 | 1 | 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | +|-------|---|----|----|----|----|----|----|----|----|----| +| 耗时 | 0.274s | 0.287s | 0.284s | 0.286s | 0.285s | 0.286s | 0.288s | 0.286s | 0.283s | 0.287s | + +| 块大小 | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 | +|-------|-----|-----|-----|-----|-----|-----|-----|-----|-----|------| +| 耗时 | 0.286s | 0.285s | 0.287s | 0.283s | 0.286s | 0.290s | 0.281s | 0.290s | 0.294s | 0.286s | + +| 块大小 | 1100 | 1200 | 1300 | 1400 | 1500 | 1600 | 1700 | 1800 | 1900 | 2000 | +|-------|------|------|------|------|------|------|------|------|------|------| +| 耗时 | 0.281s | 0.290s | 0.287s | 0.294s | 0.300s | 0.301s | 0.304s | 0.308s | 0.304s | 0.285s | + +可以观察到一定的周期现象。理论分析如下:假定每个块的工作量基本相同,考虑到块大小不是很小,虽然有随机数$b_i$的因素存在,但是工作量的方差也不大,这一条件是近似成立的。在此基础上,所有线程总共需进行$\lceil \frac{n}{blk*th} \rceil$轮计算,假设最后一轮计算中没有发生总剩余数目小于一个块的大小的现象,则总耗时正比于$\lceil \frac{n}{blk*th}\rceil*blk$。从中可以预测,当$blk*th$接近而小于一个能整除$n$的数时,耗时会较长,从表格中900,1900等处的数据可以基本验证,但是实际上它们之间的差距没有公式预言的这么大,我认为这可能是由以下原因导致的: + +1. 每个块的工作量的差距不可忽视 +2. 线程调度和同步的开销不可忽视 +3. 线程之间的计算并非完全独立,例如同一个物理核心超线程出来的两个逻辑核心共享一部分计算单元,尤其是在密集的浮点运算下,它们之间的干扰不可忽视 + +### dynamic调度 + +| 块大小 | 1 | 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | +|-------|---|----|----|----|----|----|----|----|----|----| +| 耗时 | 0.174s | 0.161s | 0.160s | 0.160s | 0.160s | 0.160s | 0.161s | 0.162s | 0.161s | 0.161s | + +| 块大小 | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 | +|-------|-----|-----|-----|-----|-----|-----|-----|-----|-----|------| +| 耗时 | 0.161s | 0.161s | 0.172s | 0.164s | 0.180s | 0.194s | 0.190s | 0.168s | 0.175s | 0.182s | + +| 块大小 | 1100 | 1200 | 1300 | 1400 | 1500 | 1600 | 1700 | 1800 | 1900 | 2000 | +|-------|------|------|------|------|------|------|------|------|------|------| +| 耗时 | 0.183s | 0.198s | 0.212s | 0.220s | 0.227s | 0.236s | 0.249s | 0.261s | 0.273s | 0.283s | + +可见dynamic调度性能普遍比static调度好一些,这证明线程间的工作分配确实是比较不均匀的。而对于$blk=2000$的情形,二者的速度几乎一样,证明dynamic调度的额外开销不是很大。综合考虑调度开销和线程间的工作分配,在$blk=100$附近达到的效果最好。 + +许多问题下,dynamic调度并不能和block+cyclic的static调度拉开差距,但是这个问题下差距比较明显,我认为这是因为在这个问题中block+cyclic的static调度与只有block的static调度没有本质的区别,每个线程的工作量仍然是随机生成的一系列整数,只是改变了一下具体负责哪一行而已。而很多其他问题中,工作量可能关于循环变量单调递增或递减,这样的block+cyclic的static调度就可以起到一定的负载均衡作用。 + +### guided调度 + +| 块大小 | 1 | 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | +|-------|---|----|----|----|----|----|----|----|----|----| +| 耗时 | 0.238s | 0.223s | 0.232s | 0.234s | 0.228s | 0.232s | 0.231s | 0.227s | 0.231s | 0.232s | + +| 块大小 | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 | +|-------|-----|-----|-----|-----|-----|-----|-----|-----|-----|------| +| 耗时 | 0.228s | 0.223s | 0.230s | 0.226s | 0.233s | 0.228s | 0.224s | 0.235s | 0.233s | 0.234s | + +| 块大小 | 1100 | 1200 | 1300 | 1400 | 1500 | 1600 | 1700 | 1800 | 1900 | 2000 | +|-------|------|------|------|------|------|------|------|------|------|------| +| 耗时 | 0.227s | 0.230s | 0.238s | 0.239s | 0.242s | 0.249s | 0.254s | 0.265s | 0.271s | 0.284s | + +可见guided调度速度普遍在dynamic调度和static调度之间,这也是符合对它的期望的。 + +三种调度方式的耗时曲线图如下: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/4_1.png} +\end{center} + +## pthread版本 + +(我实现的)pthread版本从任务分配的角度来看相当于openmp里的dynamic调度版本,也可以调节分块的大小,耗时与分块大小的关系如下: + +| 块大小 | 1 | 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | +|-------|---|----|----|----|----|----|----|----|----|----| +| 耗时 | 0.175s | 0.161s | 0.159s | 0.159s | 0.165s | 0.165s | 0.168s | 0.163s | 0.163s | 0.165s | + +| 块大小 | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 | +|-------|-----|-----|-----|-----|-----|-----|-----|-----|-----|------| +| 耗时 | 0.167s | 0.168s | 0.173s | 0.177s | 0.182s | 0.192s | 0.195s | 0.185s | 0.169s | 0.180s | + +| 块大小 | 1100 | 1200 | 1300 | 1400 | 1500 | 1600 | 1700 | 1800 | 1900 | 2000 | +|-------|------|------|------|------|------|------|------|------|------|------| +| 耗时 | 0.189s | 0.193s | 0.215s | 0.223s | 0.232s | 0.245s | 0.253s | 0.265s | 0.276s | 0.284s | + +因为它和openmp的dynamic调度非常相似,所以这里把它和openmp的dynamic调度的耗时结果汇成了一张图: + +\begin{center} +\includegraphics[width=1.0\linewidth]{pic/4_2.png} +\end{center} + +可见二者的耗时也非常相似。 diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/report.pdf" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/report.pdf" new file mode 100644 index 000000000..c89b72108 Binary files /dev/null and "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/report.pdf" differ diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/util.hpp" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/util.hpp" new file mode 100644 index 000000000..0ee7cd3da --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/hw/\346\235\216\346\231\250\346\230\2122017011466/hw5/util.hpp" @@ -0,0 +1,58 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using i8 = char; +using u8 = unsigned char; +using i32 = int; +using u32 = unsigned; +using u64 = unsigned long long; +using usize = std::size_t; +using f32 = float; +using f64 = double; + +struct Mutex { + pthread_mutex_t inner; + Mutex() { pthread_mutex_init(&inner, nullptr); } + Mutex(const Mutex &) = delete; + Mutex &operator=(const Mutex &) = delete; + void lock() { pthread_mutex_lock(&inner); } + void unlock() { pthread_mutex_unlock(&inner); } + ~Mutex() { pthread_mutex_destroy(&inner); } +}; + +struct CondVar { + pthread_cond_t inner; + CondVar() { pthread_cond_init(&inner, nullptr); } + void wait(Mutex *mu) { pthread_cond_wait(&inner, &mu->inner); } + void notify_one() { pthread_cond_signal(&inner); } + void notify_all() { pthread_cond_broadcast(&inner); } + ~CondVar() { pthread_cond_destroy(&inner); } +}; + +struct OmpMutex { + omp_lock_t inner; + OmpMutex() { omp_init_lock(&inner); } + OmpMutex(const Mutex &) = delete; + OmpMutex &operator=(const OmpMutex &) = delete; + void lock() { omp_set_lock(&inner); } + void unlock() { omp_unset_lock(&inner); } + ~OmpMutex() { omp_destroy_lock(&inner); } +}; + +template +struct Locker { + M *inner; + Locker(M *mu) : inner(mu) { inner->lock(); } + Locker(const Locker &) = delete; + Locker &operator=(const Locker &) = delete; + ~Locker() { inner->unlock(); } +}; \ No newline at end of file diff --git "a/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/readme.md" "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/readme.md" new file mode 100644 index 000000000..02e8ddb86 --- /dev/null +++ "b/\345\244\247\344\272\214\344\270\213/\351\253\230\346\200\247\350\203\275\350\256\241\347\256\227\345\257\274\350\256\272/readme.md" @@ -0,0 +1,15 @@ +总共五次作业,感觉我有点做的太认真了,应该不需要测这么多数据之类的就可以得到满分的。不过据说两个助教的给分严格程度差距较大,我没有可靠依据,仅供参考。 + +我的第一,五次作业扣了一些分: + +第一次: +``` +1.1 -0.5 n=10,p=8 +``` + +第五次: +``` +id = 2017011458 score = 95.8 1: 2.8/3.0 实现不够高效(如通讯子选择可以更好) 2: 2.8/3.0 实现不够高效(如可使用读写锁优化) 3: 3.0/3.0 4: 3.0/3.0 +``` + +考试很水,虽然没有往年题但是复习一下ppt应该完全没问题。