稠密矩阵算法实现

post on 25 Oct 2019 about 17940words require 60min
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

矩阵$A$，$B$和向量$\vec{x}$按下面定义元素：

• $a_{ij}=\frac{i-0.1\times j+1}{i+j+1}$
• $b_{ij}=\frac{(j-0.2\times i+1)(i+j+1)}{i\times i+j\times j+1}$
• $x_i=\frac{i}{i\times i+1}$

完成矩阵向量乘$A\vec{x}$的 MPI 并行算法，要求对矩阵采用 2 维划分

原理分析

• 把$n\times n$矩阵$A$在$p$（$p\le n^2$）个进程中划分，$p$个进程组成阵列$\sqrt{p}\times\sqrt{p}$，把$n\times 1$向量$\vec{x}$分布在个进程阵列的最后一列
• $p=n^2$
  • 每个进程拥有矩阵$A$的一个元素。
  • 把$n\times 1$向量$\vec{x}$分布在$n$个进程的最后一列。
• $p<n^2$
  • 每个进程$P_i$储存矩阵$A$的具有$\frac{n}{\sqrt{p}}\times\frac{n}{\sqrt{p}}$阶的子矩阵$A_{i,j}$。
  • 最后一列进程$P_i(i=\sqrt{p}−1,2\sqrt{p}−1,\dots,𝑝−1)$储存向量$\vec{x}$的$\frac{n}{\sqrt{p}}$个元素$x_k$。
• 二维划分相比一维有更低的时间复杂度和更好的可扩展性。

源代码MatVecMul.c

#include <stdlib.h>
#include <math.h>
#include <mpi.h>
typedef double lf;
lf getA(int i, int j) { return (i - 0.1 * j + 1) / (i + j + 1); }
lf getB(int i, int j) { return (j - 0.2 * i + 1) * (i + j + 1) / (i * i + j * j + 1); }
lf getX(int i) { return i / (i * i + 1.0); }
int main(int argc, char **argv)
{
	MPI_Init(&argc, &argv);
	int comSize, comRank, n = 1 << atoi(argv[1]);
	MPI_Comm_size(MPI_COMM_WORLD, &comSize);
	MPI_Comm_rank(MPI_COMM_WORLD, &comRank);
	int sqSize = round(sqrt(comSize)),	 //进程数的开方
		sqLen = (n + sqSize - 1) / sqSize, //每个进程分到矩形的边长
		sqRow = comRank / sqSize,		   //每个进程分到矩形的行编号
		sqCol = comRank % sqSize;		   //每个进程分到矩形的列编号
	MPI_Comm rowComm, colComm;			   //按行和列划分的通信域
	MPI_Comm_split(MPI_COMM_WORLD, sqRow, sqCol, &rowComm);
	MPI_Comm_split(MPI_COMM_WORLD, sqCol, sqRow, &colComm);
	lf *x = (lf *)malloc(sizeof(lf) * sqLen),
	   *y = (lf *)malloc(sizeof(lf) * sqLen);
	if (!sqCol) //第一个通信步，0列的进程各自读取待求向量的一部分，发给对角线上的进程
	{
		for (int i = 0; i < sqLen; ++i)
			x[i] = getX(sqRow * sqLen + i);
		if (sqCol != sqRow) //防止仅有一个进程时发生自己发给自己造成阻塞
			MPI_Send(
				x,
				sqLen,
				MPI_DOUBLE,
				sqRow,
				1,
				rowComm);
	}
	else if (sqCol == sqRow) //对角线上（且不在第0列）的接受
		MPI_Recv(
			x,
			sqLen,
			MPI_DOUBLE,
			0,
			1,
			rowComm,
			MPI_STATUSES_IGNORE);
	MPI_Bcast( //第二个通信步，对角线上的进程将向量分给同列的其他进程
		x,
		sqLen,
		MPI_DOUBLE,
		sqCol,
		colComm);
	for (int j = 0; j < sqLen; ++j) //每个进程计算自己的矩阵和向量相乘的结果
		for (int i = y[j] = 0; i < sqLen; ++i)
			y[j] += getA(sqRow * sqLen + j, sqCol * sqLen + i) * x[i];
	MPI_Reduce( //第二个通信步，各进程将结果同步到第0列
		y,
		x,
		sqLen,
		MPI_DOUBLE,
		MPI_SUM,
		0,
		rowComm);
	MPI_Finalize();
}

调度脚本MatVecMul.pbs

