[计算机软件及应用]矩阵运算.docVIP

矩阵运算 矩阵运算是数值计算中最重要的一类运算,特别是在线性代数和数值分析中,它是一种最基本的运算。讨论的矩阵运算包括矩阵转置、矩阵向量相乘、矩阵乘法、矩阵分解以及方阵求逆等。在讨论并行矩阵算法时分三步进行：①算法描述及其串行算法；②算法的并行化及其实现算法框架以及简单的算法分析；③算法实现的MPI源程序，以利于读者实践操作。 矩阵转置 矩阵转置及其串行算法 对于一个n阶方阵A=[aij]，将其每一下三角元素aij (ij)沿主对角线与其对称元素aji互换就构成了转置矩阵AT。假设一对数据的交换时间为一个单位时间，则下述矩阵转置(Matrix Transposing)算法18.1的运行时间为(n2-n)/2=O(n2)。 算法18.1 单处理器上矩阵转置算法 输入：矩阵An×n 输出：矩阵An×n的转置ATn×n Begin for i=2 to n do for j=1 to i-1 do 交换a[i,j]和a[j,i] endfor endfor End 图 18.1 子块转置 矩阵转置并行算法 此处主要讨论网孔上的块棋盘划分(Block-Checker Board Partitioning，又称为块状划分)的矩阵转置算法，它对循环棋盘划分(Cyclic-Checker Board Partitioning)也同样适用。另外，超立方上块棋盘划分的矩阵转置算法可参见文献[1]。 实现矩阵的转置时，若处理器个数为p，且它们的编号依次是0,1, …,p-1，则将n阶矩阵A分成p个大小为m×m的子块，。p个子块组成一个的子块阵列。记其中第i行第j列的子块为Aij，它含有A的第(i-1)m+1至第im行中的第(j-1)m+1至第jm列的所有元素。对每一处理器按行主方式赋以二维下标，记编号为i的处理器的二维下标为(v,u)，其中，，将A的子块存入下标为(v,u)表示的对应处理器中。 这样，转置过程分两步进行：第一步，子块转置，具体过程如图 18.1所示；第二步，处理器内部局部转置。 为了避免对应子块交换数据时处理器发生死锁，可令下三角子块先向与之对应的上三角子块发送数据，然后从上三角子块接收数据；上三角子块先将数据存放在缓冲区buffer中，然后从与之对应的下三角子块接收数据；最后再将缓冲区中的数据发送给下三角子块。具体并行算法框架描述如下： 算法18.2 网孔上的矩阵转置算法 输入：矩阵An×n 输出：矩阵An×n的转置ATn×n Begin 对所有处理器my_rank(my_rank=0,…,p-1)同时执行如下的算法: (1)计算子块的行号 v=my_rank/ sqrt(p), 计算子块的列号u=my_rank mod sqrt(p) (2)if (uv) then /*对存放下三角块的处理器*/ (2.1)将所存的子块发送到其对角块所在的处理器中 (2.2)接收其对角块所在的处理器中发来的子块 else /*对存放上三角块的处理器*/ (2.3)将所存的子块在缓冲区buffer中做备份 (2.4)接收其对角块所在的处理器中发来的子块 (2.5)将buffer中所存的子块发送到其对角块所在的处理器中 end if (3)for i=1 to m do /*处理器内部局部转置*/ for j=1 to i do 交换a[i,j]和a[j,i] end for end for End 若记ts为发送启动时间, tw为单位数据传输时间，th为处理器间的延迟时间，则第一步由于每个子块有n2/p个元素，又由于通信过程中为了避免死锁，错开下三角子块与上三角子块的发送顺序，因此子块的交换时间为；第二步，假定一对数据的交换时间为一个单位时间，则局部转置时间为。因此所需的并行计算时间。 MPI源程序请参见所附光盘。 矩阵-向量乘法 矩阵-向量乘法及其串行算法 矩阵-向量乘法(Matrix-Vector Multiplication)是将一个n×n阶方阵A=[aij]乘以n×1的向量B=[b1,b2, …,bn]T得到一个具有n个元素的列向量C=[c1,c2, …,cn]T。假设一次乘法和加法运算时间为一个单位时间，则下述矩阵向量乘法算法18.3的时间复杂度为O(n2)。 算法18.3 单处理器上矩阵-向量乘算法 输入：An×n，Bn×1 输出：Cn×1 Begin for i=0 to n-1 do c[i]= 0 for j=0 to n-1 do c[i]=c[i] + a[i,j]*b[j] end for end for End 矩阵-向量乘法的并行算法 矩阵-向量乘法同样可以有带状划分(Striped Partitioning)和棋盘划分(Checker Board Part