矩阵相乘 并行算法

并行处理技术课程设计分析报告课程设计题目矩阵相乘并行算法设计姓名廖杰学号M*********专业计算机技术任课教师金海石宣化所在学院计算机科学与技术学院报告提交日期2014-01-13一、实验目的1、学习使用集群；2、掌握并行处理或分布计算的编程方法；3、学会以并行处理的思想分析问题。

二、实验要求1、自行生成矩阵作为算法的输入；2、使用并行处理技术编程，例如：MPI、OpenMP、MR；3、矩阵大小至少为1000*1000；4、加速比越大成绩越高。

三、实验内容3.1、矩阵的划分：对于矩阵相乘的并行算法，可以有三种：对矩阵按行划分、按列划分和棋盘式分块划分。

和按行或列划分相比，棋盘式划分可以开发出更高的并行度。

对于一个n×n的方阵，棋盘划分最多可以使用n^2个处理器进行并行计算，但使用按行或列分解最多可以使用n个。

对矩阵相乘采用棋盘式划分的算法通常称作Cannon算法。

A）行列划分又叫带状划分（Striped Partitioning），就是将矩阵整行或者整列分成若干个组，每个组指派给一个处理器。

下图所例为4个CPU，8×8矩阵的带状划分。

在带状划分情况下，每个CPU将会均匀分配到2行(列)数据。

8×8矩阵变成了一个1×4或4×1的分块矩阵，每个CPU所属的分块矩阵大小为8×2或2×8。

B）棋盘划分就是将矩阵分成若干个子矩阵，每个子矩阵指派给一个处理器，此时任一处理器均不包含整行或者整列。

下图所示即为4个处理器情况下8×8矩阵的棋盘划分，其中处理器阵列为2×2，每个处理器分配到的子矩阵大小为4×4。

矩阵划分成棋盘状可以和处理器连成二维网孔相对应。

对于一个n×n维矩阵和p×p的二维处理器阵列，每个处理器均匀分配有（n/p）×(n/p)=n^2/p^2个元素。

使用棋盘式划分的矩阵相乘算法一般有两种，Cannon算法和Summa算法。

SUMMA算法能够计算m*l的A矩阵和l*n的B矩阵相乘（m、l、n可不相等），而cannon算法只能实现n*n的A矩阵和n*n 的B矩阵相乘，具有很大的局限性。

3.2、算法原理A) 行划分法假设是M*N，计算前，将矩阵N发送给所有从进程，然后将矩阵M分块，将M中数据按行分给各从进程，在从进程中计算M中部分行数据和N的乘积，最后将结果发送给主进程。

这里为了方便，有多少进程，就将M分了多少块，除最后一块外的其他数据块大小都相等，最后一块是剩下的数据，大小大于等于其他数据块大小，因为矩阵行数不一定整除进程数。

最后一块数据在主进程中计算，其他的在从进程中计算。

定义两个矩阵M和N，N所有进程都需要，M可以只在主进程中定义。

其他的变量视主进程和从进程需要按要求定义在合适的位置。

代码参见附录部分。

B) Cannon算法Cannon算法的基本思想可以如下表示：假设两个矩阵A和B相乘，把A和B矩阵划分成p 个方块，进程的编号从到，并在最初把子矩阵和分配给。

虽然第i行的每个进程需要全部的个子矩阵，但我们还是能调度第i行个进程的计算，使得每个进程在任何时刻都是用不同的。

每完成一次矩阵乘法，这些块在各进程之间被轮流使用，似的每次轮流之后每个进程都可以得到新的。

对列使用同样的调度，则在任何时刻，任何进程至多拥有每个矩阵的一个块，在所有进程中，改算法需要的总内存量为。

下图为此算法中不同进程上子矩阵乘法的调度过程。

假如矩阵C=A*B,则C的的计算公式如下：进程P 存储分块矩阵这一部分。

块矩阵乘法要计算所有匹配的和，然而只有在主对角线的才是匹配的。

因此需要采用循环移动分块矩阵的方法来使每个进程都有一对可以直接相乘的匹配的块，具体方法如下：(1)将排第i行的块循环左移i个位置，将第列．块循环上移j个位置；(2) 进程执行乘一加运算，然后将移动得到的块循环左移1个位置，将移动得到的块循环上移1个位置；(3)重复第2步(一1)次，每次移动后进程执行乘一加运算。

经过以上操作后就可以得到矩阵C的解。

代码请参见附录部分C) Summa算法SUMMA 算法首先将A , B 和C 划分为相同大小的矩阵，对应放在mesh_r × mesh_c 的二维mesh 上。

但SUMMA 算法将矩阵乘法分解为一系列的秩nb 修正, 即各处理器中的A 和B 分别被分解为nb 大小的列块和行块进行相乘，前面所说的分块尺寸就是指nb 的大小。

算法中, 广播实现为逻辑处理器行环或列环上的流水线传送, 达到了计算与通信的重叠. 具体描述如算法1所示。