这里只放出运行 256 进程时的调度脚本，其他同理。

#PBS -N MatVecMul
#PBS -l nodes=8:ppn=32
#PBS -j oe

source /public/software/profile.d/mpi_openmpi-intel-2.1.2.sh

mpicc $PBS_O_WORKDIR/MatVecMul.c -o $PBS_O_WORKDIR/MatVecMul -std=c11 -lm
for logN in 6 7 8 9 10 11 12 13 14 15 16 17
do
time mpiexec -machinefile $PBS_NODEFILE $PBS_O_WORKDIR/MatVecMul $logN
done

计算算法的运行时间，填表

由于要求矩阵二维分划，为方便并行的进程数取平方数。

详见附件中的MatVecMul.o6610~MatVecMul.o6614，分别对应了进程数 256、64、16、4、1 的运行时间。

可以看到，由于矩阵向量乘的时间复杂度是$O(n^2)$的，按照现代计算机单核1e9次每秒的运算速度来计算的话，除了进程数最小且问题规模最大的时候跑了15.105s，其他大部分都在1s上下，大部分都是并行的额外开销。而且可以看到这个额外开销随着进程数量增加而增加，在进程数$16\to 64$的时候有一个比较大的增幅，因为这时候从单机 16 核的配置换成了双机 64 核，通信开销增加了很多。

运行时间

问题规模\进程数	1	4	16	64	256
$2^6$	0.378	0.480	0.548	1.005	1.397
$2^7$	0.454	0.472	0.555	1.085	1.315
$2^8$	0.432	0.464	0.542	1.037	1.292
$2^9$	0.430	0.461	0.571	1.056	1.337
$2^{10}$	0.425	0.431	0.521	1.023	1.268
$2^{11}$	0.453	0.478	0.514	1.084	1.291
$2^{12}$	0.469	0.494	0.510	1.049	1.306
$2^{13}$	0.528	0.496	0.540	1.051	1.348
$2^{14}$	0.690	0.562	0.559	1.052	1.287
$2^{15}$	1.525	0.771	0.622	1.124	1.307
$2^{16}$	4.215	1.554	0.881	1.146	1.326
$2^{17}$	15.105	4.090	1.674	1.348	1.380

加速比

问题规模\进程数	1	4	16	64	256
$2^6$	1	0.7875	0.689781022	0.376119403	0.270579814
$2^7$	1	0.961864407	0.818018018	0.41843318	0.345247148
$2^8$	1	0.931034483	0.79704797	0.416586307	0.334365325
$2^9$	1	0.932754881	0.753064799	0.40719697	0.321615557
$2^{10}$	1	0.986078886	0.815738964	0.41544477	0.335173502
$2^{11}$	1	0.947698745	0.881322957	0.417896679	0.350890782
$2^{12}$	1	0.949392713	0.919607843	0.447092469	0.359111792
$2^{13}$	1	1.064516129	0.977777778	0.502378687	0.391691395
$2^{14}$	1	1.227758007	1.234347048	0.655893536	0.536130536
$2^{15}$	1	1.977950713	2.451768489	1.356761566	1.166794185
$2^{16}$	1	2.712355212	4.784335982	3.678010471	3.178733032
$2^{17}$	1	3.693154034	9.023297491	11.20548961	10.94565217

运算效率

问题规模\进程数	1	4	16	64	256
$2^6$	1	0.196875	0.043111314	0.005876866	0.001056952
$2^7$	1	0.240466102	0.051126126	0.006538018	0.001348622
$2^8$	1	0.232758621	0.049815498	0.006509161	0.001306115
$2^9$	1	0.23318872	0.04706655	0.006362453	0.001256311
$2^{10}$	1	0.246519722	0.050983685	0.006491325	0.001309271
$2^{11}$	1	0.236924686	0.055082685	0.006529636	0.001370667
$2^{12}$	1	0.237348178	0.05747549	0.00698582	0.00140278
$2^{13}$	1	0.266129032	0.061111111	0.007849667	0.001530045
$2^{14}$	1	0.306939502	0.077146691	0.010248337	0.00209426
$2^{15}$	1	0.494487678	0.153235531	0.021199399	0.00455779
$2^{16}$	1	0.678088803	0.299020999	0.057468914	0.012416926
$2^{17}$	1	0.923288509	0.563956093	0.175085775	0.042756454

验证等效率函数

见「运算效率」中标黑部分。

