在科学与工程计算的许多问题中， 矩阵乘积是最基本的算法之一。在分布存储并行机上的经典矩阵乘积算法主要有1969年Cannon提出的二维mesh 上的矩阵乘积算法和1987年Fox等提出的“广播-乘积-滚动”算法。 1994年Choi 等提出的PUMMA 算法将Fox 算法推广到二维块卷帘数据分布上。同年，Huss-Lederman等又将Fox 算法推广到虚二维环绕数据分布。1995年van de Geijn 和Watts提出了一个简单而高效的矩阵乘积算法， 称为SUMMA 算法。
曙光4000并行机是国家智能计算机研究开发中心研制的一种基于分布存储结构和消息传送机制的大规模并行计算机系统。本文讨论了典型矩阵乘积算法SUMMA 算法在类似于曙光4000这样的大规模并行机上的实现。在后续阶段，作者将继续分析并比较它和其他矩阵乘法并行算法的性能。
SUMMA 算法首先将A , B 和C 划分为相同大小的矩阵，对应放在mesh_r × mesh_c 的二维mesh 上。 但SUMMA 算法将矩阵乘法分解为一系列的秩nb 修正, 即各处理器中的A 和B 分别被分解为nb 大小的列块和行块进行相乘，前面所说的分块尺寸就是指nb 的大小。算法中, 广播实现为逻辑处理器行环或列环上的流水线传送, 达到了计算与通信的重叠. 具体描述如算法1所示。
算法1. 二维mesh 上的SUMMA 算法
C= 0 for i= 0 t o k-1 step nb do cur col = i×c/ n cur row = i×r / m if my col = cur rol 向本行广播A 从i mod (k/c) 列开始的nb 列, 存于A′ if my row = cur row 向本列广播B 从i mod (k/r) 行开始的nb 行, 存于B ′ C= A′×B ′ end for
SUMMA算法的核心思想是：各处理器收集来自同一行处理器中A矩阵子块的所有列和同一列处理器中B矩阵子块的所有行，然后将行列相乘后累加，形成一个C矩阵的分块矩阵。最后由rank=0的处理器将其他处理器的数据收集起来，形成最终的矩阵C。
SUMMA算法相较于cannon算法的优势只要体现在SUMMA算法能够计算m*l的A矩阵和l*n的B矩阵相乘（m、l、n可不相等），而cannon算法只能实现n*n的A矩阵和n*n的B矩阵相乘，具有很大的局限性。不过由于水平有限，并且为了矩阵分块比较好分，此MPI实现的SUMMA算法要求输入的矩阵维数（m、l、n）必须为处理器个数（由用户自己定义）的整数倍，否则将提示错误。至于如何实现随意维数的矩阵相乘，还要请各位伟大的程序猿多多指教。
由于源码太长，这里只附上一部分MPI源码，需要全部源码的看客可邮件与我联系，联系方式在最下方。
/*summa算法 */ void My_summa(const int my_rank, int *a, int *b) { int k, i, j; int *c; int len_c; int dest, tag; int my_rank_row, my_rank_col; int *a_col_send, *b_row_send, *a_col_recv, *b_row_recv; MPI_Status status_a, status_b; my_rank_row = my_rank / mesh_c; my_rank_col = my_rank % mesh_c; a_col_send = (int *)malloc((m / mesh_r) * sizeof(int)); b_row_send = (int *)malloc((n / mesh_c) * sizeof(int)); a_col_recv = (int *)malloc((m / mesh_r) * sizeof(int)); b_row_recv = (int *)malloc((n / mesh_c) * sizeof(int)); /*各处理器将自己所拥有的A矩阵分块的列发送给跟自己同一行的其他处理器*/ for (k = 0; k < l / mesh_c; ++k) { for (i = 0; i < m / mesh_r; ++i) { a_col_send[i] = a[i * (l / mesh_c) + k]; } for (i = 0; i < mesh_c; ++i) { dest = my_rank_row * mesh_c + i; tag = my_rank_col * (l / mesh_c) + k; MPI_Send(a_col_send, m / mesh_r, MPI_INT, dest, tag, MPI_COMM_WORLD); } } /*各处理器将自己所拥有的B矩阵分块的行发送给跟自己同一列的其他处理器*/ for (k = 0; k < l / mesh_r; ++k) { for (i = 0; i < n / mesh_c; ++i) { b_row_send[i] = b[k * (n / mesh_c) + i]; } for (i = 0; i < mesh_r; ++i) { dest = my_rank_col + i * mesh_c; tag = my_rank_row * (l / mesh_r) + k; MPI_Send(b_row_send, n / mesh_c, MPI_INT, dest, tag, MPI_COMM_WORLD); } } /*初始化C分块 */ len_c = (m / mesh_r) * (n / mesh_c); c = (int *)malloc(len_c * sizeof(int)); for (i = 0; i < len_c; ++i) { c[i] = 0; } for (k = 0; k < l; ++k) { /*各个处理器接收来自同一行处理器中的所有A分块的列和同一列处理器中的所有B分块的行*/ MPI_Recv(a_col_recv, m / mesh_r, MPI_INT, my_rank_row * mesh_c + k / (l / mesh_c), k, MPI_COMM_WORLD, &status_a); MPI_Recv(b_row_recv, n / mesh_c, MPI_INT, k / (l / mesh_r) * mesh_c + my_rank_col, k, MPI_COMM_WORLD, &status_b); /*各个处理器已经拥有了计算对应C子块的所有信息，相乘累加后即可得到对应C子块*/ for (i = 0; i < m / mesh_r; ++i) { for (j = 0; j < n / mesh_c; ++j) { c[i * (n / mesh_c) + j] += a_col_recv[i] * b_row_recv[j]; } } } /*各处理器将C子块发送给0号进程*/ MPI_Send(c, len_c, MPI_INT, 0, 99, MPI_COMM_WORLD); free(c); free(a_col_send), free(b_row_send), free(a_col_recv), free(b_row_recv); } /*0号进程接收所有C子块信息*/ void Collect_C() { int k, i, j; int len_c, id_recv; int *c_recv; int c_imin, c_imax, c_jmin, c_jmax; MPI_Status status; len_c = (m / mesh_r) * (n / mesh_c); c_recv = (int *)malloc(len_c * sizeof(int)); /*将C子块填到相应的C矩阵坐标下 */ for (k = 0; k < num_procs; ++k) { MPI_Recv(c_recv, len_c, MPI_INT, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status); id_recv = status.MPI_SOURCE; c_imin = (id_recv / mesh_c) * (m / mesh_r); c_imax = c_imin + m / mesh_r - 1; c_jmin = (id_recv % mesh_c) * (n / mesh_c); c_jmax = c_jmin + n / mesh_c - 1; for (i = c_imin; i <= c_imax; ++i) { for (j = c_jmin; j <= c_jmax; ++j) { C[i][j] = c_recv[(i - c_imin) * (n / mesh_c) + (j - c_jmin)]; } } } free(c_recv); } /*打印矩阵 */ void Print(char *str, int **matrix, int m, int n) { int i, j; printf("%s", str); for (i = 0; i < m; ++i) { for (j = 0; j < n; ++j) { printf("%-6d", matrix[i][j]); } printf("\n"); } printf("\n"); }
编译： mpicc summa.c -o summa -lm 注：如出现p0_18495: p4_error: Child process exited while making connection to remote process on node2: 0
Killed by signal 2. 这种错误，解决办法：在主目录下直接输入：source /home/guest/.bashrc
运行： mpirun -np 9 summa 9 6 12
运行结果：
Allocate_A_B_C cheers! Random_A_B cheers! Scatter_A_B cheers! Collect_C cheers! random matrix A: 54 9 90 76 54 45 52 70 1 99 64 31 40 34 34 31 7 16 82 26 18 70 29 44 64 58 29 71 98 89 42 5 50 84 81 4 29 86 26 83 85 42 66 77 76 0 8 35 69 91 13 87 13 43 random matrix B: 83 77 1 12 0 51 53 94 56 3 78 90 60 8 28 38 43 65 81 9 42 9 61 50 97 30 41 10 17 54 53 52 84 6 17 84 58 70 79 66 26 57 8 38 17 36 76 60 1 9 21 43 19 83 94 16 13 87 78 83 94 32 35 78 90 4 62 48 75 93 67 53 Matrix C = A * B: 22444 14176 12709 12738 8969 17193 16835 15749 16331 12402 19294 24297 17333 13092 12303 14998 9607 18335 17209 11844 10776 12807 22936 21159 11967 7117 5542 5707 4419 8498 8574 7892 8342 3843 9746 11445 18337 13631 9227 11451 7755 13417 13420 14114 12063 9723 18818 19131 24187 14962 13659 19104 14705 21137 24925 16584 17612 20247 28026 28207 13965 11511 10709 10533 5148 16694 13815 11273 9543 10914 17401 20205 18936 11620 13315 16285 11693 20427 21139 11382 13086 15306 23702 23355 20768 9170 6731 7552 7905 13279 16685 12657 16043 5298 14106 18693 21549 14014 11801 14071 10513 16346 16301 13559 13651 9366 21661 20430 Print cheers! successful run time is 0.109375