1. 实验内容与方法
1) 对矩阵进行分块,之后按块进行乘法,使用的是cannon算法。Cannon算法是一种存储有效的算法。为了使两矩阵下标满足相乘的要求,它和并行分块乘法不同,不是仅仅让B矩阵的各列块循环移动,而是有目的地让A的各行块以及B的各列块皆施行循环移位,从而实现对C的子块的计算。将矩阵A和B分成p个方块Aij和Bij,,每块大小为,并将它们分配给个处理器。开始时处理器Pij存放块Aij和Bij,并负责计算块Cij,然后算法开始执行:
⑴将块Aij向左循环移动i步;将块Bij向上循环移动j步;
⑵Pij执行乘加运算后将块Aij向左循环移动1步,块Bij向上循环移动1步;
⑶重复第⑵步,总共执行次乘加运算和次块Aij和Bij的循环单步移位。
2) 初始化数组。使用vector定义矩阵A、B和C,设置矩阵大小均为6000 * 6000。使用srand((unsigned) time(NULL))设置使用系统定时器的值作为随机种子,随机初始化矩阵A和矩阵B的值。
3) 程序计时。使用MPI自带的MPI_Wtime()函数,MPI_Wtime()记录的是墙上时间,在CPU时间的基础上包括了空闲等待时间。在各进程计算过程分块矩阵C之前和之后使用MPI_Wtime()函数计时,不计算最开始的初始化过程MPI发送的时间和最后接收数据汇总成矩阵C的时间。
2. 运行时间
矩阵A、B和C的大小均为6000 * 6000,且。经过代码运行计时,得到以下的实验结果表格。
处理器个数 | 用时/秒 |
串行 | 1734.048241 |
4核 | 321.548999 |
9核 | 194.883003 |
16核 | 111.360064 |
25核 | 61.248356 |
36核(3 * 12) | 35.795692 |
3. 加速比
处理器个数 | 加速比 |
串行 | 1 |
4核 | 5.39279627 |
9核 | 8.89789368 |
16核 | 15.5715449 |
25核 | 28.3117516 |
36核(3 * 12) | 48.4429311 |
4. 实验分析
从实验结果和加速比可以看出:
1 随着核数的增加,矩阵乘法所需要的时间在相应缩短;
2 矩阵乘法的加速比与运行的核心数近似成正比;
3 核心数越多,每个核分到的分块矩阵就越小,每个核计算速度就越快,加速比比核心数更高。
5. 源代码
#include "iostream"
#include "cstdlib"
#include "cstring"
#include "mpi.h"
#include "ctime"
#include "cmath"
#include "algorithm"
#include "vector"
using namespace std;
#define N 6000 //n是矩阵大小
// 全局变量声明
int block, block_size, p, sp;//, block是块的大小,block_size是块元素个数,等于block*block,p是处理器个数,sp是sqrt(p)
vector<vector<double>> A(N, vector<double>(N, 0)), B(N, vector<double>(N, 0)), C(N, vector<double>(N, 0));
double *a, *b, *c, *temp_a, *temp_b;
int my_rank, my_row, my_col;// my_rank是处理器的id,(my_row, my_col)是处理器逻辑阵列坐标
MPI_Status status;
double start_time, end_time;
int get_index(int row, int col, int sp);//处理器逻辑阵列坐标到rank的转换
void init_Matrix();//随机生成矩阵A和B
void scatter_A_B();//id为0的处理器向其他处理器发送矩阵A、B的分块
void init_alignment();//矩阵A和B的初始化,即A每行不同程度的左移,B每行不同程度的上移
void shift();//矩阵A每行左移一位,矩阵B每行上移一位,同时计算分块的中间结果C
void collect_C();//id为0的处理器从其余处理器收集分块矩阵C
void serial_Multiplication();//矩阵串行的计算
int main() {
MPI_Init(NULL, NULL);
MPI_Comm_size(MPI_COMM_WORLD, &p);
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
sp = sqrt(p);
if (sp * sp != p) {
if (my_rank == 0)
cout << "not a quadratic number" << endl;
MPI_Finalize();
exit(1);
}
block = N / sp;
block_size = block * block;
my_col = my_rank % sp;
my_row = my_rank / sp;
a = (double *) malloc(block_size * sizeof(double));
b = (double *) malloc(block_size * sizeof(double));
c = (double *) malloc(block_size * sizeof(double));;
for (int i = 0; i < block_size; ++i) {
c[i] = 0;
}
temp_a = (double *) malloc(block_size * sizeof(double));
temp_b = (double *) malloc(block_size * sizeof(double));
if (my_rank == 0) {
init_Matrix();
scatter_A_B();
} else {
MPI_Recv(a, block_size, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
MPI_Recv(b, block_size, MPI_FLOAT, 0, 2, MPI_COMM_WORLD, &status);
}
start_time = MPI_Wtime();
init_alignment();
shift();
end_time = MPI_Wtime();
if (my_rank == 0) {
cout << "matrix size : " << N << " * " << N << endl;
printf("time = %lf\n", end_time - start_time);
}
if (my_rank == 0) {
collect_C();
} else {
MPI_Send(c, block_size, MPI_FLOAT, 0, 1, MPI_COMM_WORLD);
}
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
return 0;
}
/*
处理器逻辑阵列坐标到rank的转换
输入:坐标(row, col),sqrt(p)
输出:rank
*/
int get_index(int row, int col, int sp) {
return ((row + sp) % sp) * sp + (col + sp) % sp;
}
/*
随机生成矩阵A和B
*/
void init_Matrix() {
srand((unsigned int) time(NULL)); //设置随机数种子
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
A[i][j] = rand();
B[i][j] = rand();
}
}
}
/*
id为0的处理器向其他处理器发送矩阵A、B的分块
*/
void scatter_A_B() {
int p_i_min, p_j_min, p_i_max, p_j_max;
for (int k = 0; k < p; k++) {
p_j_min = (k % sp) * block;
p_j_max = (k % sp + 1) * block - 1;
p_i_min = k / sp * block;
p_i_max = (k / sp + 1) * block - 1;
// id为0的处理器将矩阵A和B中的数据转换成一位数组传到temp中
for (int i = p_i_min; i <= p_i_max; ++i) {
for (int j = p_j_min; j <= p_j_max; j++) {
temp_a[(i - p_i_min) * block + j - p_j_min] = A[i][j];
temp_b[(i - p_i_min) * block + j - p_j_min] = B[i][j];
}
}
if (k == 0) { // id为0的处理器直接拷贝过去,其他处理器则发送过去
memcpy(a, temp_a, block_size * sizeof(double));
memcpy(b, temp_b, block_size * sizeof(double));
} else {
MPI_Send(temp_a, block_size, MPI_FLOAT, k, 1, MPI_COMM_WORLD);
MPI_Send(temp_b, block_size, MPI_FLOAT, k, 2, MPI_COMM_WORLD);
}
}
}
/*
将A中坐标为(i, j)的分块循环左移i步
将B中坐标为(i, j)的分块循环上移j步
*/
void init_alignment() {
MPI_Sendrecv(a, block_size, MPI_FLOAT, get_index(my_row, my_col - my_row, sp), 1, temp_a, block_size, MPI_FLOAT,
get_index(my_row, my_col + my_row, sp), 1, MPI_COMM_WORLD, &status);
memcpy(a, temp_a, block_size * sizeof(double));
MPI_Sendrecv(b, block_size, MPI_FLOAT, get_index(my_row - my_col, my_col, sp), 1, temp_b, block_size, MPI_FLOAT,
get_index(my_row + my_col, my_col, sp), 1, MPI_COMM_WORLD, &status);
memcpy(b, temp_b, block_size * sizeof(double));
}
/*
矩阵A每行左移一位,矩阵B每行上移一位,同时计算分块的中间结果C
*/
void shift() {
double total = 0;
for (int l = 0; l < sp; l++) {
for (int i = 0; i < block; i++) {
for (int j = 0; j < block; j++) {
total = c[i * block + j];
for (int k = 0; k < block; k++) {
total += a[i * block + k] * b[k * block + j];
}
c[i * block + j] = total;
}
}
// 将A全部循环左移1位
MPI_Sendrecv(a, block_size, MPI_FLOAT, get_index(my_row, my_col - 1, sp), 1, temp_a, block_size, MPI_FLOAT,
get_index(my_row, my_col + 1, sp), 1, MPI_COMM_WORLD, &status);
memcpy(a, temp_a, block_size * sizeof(double));
// 将B全部循环上移1位
MPI_Sendrecv(b, block_size, MPI_FLOAT, get_index(my_row - 1, my_col, sp), 1, temp_b, block_size, MPI_FLOAT,
get_index(my_row + 1, my_col, sp), 1, MPI_COMM_WORLD, &status);
memcpy(b, temp_b, block_size * sizeof(double));
}
}
/*
id为0的处理器从其余处理器收集分块矩阵C
*/
void collect_C() {
int p_i_min, p_j_min, p_i_max, p_j_max;
// 将id为0的处理器的分块矩阵c结果赋给C对应位置
for (int i = 0; i < block; i++) {
for (int j = 0; j < block; j++) {
C[i][j] = c[i * block + j];
}
}
for (int k = 1; k < p; k++) {
MPI_Recv(c, block_size, MPI_FLOAT, k, 1, MPI_COMM_WORLD, &status);
p_j_min = (k % sp) * block;
p_j_max = (k % sp + 1) * block - 1;
p_i_min = k / sp * block;
p_i_max = (k / sp + 1) * block - 1;
for (int i = p_i_min; i <= p_i_max; i++) {
for (int j = p_j_min; j <= p_j_max; ++j) {
C[i][j] = c[(i - p_i_min) * block + j - p_j_min];
}
}
}
}
void serial_Multiplication();//矩阵串行的计算 使用1个线程即串行的计算