高性能矩阵乘法

1、高性能矩阵乘法第1页，共27页，2022年，5月20日，20点24分，星期四2022/9/82并行算法优化研究相对于传统面向对象串行算法的4个挑战：同步：两个或者多个线程协调其行为的过程通信：与线程之间交换数据相关的带宽和延迟问题负载均衡：多个线程之间工作量分布的情况，给各个线程（执行核）分配均匀的工作可扩展性：衡量在性能更加强劲的系统上运行软件时能否有效利用更多线程的指标， 观察应用程序在更高级的平台上运行4核到8核线性增长第2页，共27页，2022年，5月20日，20点24分，星期四2022/9/83多线程（核）设计主要分解模式任务分解：对程序根据其执行的功能进行分解的过程数据分解：将应用

2、程序根据各任务所处理的数据而非按任务的天然特性来进行分解数据流分解：研究数据在诸任务之间如何流动，根据任务之间的数据流关系对问题进行分解模式分解方式任务级并行模式任务分解Divide and Conquer任务/数据分解几何分解模式数据分解流水线模式数据流分解波峰（wavefront）模式数据流分解第3页，共27页，2022年，5月20日，20点24分，星期四2022/9/84多线程（核）设计主要分解模式任务分解：对程序根据其执行的功能进行分解的过程数据分解：将应用程序根据各任务所处理的数据而非按任务的天然特性来进行分解数据流分解：研究数据在诸任务之间如何流动，根据任务之间的数据流关系对问题进

3、行分解分解方式设计说明任务分解不同的程序行为采用不同的线程实现常用于GUI应用程序数据分解多个线程对不同的数据块执行相同的操作常用于音频、图像处理和科学计算应用程序数据流分解一个线程的输出作为另一个线程的输入尤其应注意尽量消除启动和排空延迟第4页，共27页，2022年，5月20日，20点24分，星期四2022/9/85矩阵乘法算法探讨在工程科学计算中，矩阵乘积是最基本的运算典型的n阶稠密方阵乘积算法的时间复杂度是O(n3) 。目前对大型矩阵乘积运算的处理主要是采用分治思想，将矩阵分布在多个节点上，但每个结点上的小矩阵仍要立方级乘法次数。 基于分之思想的两种划分策略：条形划分和块状（棋盘）划分的

4、6种常见分布式矩阵乘法并行算法。第5页，共27页，2022年，5月20日，20点24分，星期四2022/9/86基于不同划分策略的矩阵乘法算法探讨1、条形（striped partitioning）划分的矩阵乘法并行算法 行条划分列条划分两两组合：行列、行行、列列、列行 第6页，共27页，2022年，5月20日，20点24分，星期四2022/9/87基于不同划分策略的矩阵乘法算法探讨2、块状划分(checkerboard partitioning)的矩阵乘法并行算法 称为棋盘划分第7页，共27页，2022年，5月20日，20点24分，星期四Cannon Description for impl

5、ementation of MPI program to compute Matrix Matrix Multiplication using block checkerboard partitioning and Cannon Algorithm 2022/9/88第8页，共27页，2022年，5月20日，20点24分，星期四Cannon Objective Computing the matrix-matrix multiplication on SMP System. Use block checkerboard partitioning of the matrices and Cann

6、ons Algorithm. AssumptionSize of the square matrices p= q2 and the size of square matrices A and B is evenly divisible by q. It is assumed that the number of blocks are equal to the number of processors.2022/9/89第9页，共27页，2022年，5月20日，20点24分，星期四Cannon Cannons algorithm is based on cartesian virtual to

7、pologyA and B are square matrices of size n and C be the output matrix.These matrices are dived into blocks or submatrices to perform matrix-matrix operations in paralleln x n matrix A can be regarded as q x q array of blocks Ai, j (0=i q, 0= j q) such that each block is an (n/q) x (n/q) submatrixWe

8、 use p processors to implement the block version of matrix multiplication in parallel by choosing q as a square root of p and compute a distinct block Ci, j on each processor.2022/9/810第10页，共27页，2022年，5月20日，20点24分，星期四传统并行 2022/9/811第11页，共27页，2022年，5月20日，20点24分，星期四传统并行 The matrices A and B are partit

9、ioned into p blocks, A i, j and B i, j (0=i q, 0= j q) of size (n/q x n/q) on each process. These blocks are mapped onto a q x q logical mesh of processes. The processes are labeled from P0,0 to Pq-1,q-1.2022/9/812第12页，共27页，2022年，5月20日，20点24分，星期四传统并行 Process Pi, j initially store block matrices Ai,

10、j and Bi, j and computes block Ci, j of result matrix. To compute submatrix Ci, j , we need all submatrices , Ai, k and Bk, j ( 0 = k q ) . To acquire all the required blocks, an all-to-all broadcast of matrix Ai, j s is performed in each row and similarly in each column of matrix Bi, js. MPI collec

