在科學與工程計算的許多問題中， 矩陣乘積是最基本的算法之一。在分布存儲並行機上的經典矩陣乘積算法主要有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