由于并行开销较大而计算量小，这只取问题规模较大的部分做一下验证。不考虑并行同步开销时$W(n)=\frac{n^2}{p^2}$，因此规模较大的时候要保持 E，则 n 与 p 呈线性关系。

简单矩阵乘法$C=AB$（矩阵采用 2 维划分）

原理分析

• 考虑两个$n\times n$矩阵$A$和$B$，分别划分成$p$个大小为$\frac{n}{\sqrt{p}}\times\frac{n}{\sqrt{p}}$的块$A_{i,j}$和$B_{i,j}$。
• 进程$P_{i,j}$最初储存$A_{i,j}$和$B_{i,j}$，并计算结果矩阵的块$C_{i,j}$。
• 计算子矩阵$C_{i,j}$需要所有子矩阵$A_{i,k}$和$B_{k,j}$（$0\le k\le\sqrt{p}$）。
• 每行进行矩阵$A$的块的全收集通信，同时每列进行矩阵$B$的块的全收集通信。
• 最后执行子矩阵的乘加运算。

源代码MatMatMul.c

#include <stdlib.h>
#include <math.h>
#include <mpi.h>
typedef double lf;
lf getA(int i, int j) { return (i - 0.1 * j + 1) / (i + j + 1); }
lf getB(int i, int j) { return (j - 0.2 * i + 1) * (i + j + 1) / (i * i + j * j + 1); }
lf getX(int i) { return i / (i * i + 1.0); }
int main(int argc, char **argv)
{
	MPI_Init(&argc, &argv);
	int comSize, comRank, n = 1 << atoi(argv[1]);
	MPI_Comm_size(MPI_COMM_WORLD, &comSize);
	MPI_Comm_rank(MPI_COMM_WORLD, &comRank);
	int sqSize = round(sqrt(comSize)),	 //进程数的开方
		sqLen = (n + sqSize - 1) / sqSize, //每个进程分到矩形的边长
		sqRow = comRank / sqSize,		   //每个进程分到矩形的行编号
		sqCol = comRank % sqSize;		   //每个进程分到矩形的列编号
	MPI_Comm rowComm, colComm;			   //按行和列划分的通信域
	MPI_Comm_split(MPI_COMM_WORLD, sqRow, sqCol, &rowComm);
	MPI_Comm_split(MPI_COMM_WORLD, sqCol, sqRow, &colComm);
	lf *c = (lf *)malloc(sizeof(lf) * sqLen * sqLen),
	   *ra = (lf *)malloc(sizeof(lf) * sqLen * sqLen),
	   *rb = (lf *)malloc(sizeof(lf) * sqLen * sqLen);
	for (int i = 0; i < sqLen; ++i)
		for (int j = 0; j < sqLen; ++j)
			c[i * sqLen + j] = getA(sqRow * sqLen + i, sqCol * sqLen + j);
	MPI_Alltoall( //按行分发a
		c,
		sqLen * sqLen / sqSize,
		MPI_DOUBLE,
		ra,
		sqLen * sqLen / sqSize,
		MPI_DOUBLE,
		rowComm);
	for (int i = 0; i < sqLen; ++i)
		for (int j = 0; j < sqLen; ++j)
			c[j * sqLen + i] = getB(sqRow * sqLen + i, sqCol * sqLen + j); //保存的是b的转置
	MPI_Alltoall(														   //按列分发b（转置）
		c,
		sqLen * sqLen / sqSize,
		MPI_DOUBLE,
		rb,
		sqLen * sqLen / sqSize,
		MPI_DOUBLE,
		colComm);
	for (int i = 0; i < sqLen; ++i)
		for (int j = 0; j < sqLen; ++j)
			for (int k = c[i * sqLen + j] = 0; k < sqLen; ++k)
				c[i * sqLen + j] += ra[i * sqLen + k] * rb[j * sqLen + k]; //b中是转置
	free(c), free(ra), free(rb);
	MPI_Finalize();
}

调度脚本MatMatMul.pbs

#PBS -N MatMatMul
#PBS -l nodes=8:ppn=32
#PBS -j oe

source /public/software/profile.d/mpi_openmpi-intel-2.1.2.sh

mpicc $PBS_O_WORKDIR/MatMatMul.c -o $PBS_O_WORKDIR/MatMatMul -std=c11 -lm
for logN in 6 7 8 9 10 11 12 13 14 15 16 17
do
time mpiexec -machinefile $PBS_NODEFILE $PBS_O_WORKDIR/MatMatMul $logN
done

计算算法的运行时间，填表