C= 0for i= 0 t o k-1 step nb docur col = i×c/ ncur row = i×r / mif 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 forSUMMA算法的核心思想是：各处理器收集来自同一行处理器中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矩阵相乘，具有很大的局限性。

代码参见附录部分。

3.3、程序运行结果对比分析A) 统一的实验条件矩阵大小：1000*1000；矩阵数字范围：0~10；矩阵数字分布是否随机：是；分配的进程数：9；B) 实验进程数解释由于Cannon算法本身局限性，要使用Cannon算法，必须满足进程数为整数的平方，比如1、4、9、16等。

在本次的实验环境之下，经过多次对比分析，发现对于分行还是分块算法，进程数安排在8~15可以得到最理想的运行速度：进程数目过小则每个进程单独运算的时间过多，进程数过大则选路时间（进程与进程之间的通信时间）过长。

而对比要求每个算法的进程相同，故此处选择进程数目为9.C) 算法运行时间对比Cannon算法运行时间如下：分行法运行时间如下：串行算法运行时间如下：由于Summa算法与Cannon算法思路几乎相同，而且在算法预处理阶段要比Cannon算法更耗时，故没有做多余的实验。

CANNON显而易见，单纯的运用分行算法所花费的时间是最短的。

D) 结果分析Cannon算法相对于简单的行划分并行处理算法，其优势仅仅在于并行度可以更高（可达到N*N个进程，N为矩阵宽），但在并行度相同的情况下，其多出的预处理过程、矩阵发送与结果回收机制会占用更多的时间。

3.4、程序调优A) 行划分算法优化1. 循环优化在预估计矩阵大小为10的倍数的基础上，对每一个步长为1的循环做处理，改为步长为10的循环，将十次循环体全部压缩在一次循环中，从而大量减少了循环的判别时间，提升循环运算速度。

例如在单个线程在计算部分结果时，采用的循环为：for(i=0;i<line;i++){for(j=0;j<width;j++){DATA temp=0;for(k=0;k<width;k+=10){temp += buffer[i*width+k]*n[j*width+k];temp += buffer[i*width+k+1]*n[j*width+k+1];temp += buffer[i*width+k+2]*n[j*width+k+2];temp += buffer[i*width+k+3]*n[j*width+k+3];temp += buffer[i*width+k+4]*n[j*width+k+4];temp += buffer[i*width+k+5]*n[j*width+k+5];temp += buffer[i*width+k+6]*n[j*width+k+6];temp += buffer[i*width+k+7]*n[j*width+k+7];temp += buffer[i*width+k+8]*n[j*width+k+8];temp += buffer[i*width+k+9]*n[j*width+k+9];}ans[i*width+j] = temp;}}在将循环次数压缩的同时，为了进一步减少循环的运算量，在每一个步长为10的循环之前做预处理，避免循环体中的重复运算。

例如在主进程在接受其他进程时，将结果矩阵整合的过程：for(k=1;k<numprocs;k++){MPI_Recv(ans,line*width,MPI_INT,k,2,MPI_COMM_WORLD,&status);for(i=0;i<line;i++){count=i*k*width; //将i*k*width提前算好，减少了下一步循环的重复运算count1=i*width;for(j=0;j<width;j+=10){p[count+j] = ans[count1+j];p[count+j+1] = ans[count1+j+1];p[count+j+2] = ans[count1+j+2];p[count+j+3] = ans[count1+j+3];p[count+j+4] = ans[count1+j+4];p[count+j+5] = ans[count1+j+5];p[count+j+6] = ans[count1+j+6];p[count+j+7] = ans[count1+j+7];p[count+j+8] = ans[count1+j+8];p[count+j+9] = ans[count1+j+9];}}}2. 节省空间在进行矩阵工作量划分并传送的时候，为每一个进程开辟仅仅是自己所需要大小的空间，例如在9进程的环境下，每个进程所需要接受的缓存空间为B矩阵大小以及大约1/9大小A 矩阵。

内存开辟：buffer = (DATA *)malloc(sizeof(DATA)*width*line);矩阵A分块传输：for(k=1;k<numprocs;k++){for(i=k;i<width;i+=numprocs){count=i/numprocs*width;count1=i*width;for(j=0;j<width;j+=10){buffer[count+j]=m[count1+j];buffer[count+j+1]=m[count1+j+1];buffer[count+j+2]=m[count1+j+2];buffer[count+j+3]=m[count1+j+3];buffer[count+j+4]=m[count1+j+4];buffer[count+j+5]=m[count1+j+5];buffer[count+j+6]=m[count1+j+6];buffer[count+j+7]=m[count1+j+7];buffer[count+j+8]=m[count1+j+8];buffer[count+j+9]=m[count1+j+9];}}MPI_Send(buffer,line*width,MPI_INT,k,1,MPI_COMM_WORLD);同样的方式也运用在运行空间的开辟上。