稠密矩阵算法实现 LeeRinji

矩阵$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 维划分

原理分析

源代码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 维划分)

原理分析

源代码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 算法

原理分析

源代码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 算法

原理分析

源代码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 倍,运算效率处于同一个数量级)。