由于要求矩阵二维分划，因此并行的进程数必须是平方数。

详见附件中的MatMatMul.o6633~MatMatMul.o6638，分别对应了进程数 256、64、16、4、1 的运行时间。下表中的x代表没有在限定时间中求解完毕。

可以看到，由于矩阵向量乘的时间复杂度是$O(n^3)$的，按照单核1e9次每秒的运算速度来计算的话，问题规模稍大一点的时候，在一个作业脚本的一个小时限制内基本上是跑不完的。当然，可以看到这里同问题规模且足够大的时候，运行时间会随进程数增加而减少，并行化取得了一定的加速结果。

运行时间

问题规模\进程数	1	4	16	64	256
$2^6$	0.543	0.448	0.587	1.245	1.407
$2^7$	0.496	0.433	0.572	1.165	1.407
$2^8$	0.530	0.477	0.575	1.136	1.375
$2^9$	0.520	0.507	0.562	1.181	1.597
$2^{10}$	1.208	0.539	0.569	1.096	1.664
$2^{11}$	7.368	1.140	0.685	1.151	1.665
$2^{12}$	61.114	7.279	1.473	1.223	1.657
$2^{13}$	481.676	61.072	7.449	2.231	1.712
$2^{14}$	x	485.722	60.931	11.028	3.106
$2^{15}$	x	x	487.180	80.777	15.511
$2^{16}$	x	x	x	636.016	111.355
$2^{17}$	x	x	x	x	875.914

加速比

问题规模\进程数	1	4	16	64	256
$2^6$	1	1.212053571	0.925042589	0.436144578	0.385927505
$2^7$	1	1.145496536	0.867132867	0.425751073	0.352523099
$2^8$	1	1.111111111	0.92173913	0.466549296	0.385454545
$2^9$	1	1.025641026	0.925266904	0.440304826	0.32561052
$2^{10}$	1	2.241187384	2.123022847	1.102189781	0.725961538
$2^{11}$	1	6.463157895	10.75620438	6.401390096	4.425225225
$2^{12}$	1	8.395933507	41.48947726	49.97056419	36.88231744
$2^{13}$	1	7.887018601	64.66317627	215.9013895	281.3528037
$2^{14}$	x	x	x	x	x
$2^{15}$	x	x	x	x	x
$2^{16}$	x	x	x	x	x
$2^{17}$	x	x	x	x	x

运算效率

问题规模\进程数	1	4	16	64	256
$2^6$	1	0.303013393	0.057815162	0.006814759	0.001507529
$2^7$	1	0.286374134	0.054195804	0.006652361	0.001377043
$2^8$	1	0.277777778	0.057608696	0.007289833	0.001505682
$2^9$	1	0.256410256	0.057829181	0.006879763	0.001271916
$2^{10}$	1	0.560296846	0.132688928	0.017221715	0.002835787
$2^{11}$	1	1.615789474	0.672262774	0.10002172	0.017286036
$2^{12}$	1	2.098983377	2.593092329	0.780790065	0.144071553
$2^{13}$	1	1.97175465	4.041448517	3.373459211	1.09903439
$2^{14}$	x	x	x	x	x
$2^{15}$	x	x	x	x	x
$2^{16}$	x	x	x	x	x
$2^{17}$	x	x	x	x	x

验证等效率函数

不考虑并行开销时（根据古斯塔夫森定律，问题规模大的时候可忽略），等效率函数$W(n)=\frac{n^3}{p^2}$，即$n\to p^{\frac{2}{3}}$。见「运算效率」中标黑部分（进程数增加了 64 倍，问题规模增加 16 倍，运算效率处于同一个数量级）。

矩阵乘法$C=AB$的 CANNON 算法

原理分析

• 最初子矩阵$A_{i,j}$和$B_{i,j}$分配给进程$P_{i,j}$。
• 我们调度第$i$行$\sqrt{p}$个进程的计算，使得每个进程在不同时刻使用不同的$A_{i,k}$。
• 每完成一次子矩阵乘法，这些块在各进程之间进行一次移位，使得每个进程获得下一步需要的新的$A_{i,k}$。
• 对列使用同样的调度。
• 该算法需要的内存总量为$O(n^2)$。
• Cannon 算法的成本最优条件和等效率函数与前面提到的简单方法一致，但更节省空间。

源代码Cannon.c

