矩阵相乘-并行算法

并行处理技术并行处理技术 课程设计分析报告课程设计分析报告 课程设计题目矩阵相乘并行算法设计 姓名廖杰 学号M201372880 专业计算机技术 任课教师金海 石宣化 所在学院计算机科学与技术学院 报告提交日期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 行 0 1 的每个进程需要全部的个子矩阵 但我们还是能调度第 i 行个进程的计算 使得 每个进程在任何时刻都是用不同的 每完成一次矩阵乘法 这些块在各进程之间被轮 流使用 似的每次轮流之后每个进程都可以得到新的 对列使用同样的调度 则在任 何时刻 任何进程至多拥有每个矩阵的一个块 在所有进程中 改算法需要的总内存量为 下图为此算法中不同进程上子矩阵乘法的调度过程 2 假如矩阵 C A B 则 C 的 的计算公式如下 1 0 0 1 0 1 进程 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 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 矩阵相乘 具有很大的局限性 代码参见附录部分 3 3 程序运行结果对比分析 A 统一的实验条件 矩阵大小 1000 1000 矩阵数字范围 0 10 矩阵数字分布是否随机 是 分配的进程数 9 B 实验进程数解释 由于 Cannon 算法本身局限性 要使用 Cannon 算法 必须满足进程数为整数的平方 比如 1 4 9 16 等 在本次的实验环境之下 经过多次对比分析 发现对于分行还是分块算 法 进程数安排在 8 15 可以得到最理想的运行速度 进程数目过小则每个进程单独运算的 时间过多 进程数过大则选路时间 进程与进程之间的通信时间 过长 而对比要求每个 算法的进程相同 故此处选择进程数目为 9 C 算法运行时间对比 Cannon 算法运行时间如下 分行法运行时间如下 串行算法运行时间如下 由于 Summa 算法与 Cannon 算法思路几乎相同 而且在算法预处理阶段要比 Cannon 算法 更耗时 故没有做多余的实验 算法算法分行分行CANNON串行串行 时间时间1 218810s0 76s10 420s 显而易见 单纯的运用分行算法所花费的时间是最短的 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 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 同样的方式也运用在运行空间的开辟上 这样做不仅仅是内存空间的节约 同时也减少了进程之间的数据传输量 大大节省了进程 之间的协作时间 B 稀疏矩阵处理 虽然程序并未对稀疏矩阵进行优化 但是还是试着对程序的输入数据模式进行更改 体验 一下稀疏矩阵运算的速度提升有多快 void magten DATA a int width int i j srand unsigned time NULL for i 0 i width i for j 0 j width j 10 a i width j rand MUL rand 2 rand 3 a i width j 1 rand MUL rand 2 rand 3 a i width j 2 rand MUL rand 2 rand 3 a i width j 3 rand MUL rand 2 rand 3 a i width j 4 rand MUL rand 2 rand 3 a i width j 5 rand MUL rand 2 rand 3 a i width j 6 rand MUL rand 2 rand 3 a i width j 7 rand MUL rand 2 rand 3 a i width j 8 rand MUL rand 2 rand 3 a i width j 9 rand MUL rand 2 rand 3 如上面所示 对于每一个生成的数据都再一次进行乘法运算 其中乘数是 0 的概率为 2 3 运行结果如下 可以看出 稀疏矩阵的乘法由 0 76s 变为 0 75s 仅仅是短暂的提升 提升计时显示的精度 可以看到 对于稀疏矩阵的处理要比普通矩阵快 0 015s 提速约为 2 3 5 算法最优速度 对算法运用不同的进程数目运算进行了大量重复试验 最终得出在进程数大概为 12 的时候 本算法的运行速度最快 最终结果如下 发出工作量时间为 0 138184s 运算时间为 0 495569s 接收答案时间为 0 025461s 总运算时间 0 659240s 四 总结分析 串行串行 并行并行行划分行划分串行串行 运行时间运行时间0 659240s10 420s 进程数目进程数目121 加速比加速比15 801 效率效率15 80 12 1 3171 从效率大于 1 上可以看出 本次课程设计做出的算法为超线性加速 这主要得益于对循环 体的优化 附录 Cannon 算法 include include mpi h include include include include define MUL 10 MPI Status status double A B C C A B double a b c 各个进程的缓冲区 int n 矩阵的行列数 int np 每个进程控制的小矩阵的行列数 int p rank 进程个个数 当前进程的编号 笛卡尔进程编号 double tempa tempb void ProduceABC 在根处理器中生成矩阵 AB 初始化矩阵 C void PrintABC 输出结果 void ScatterAB 分发矩阵 AB 中的元素到各个进程中 void MainProcess cannon 算法的主过程 void collectC 收集结果矩阵 C void Mutiply 矩阵相乘 void Printab void Printc int main int argc char argv int i double starttime endtime MPI Init MPI Comm size MPI COMM WORLD MPI Comm rank MPI COMM WORLD if rank 0 printf please input the raw number n fflush stdout scanf d printf n MPI Bcast n atoi argv 1 np n int sqrt p a double malloc np np sizeof double b double malloc np np sizeof double c double malloc np np sizeof double memset c 0 np np sizeof double tempa double malloc np np sizeof double tempb double malloc np np sizeof double if rank 0 在根处理器中为矩阵 ABC 分配空间 A double malloc n sizeof double B double malloc n sizeof double C double malloc n sizeof double for i 0 i n i A i double malloc n sizeof double B i double malloc n sizeof double C i double malloc n sizeof double ProduceABC 在根处理器中随机生成矩阵 AB 初始化矩阵 C ScatterAB 分发矩阵 AB 中的元素到各个进程中 else MPI Recv a np np MPI DOUBLE 0 1 MPI COMM WORLD MPI Recv b np np MPI DOUBLE 0 2 MPI COMM WORLD starttime MPI Wtime 开始时间 MainProcess cannon 算法的主过程 if rank 0 collectC 收集结果矩阵 C PrintABC 输出结果 endtime MPI Wtime printf time used lf n endtime starttime for i 0 i n i free A i free B i free C i free A free B free C else MPI Send c np np MPI DOUBLE 0 1 MPI COMM WORLD free a free b free c free tempa free tempb MPI Finalize return 0 void ProduceABC 在根处理器中生成矩阵 AB int i j srand unsigned time NULL for i 0 i n i for j 0 j n j A i j rand MUL B i j rand MUL C i j 0 0 void PrintABC 输出结果 printf A 0 0 f nB 0 0 f nC 0 0 f n A 0 0 B 0 0 C 0 0 void ScatterAB 分发矩阵 AB 中的元素到各个进程中 int imin imax jmin jmax int sp int i j k m for k 0 k p k 计算相应处理器所分得的矩阵块在总矩阵中的坐标范围 sp int sqrt p imin k sp np imax imin np 1 jmin k sp np jmax jmin np 1 rank 0 的处理器将 A B 中的相应块拷至 tempa tempb 准备向其他处理器发送 m 0 for i imin i imax i for j jmin j jmax j tempa m A i j tempb m B j i 矩阵 B 按列优先存储 m 根处理器将自己对应的矩阵块从 tempa tempb 拷至 a b if k 0 memcpy a tempa np np sizeof double memcpy b tempb np np sizeof double else 根处理器向其他处理器发送 tempa tempb 中相关的矩阵块 MPI Send tempa np np MPI DOUBLE k 1 MPI COMM WORLD MPI Send tempb np np MPI DOUBLE k 2 MPI COMM WORLD void MainProcess cannon 算法的主过程 MPI Comm comm 笛卡尔结构通讯器 int crank int dims 2 periods 2 coords 2 int source dest up down right left int i dims 0 dims 1 int sqrt p periods 0 periods 1 1 MPI Cart create MPI COMM WORLD 2 dims periods 1 MPI Comm rank comm MPI Cart coords comm crank 2 coords MPI Cart shift comm 1 1 MPI Cart shift comm 0 1 MPI Cart shift comm 1 coords 0 MPI Sendrecv replace a np np MPI DOUBLE dest 1 source 1 comm MPI Cart shift comm 0 coords 1 MPI Sendrecv replace b np np MPI DOUBLE dest 1 source 1 comm Mutiply 矩阵相乘 for i 1 i dims 0 i MPI Sendrecv replace a np np MPI DOUBLE left 1 right 1 comm MPI Sendrecv replace b np np MPI DOUBLE up 1 down 1 comm Mutiply 矩阵相乘 MPI Comm free void collectC 收集结果矩阵 C int i j k s m int imin imax jmin jmax int sp int sqrt p 根处理器中的 c 赋给总矩阵 C for i 0 i np i for j 0 j np j C i j c i np j for k 1 k p k 根处理器从其他处理器接收相应的分块 c MPI Recv c np np MPI DOUBLE k 1 MPI COMM WORLD printf rank d n k Printc imin k sp np imax imin np 1 jmin k sp np jmax jmin np 1 将接收到的 c 拷至 C 中的相应位置 从而构造出 C for i imin m 0 i imax i m for j jmin s 0 j jmax j s C i j c m np s void Mutiply 矩阵相乘 int i j k for i 0 i np i for j 0 j np j for k 0 k np k c i np j a i np k b j np k b 按列优先来搞 优化后的行划分算法 带稀疏 include mpi h include include include define DATA int define MUL 3 void magten DATA a int width int i j srand unsigned time NULL for i 0 i width i for j 0 j width j 10 a i width j rand MUL rand 2 rand 3 a i width j 1 rand MUL rand 2 rand 3 a i width j 2 rand MUL rand 2 rand 3 a i width j 3 rand MUL rand 2 rand 3 a i width j 4 rand MUL rand 2 rand 3 a i width j 5 rand MUL rand 2 rand 3 a i width j 6 rand MUL rand 2 rand 3 a i width j 7 rand MUL rand 2 rand 3 a i width j 8 rand MUL rand 2 rand 3 a i width j 9 rand MUL rand 2 rand 3 void matrix c DATA a int width int i j DATA temp for i 0 i width i for j i j width j temp a i width j a i width j a j width i a j width i temp void print matrix DATA a int width int i j for i 0 i width i for j 0 j width j printf d a i width j printf n int main int argc char argv DATA m n p buffer ans int width 1000 int count count1 int myid numprocs MPI Status status int i j k MPI Init MPI Comm rank MPI COMM WORLD MPI Comm size MPI COMM WORLD int line width numprocs if myid 0 m DATA malloc sizeof DATA width width n DATA malloc sizeof DATA width width p DATA malloc sizeof DATA width width buffer DATA malloc sizeof DATA width line ans DATA malloc sizeof DATA width line magten m width magten n width matrix c n width doubleclockstart MPI Wtime for i 1 i numprocs i MPI Send n width width MPI INT i 0 MPI COMM WORLD 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 double time1 MPI Wtime printf send B and A i time 6f n time1 clockstart for i 0 i width i numprocs for j 0 j width j DATA temp 0 for k 0 k width k 10 temp m i width k n j width k temp m i width k 1 n j width k 1 temp m i width k 2 n j width k 2 temp m i width k 3 n j width k 3 temp m i width k 4 n j width k 4 temp m i width k 5 n j width k 5 temp m i width k 6 n j width k 6 temp m i width k 7 n j width k 7 temp m i width k 8 n j width k 8 temp m i width k 9 n j width k 9 p i width j temp double time2 MPI Wtime printf compute the last A i time 6f n time2 time1 for k 1 k numprocs k MPI Recv ans line width MPI INT k 2 MPI COMM WORLD for i 0 i line i count 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 double time3 MPI Wtime printf receive the results from slave time 6f n time3 time2 double clockend MPI Wtime print matrix p width free m free n free p free ans free buffer printf myid d time 6fsecs n myid clockend clockstart else n DATA malloc sizeof DATA width width ans DATA malloc sizeof DATA width line buffer DATA malloc sizeof DATA width line MPI Recv n width width MPI INT 0 0 MPI COMM WORLD MPI Recv buffer width line MPI INT 0 1 MPI COMM WORLD 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