矩阵$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}$
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 呈线性关系。
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 倍,运算效率处于同一个数量级)。
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 倍,运算效率处于同一个数量级)。
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 倍,运算效率处于同一个数量级)。