#include <stdlib.h>
#include <math.h>
#include <mpi.h>
typedef double lf;
lf getA(int i, int j) { return (i - 0.1 * j + 1) / (i + j + 1); }
lf getB(int i, int j) { return (j - 0.2 * i + 1) * (i + j + 1) / (i * i + j * j + 1); }
lf getX(int i) { return i / (i * i + 1.0); }
int main(int argc, char **argv)
{
	MPI_Init(&argc, &argv);
	int comSize, comRank;
	MPI_Comm_size(MPI_COMM_WORLD, &comSize);
	MPI_Comm_rank(MPI_COMM_WORLD, &comRank);

	int dims[2] = {round(sqrt(comSize)), round(sqrt(comSize))}, //进程矩阵长宽
		periods[2] = {1, 1},
		n = 1 << atoi(argv[1]), //矩阵阶数
		sqLen = n / dims[0],	//每个进程分配矩阵长度
		sqRank,
		sqCoord[2];

	MPI_Comm cartComm;
	MPI_Cart_create( //创建笛卡尔拓扑
		MPI_COMM_WORLD,
		2,
		dims,
		periods, //每个维度上循环
		0,
		&cartComm);
	MPI_Comm_rank(cartComm, &sqRank); //获得笛卡尔拓扑下的编号
	MPI_Cart_coords(				  //获得笛卡尔坐标
		cartComm,
		sqRank,
		2,
		sqCoord);

	lf *a = (lf *)malloc(sqLen * sqLen * sizeof(lf)),
	   *b = (lf *)malloc(sqLen * sqLen * sizeof(lf)),
	   *c = (lf *)malloc(sqLen * sqLen * sizeof(lf));

	for (int i = 0; i < sqLen; ++i)
		for (int j = 0; j < sqLen; ++j)
		{ // 初始化矩阵内容
			a[i * sqLen + j] = getA(sqCoord[0] * sqLen + i, sqCoord[1] * sqLen + j);
			b[i * sqLen + j] = getB(sqCoord[0] * sqLen + i, sqCoord[1] * sqLen + j);
			c[i * sqLen + j] = 0;
		}

	int shiftsource, shiftdest;
	MPI_Cart_shift( // 将a、b分发给相应进程
		cartComm,
		0,
		-sqCoord[1],
		&shiftsource,
		&shiftdest);
	MPI_Sendrecv_replace(
		a,
		sqLen * sqLen,
		MPI_DOUBLE,
		shiftdest,
		1,
		shiftsource,
		1,
		cartComm,
		MPI_STATUSES_IGNORE);
	MPI_Cart_shift(
		cartComm,
		1,
		-sqCoord[0],
		&shiftsource,
		&shiftdest);
	MPI_Sendrecv_replace(
		b,
		sqLen * sqLen,
		MPI_DOUBLE,
		shiftdest,
		1,
		shiftsource,
		1,
		cartComm,
		MPI_STATUSES_IGNORE);

	int uRank, dRank, lRank, rRank; //上下左右，其中右、下为源
	MPI_Cart_shift(
		cartComm,
		0,
		-1,
		&rRank,
		&lRank);
	MPI_Cart_shift(
		cartComm,
		1,
		-1,
		&dRank,
		&uRank);

	for (int i = 0; i < dims[0]; ++i)
	{
		for (int i = 0; i < sqLen; ++i) //本地乘加
			for (int j = 0; j < sqLen; ++j)
				for (int k = 0; k < sqLen; ++k)
					c[i * sqLen + j] += a[i * sqLen + k] * b[k * sqLen + j];
		MPI_Sendrecv_replace(
			a,
			sqLen * sqLen,
			MPI_DOUBLE,
			lRank,
			1,
			rRank,
			1,
			cartComm,
			MPI_STATUSES_IGNORE);
		MPI_Sendrecv_replace(
			b,
			sqLen * sqLen,
			MPI_DOUBLE,
			uRank,
			1,
			dRank,
			1,
			cartComm,
			MPI_STATUSES_IGNORE);
	}

	MPI_Cart_shift( // 从相应进程处恢复发送的a、b子块
		cartComm,
		0,
		-sqCoord[1],
		&shiftsource,
		&shiftdest);
	MPI_Sendrecv_replace(
		a,
		sqLen * sqLen,
		MPI_DOUBLE,
		shiftdest,
		1,
		shiftsource,
		1,
		cartComm,
		MPI_STATUSES_IGNORE);
	MPI_Cart_shift(
		cartComm,
		1,
		-sqCoord[0],
		&shiftsource,
		&shiftdest);
	MPI_Sendrecv_replace(
		b,
		sqLen * sqLen,
		MPI_DOUBLE,
		shiftdest,
		1,
		shiftsource,
		1,
		cartComm,
		MPI_STATUSES_IGNORE);
	free(a), free(b), free(c);
	MPI_Comm_free(&cartComm);
	MPI_Finalize();
}