11、tive communication is used to perform this operations. 2022/9/813第13页，共27页，2022年，5月20日，20点24分，星期四传统并行 After Pi, j acquires, A i,0 , A i,1 , A i,2 , A i, q-1 and B0, j , B1, j , B2, j , Bq-1, j , it performs the serial block matrix to matrix multiplication and accumulates the partial block matrix Ci,

12、 j of matrix C . To obtain the resultant product matrix C, processes with rank 0 gathers all the block matrices by using MPI_Gather collective communication operation.2022/9/814第14页，共27页，2022年，5月20日，20点24分，星期四Cannon p processors arranged in q x q square grid of processors and the input matrices.A an

13、d B are distributed among the processes in checkerboard fashion. It results in constructing p block matrices of A and B. It uses only point-to-point communication for circularly shifting blocks of matrix A and matrix B among p processes.2022/9/815第15页，共27页，2022年，5月20日，20点24分，星期四Cannon-inital The alg

14、orithm proceeds in q stages. The first step in this algorithm is to perform initial alignment of the block matrix A and block matrix B. The blocks of matrix A are circularly shifted to the i positions to left in the row of the square grid of processes, where i is the row number of the process in the

15、 mesh. The blocks of matrix B are circularly shifted j positions upwards, where j is the column number of the process in the processes mesh.2022/9/816第16页，共27页，2022年，5月20日，20点24分，星期四Cannon-inital2022/9/817第17页，共27页，2022年，5月20日，20点24分，星期四Cannon-runningThe algorithm performs the following steps in eac

16、h stage : 1. Multiply the block of matrix A and matrix B and add the resultant matrix to get the block matrix C, which is initially set to zero.2. Circularly shift the blocks of matrix A to left in the rows of the processes and the blocks of matrix B upwards in the columns of the square grid of proc

17、esses in a wrap around manner.2022/9/818第18页，共27页，2022年，5月20日，20点24分，星期四Cannon-running2022/9/819第19页，共27页，2022年，5月20日，20点24分，星期四Cannon-running2022/9/820第20页，共27页，2022年，5月20日，20点24分，星期四书中Cannon-bug2022/9/821MPI_Send and MPI_Recv is not used for point-to-point communicationbecause if all the processes

18、 call MPI_Send or MPI_Recv in different order the deadlocked situation may arise .How to fix?指派一个缓冲区，使用MPI_Irecv/MPI_Isend 非阻塞式通讯函数，MPI_wait.MPI_Sendrecv. 第21页，共27页，2022年，5月20日，20点24分，星期四2022/9/822Cannon-bug死锁的问题问题来源于main_shift()这个函数中MPI函数的使用。在Cannon-mpi代码的main_shift()模块中，文献中算法使用的是MPI的阻塞通信函数:MPI_Sen

19、d/MPI_Recv,这就使得Cannon算法在执行循环左移和循环上移时，矩阵规模超过共享buff的容量时出现循环等待的死锁状况。在曙光4000集群系统上，该算法的发生死锁的矩阵下限规模是200200的浮点型矩阵。第22页，共27页，2022年，5月20日，20点24分，星期四2022/9/823Cannon-bug原始（阻塞式）的main_shift模块： void main_shift() /* 将分块b左移位*/ MPI_Send(a , dl2, MPI_FLOAT, get_index(my_row, my_col-1, sp), 1, MPI_COMM_WORLD); MPI_Re

20、cv(a , dl2, MPI_FLOAT, get_index(my_row, my_col+1, sp), 1, MPI_COMM_WORLD, &status); /* 将分块b上移位*/ MPI_Send(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD); MPI_Recv(b , dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status); 第23页，共27页，2022年，5月20日，20点24分，

21、星期四2022/9/824Cannon-bug改进（非阻塞式）的main_shift模块 ci*dl+j += ai*dl+k*bj*dl+k;/改进了的Cannon按行存取 /* 将分块a左移位*/ MPI_Isend(a , dl2, MPI_FLOAT, get_index(my_row, my_col-1, sp), 1, MPI_COMM_WORLD, &myrequest_s); MPI_Irecv(buf , dl2, MPI_FLOAT, get_index(my_row, my_col+1, sp),1,MPI_COMM_WORLD,&myrequest_r );MPI_Wa

22、it(&myrequest_s,&status);MPI_Wait(&myrequest_r,&status);memcpy(a,buf,sizeof(float) *dl2); /* 将分块b上移位*/ MPI_Isend(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD,& myrequest_s); MPI_Irecv(buf , dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD,&myrequest_r);MPI_Wait(&myrequest_s,&status);MPI_Wait(&myrequest_r,&status);memcpy(b,buf,sizeof(float) *dl2);第24页，共27页，2022年，5月20日，20点24分，星期四2022/9/825Cannon-bugMPI_Irecv仅仅初始化接受操作,在与之对应的MPI_Wait函数的调用返回之前，将不能访问bufferMPI_Irecv函数返回时，handle指向一个MPI_Request对象，它代表了一个已近初始化了的通信操作。这个函数并不返回一个指向MPI_St