调度脚本Cannon.pbs

这里只放出运行 256 进程时的调度脚本，其他同理。

#PBS -N Cannon
#PBS -l nodes=8:ppn=32
#PBS -j oe

source /public/software/profile.d/mpi_openmpi-intel-2.1.2.sh

mpicc $PBS_O_WORKDIR/Cannon.c -o $PBS_O_WORKDIR/Cannon -std=c11 -lm
for logN in 6 7 8 9 10 11 12 13 14 15 16 17
do
time mpiexec -machinefile $PBS_NODEFILE $PBS_O_WORKDIR/Cannon $logN
done

计算算法的运行时间，填表

由于要求矩阵二维分划，因此并行的进程数必须是平方数。

详见附件中的Cannon.o6645~Cannon.o6649，分别对应了进程数 256、64、16、4、1 的运行时间。下表中的x代表没有在限定时间中求解完毕。

由于 Cannon 算法相对于简单并行矩阵乘法相比减少的只是空间占用，时间复杂度没有改变，并且由于算法变得更加复杂，时间常数变的比较大，当问题规模达到$2^{17}$的时候均不能在规定时间内计算出结果。其他情况下的运行结果和简单矩阵乘法类似，不再赘述。

运行时间

问题规模\进程数	1	4	16	64	256
$2^6$	0.499	0.476	0.791	1.600	1.794
$2^7$	0.457	0.452	0.583	1.167	1.415
$2^8$	0.491	0.483	0.563	1.171	1.347
$2^9$	0.580	0.562	0.496	1.142	1.347
$2^{10}$	1.143	0.744	0.640	1.027	1.339
$2^{11}$	8.048	1.846	0.883	1.144	1.394
$2^{12}$	61.985	15.506	4.000	1.801	1.495
$2^{13}$	490.234	124.022	30.130	10.659	2.781
$2^{14}$	x	980.345	247.967	81.121	25.447
$2^{15}$	x	x	1972.293	648.052	217.948
$2^{16}$	x	x	x	x	1751.185
$2^{17}$	x	x	x	x	x

加速比

问题规模\进程数	1	4	16	64	256
$2^6$	1	1.048319328	0.630847029	0.311875	0.278149387
$2^7$	1	1.011061947	0.783876501	0.391602399	0.322968198
$2^8$	1	1.016563147	0.872113677	0.419299744	0.364513734
$2^9$	1	1.03202847	1.169354839	0.507880911	0.430586488
$2^{10}$	1	1.536290323	1.7859375	1.112950341	0.853622106
$2^{11}$	1	4.359696641	9.114382786	7.034965035	5.773314204
$2^{12}$	1	3.997484845	15.49625	34.41699056	41.46153846
$2^{13}$	1	3.952798697	16.27062728	45.99249461	176.2797555
$2^{14}$	x	x	x	x	x
$2^{15}$	x	x	x	x	x
$2^{16}$	x	x	x	x	x
$2^{17}$	x	x	x	x	x

运行效率

问题规模\进程数	1	4	16	64	256
$2^6$	1	0.262079832	0.039427939	0.004873047	0.001086521
$2^7$	1	0.252765487	0.048992281	0.006118787	0.001261595
$2^8$	1	0.254140787	0.054507105	0.006551558	0.001423882
$2^9$	1	0.258007117	0.073084677	0.007935639	0.001681978
$2^{10}$	1	0.384072581	0.111621094	0.017389849	0.003334461
$2^{11}$	1	1.08992416	0.569648924	0.109921329	0.022552009
$2^{12}$	1	0.999371211	0.968515625	0.537765478	0.161959135
$2^{13}$	1	0.988199674	1.016914205	0.718632728	0.688592795
$2^{14}$	x	x	x	x	x
$2^{15}$	x	x	x	x	x
$2^{16}$	x	x	x	x	x
$2^{17}$	x	x	x	x	x

验证等效率模型

不考虑并行开销时（根据古斯塔夫森定律，问题规模大的时候可忽略），等效率函数$W(n)=\frac{n^3}{p^2}$，即$n\to p^{\frac{2}{3}}$。见「运算效率」中标黑部分（进程数增加了 64 倍，问题规模增加 16 倍，运算效率处于同一个数量级）。

矩阵乘法$C=AB$的 DNS 算法

原理分析

• 使用三维划分。
• 我们知道结果矩阵 C 中的每个元素实际上是由$A$的行和$B$的列进行一个点乘得到的。把点乘操作的多次乘法运算展开成第三维。
• 因此，在进程立方中每个进程只负责了一次乘加操作。而二维算法中负责一次点乘操作。
• DNS 算法最多允许$n^3$个进程参与计算。
• 假设有一个$n\times n\times n$的进程阵列。
• 将$A$的第$j$列发送到第$j$层的某一列进程中，将$B$的第$i$行发送到第$i$层的某一行进程中，然后对同层进行广播。
• 每个处理器计算一个单独的乘加操作。
• 沿着第三维做加法归约，算出具体对应的矩阵$C$的每一个元素。
• 一个乘加操作运行时间是一个常数，而归约和广播复杂度为$O(\log n)$，因此总体需要的并行运行时间为$O(\log n)$。
• 使用$n^3$个进程是并不是成本最优的。
• 假设进程的数量是完全立方数$p=q^3$。
• 两个输入矩阵被划分成$\frac{n}{q}\times\frac{n}{q}$的块。
• 因此矩阵可以被看做是由这些块排成的$q\times q$二维方阵。
• 算法和前面所述是完全类似的，唯一不同的是每个进程现在对子块操作而不是对矩阵的单个元素操作。

源代码DNS.c

#include <stdlib.h>
#include <math.h>
#include <mpi.h>
typedef double lf;
lf getA(int i, int j) { return (i - 0.1 * j + 1) / (i + j + 1); }
lf getB(int i, int j) { return (j - 0.2 * i + 1) * (i + j + 1) / (i * i + j * j + 1); }
lf getX(int i) { return i / (i * i + 1.0); }
int main(int argc, char **argv)
{
	MPI_Init(&argc, &argv);
	int comSize, comRank, n = 1 << atoi(argv[1]);
	MPI_Comm_size(MPI_COMM_WORLD, &comSize);
	MPI_Comm_rank(MPI_COMM_WORLD, &comRank);

	int sqSize = round(pow(comSize, 1.0 / 3)), //进程数的开立方
		sqLen = (n + sqSize - 1) / sqSize,	 //每个进程分到矩形的边长
		hRank = comRank / sqSize / sqSize,
		rRank = comRank / sqSize % sqSize,
		cRank = comRank % sqSize,
		sqHet = rRank * sqSize + cRank,
		sqRow = cRank * sqSize + hRank,
		sqCol = hRank * sqSize + rRank;

	MPI_Comm hetComm, rowComm, colComm; //按高分的通信域是一条直线，按行、列分的通信域是一个二维平面
	MPI_Comm_split(MPI_COMM_WORLD, sqHet, hRank, &hetComm);
	MPI_Comm_split(MPI_COMM_WORLD, sqRow, rRank, &rowComm);
	MPI_Comm_split(MPI_COMM_WORLD, sqCol, cRank, &colComm);

	lf *c = (lf *)malloc(sizeof(lf) * sqLen * sqLen),
	   *a = (lf *)malloc(sizeof(lf) * sqLen * sqLen),
	   *b = (lf *)malloc(sizeof(lf) * sqLen * sqLen);
	if (!hRank)
		for (int i = 0; i < sqLen; ++i)
			for (int j = 0; j < sqLen; ++j)
			{
				a[i * sqLen + j] = getA(rRank * sqLen + i, cRank * sqLen + j);
				b[j * sqLen + i] = getB(rRank * sqLen + i, cRank * sqLen + j); //b中存的是转置
			}
	MPI_Bcast( //按行广播A
		a,
		sqLen * sqLen,
		MPI_DOUBLE,
		0,
		hetComm);
	MPI_Bcast( //按列广播B
		b,
		sqLen * sqLen,
		MPI_DOUBLE,
		0,
		hetComm);
	/*//相当于每个进程执行下面的代码
	for (int i = 0; i < sqLen; ++i)
		for (int j = 0; j < sqLen; ++j)
		{
			a[i * sqLen + j] = getA(sqRow * sqLen + i, sqHet * sqLen + j);
			b[i * sqLen + j] = getB(sqHet * sqLen + i, sqCol * sqLen + j);
		}
	*/
	for (int i = 0; i < sqLen; ++i)
		for (int j = 0; j < sqLen; ++j)
			for (int k = c[i * sqLen + j] = 0; k < sqLen; ++k)
				c[i * sqLen + j] += a[i * sqLen + k] * b[j * sqLen + k];

	MPI_Reduce( //规约答案到最下面一层
		hRank ? c : MPI_IN_PLACE,
		c,
		sqLen * sqLen,
		MPI_DOUBLE,
		MPI_SUM,
		0,
		hetComm);
	free(a), free(b), free(c);
	MPI_Finalize();
}

调度脚本DNS.pbs

这里只放出运行 512 进程时的调度脚本，其他同理。

#PBS -N DNS
#PBS -l nodes=16:ppn=32
#PBS -j oe

source /public/software/profile.d/mpi_openmpi-intel-2.1.2.sh

mpicc $PBS_O_WORKDIR/DNS.c -o $PBS_O_WORKDIR/DNS -std=c11 -lm
for logN in 6 7 8 9 10 11 12 13 14 15 16 17
do
time mpiexec -machinefile $PBS_NODEFILE $PBS_O_WORKDIR/DNS $logN
done

计算算法的运行时间，填表

由于算法的前提要求，并行的进程数必须是立方数。

详见附件中的DNS.o6678~DNS.o6681，分别对应了进程数 512、64、8、1 的运行时间。下表中的x代表没有在限定时间中求解完毕。

DNS 算法相对 Cannon 算法来说，虽然时间复杂度没有改变，但是常数更大，运行速度更慢（可以对比 64 进程的运行时间）。其他情况下的运行结果和简单矩阵乘法类似，不再赘述。

运行时间

问题规模\进程数	1	8	64	512
$2^6$	0.488	0.365	1.497	1.673
$2^7$	0.489	0.569	1.123	1.445
$2^8$	0.401	0.516	1.116	1.416
$2^9$	0.436	0.450	1.100	1.445
$2^{10}$	0.853	0.594	1.039	1.416
$2^{11}$	6.585	1.259	1.149	1.448
$2^{12}$	52.309	6.933	2.580	1.567
$2^{13}$	413.207	54.147	14.935	2.699
$2^{14}$	x	424.275	112.297	12.109
$2^{15}$	x	x	882.313	66.739
$2^{16}$	x	x	x	678.642
$2^{17}$	x	x	x	x

加速比

问题规模\进程数	1	8	64	512
$2^6$	1	1.336986301	0.325985304	0.291691572
$2^7$	1	0.85940246	0.435440784	0.338408304
$2^8$	1	0.777131783	0.359318996	0.28319209
$2^9$	1	0.968888889	0.396363636	0.301730104
$2^{10}$	1	1.436026936	0.820981713	0.60240113
$2^{11}$	1	5.230341541	5.731070496	4.547651934
$2^{12}$	1	7.544930045	20.2748062	33.38162093
$2^{13}$	1	7.631207638	27.66702377	153.096332
$2^{14}$	x	x	x	x
$2^{15}$	x	x	x	x
$2^{16}$	x	x	x	x
$2^{17}$	x	x	x	x

运行效率

问题规模\进程数	1	8	64	512
$2^6$	1	0.167123288	0.00509352	0.00056971
$2^7$	1	0.107425308	0.006803762	0.000660954
$2^8$	1	0.097141473	0.005614359	0.00055311
$2^9$	1	0.121111111	0.006193182	0.000589317
$2^{10}$	1	0.179503367	0.012827839	0.001176565
$2^{11}$	1	0.653792693	0.089547977	0.008882133
$2^{12}$	1	0.943116256	0.316793847	0.065198478
$2^{13}$	1	0.953900955	0.432297246	0.299016273
$2^{14}$	x	x	x	x
$2^{15}$	x	x	x	x
$2^{16}$	x	x	x	x
$2^{17}$	x	x	x	x

验证等效率模型

不考虑并行开销时（根据古斯塔夫森定律，问题规模大的时候可忽略），等效率函数$W(n)=\frac{n^3}{p^2}$，即$n\to p^{\frac{2}{3}}$。见「运算效率」中标黑部分（进程数增加了 8 倍，问题规模增加 4 倍，运算效率处于同一个数量级）。

Related posts

• Summary Of 2024 Fall 19 Dec 2024
• PrefixSpan数据挖掘算法 25 Jun 2024
• 现代偏微分方程计算方法课程项目 21 Jun 2024

Loading comments...