(n-1)(t,+mn)ogp= O(n log p),因此并行计算时间Tp=(n2-n)/3p+(n-1)(t+nhn)ogp MPI源程序请参见章末附录。 16QR分解 A=a]为一个n阶实矩阵,对A进行QR分解,就是求一个非奇异( Nonsingular)方阵Q 与上三角方阵R,使得A=QR。其中方阵Q满足:Q=g,称为正交矩阵 Orthogonal matrix) 因此QR分解又称为正交三角分解 16.1矩阵QR分解的串行算法 对A进行正交三角分解,可用一系列平面旋转矩阵左乘A,以使A的下三角元素逐个地 被消为0。设要消去的下三角元素为a(1>),则所用的平面旋转方阵为 C S 1 其中,除了(,0(/元素等于c,(u/元素与(,元素分别等于-s与s以外,其余元素皆与单位 方阵Ⅰ的对应元素相等,而c,s由下式计算: 其中,O为旋转角度。这样,在旋转后所得到的方阵A中,元素an'为0,A与A相比仅在 第i行、第j行上的元素不同: k=C×a1k+S×a 消去A的r=m(n-1)2个下三角元素后得到上三角阵R,实际上是用r个旋转方阵Q,Q2…,Q 左乘A,即 QQ-1…Q1A 设Qo=QQ-1…Q,易知Q为一正交方阵,则有 R=O04 Qo R OoR=OR 其中Q=Q01=Q0为一正交方阵,而Qo可以通过对单位阵进行同样的变换而得到这样可 得到A的正交三角分解。QR分解串行算法如下: 算法1811单处理器上矩阵的QR分解串行算法算法 输入:矩阵Anxm,单位矩阵Q 输出:矩阵Qnxm,矩阵Rnx
(n-1)(ts+ntw)logp=O(n log p)，因此并行计算时间 Tp=(n 3 -n)/3p+(n-1)(ts+ntw)logp。 MPI 源程序请参见章末附录。 1.6 QR 分 解 A=[aij] 为一个 n 阶实矩阵，对 A 进行 QR 分解，就是求一个非奇异(Nonsingular)方阵 Q 与上三角方阵 R，使得 A=QR。其中方阵 Q 满足：QT=Q−1，称为正交矩阵(Orthogonal Matrix)， 因此 QR 分解又称为正交三角分解。 1.6.1 矩阵 QR 分解的串行算法 对 A 进行正交三角分解，可用一系列平面旋转矩阵左乘 A，以使 A 的下三角元素逐个地 被消为 0。设要消去的下三角元素为 aij(i>j)，则所用的平面旋转方阵为： 其中，除了(i,i)(j,j)元素等于 c，(i,j)元素与(j,i)元素分别等于−s 与 s 以外，其余元素皆与单位 方阵 I 的对应元素相等，而 c,s 由下式计算： 2 2 cos / a jj a jj aij c = = + 2 2 sin / aij a jj aij s = = + 其中，θ 为旋转角度。这样，在旋转后所得到的方阵 A' 中，元素 ' aij 为 0，A' 与 A 相比仅在 第 i 行、第 j 行上的元素不同： jk jk aik a = c a + s ' ik jk aik a = −s a + c ' (k =1,2,..., n) 消去 A 的 r=n(n-1)/2 个下三角元素后得到上三角阵 R，实际上是用 r 个旋转方阵 Q1,Q2,…,Qr 左乘 A，即： R= Qr Qr-1…Q1A 设 Q0=Qr Qr-1…Q1，易知 Q0 为一正交方阵，则有： R= Q0A 即 A = Q0 −1R= Q0 TR= QR 其中 Q= Q0 −1 =Q0 T为一正交方阵，而 Q0 可以通过对单位阵 I 进行同样的变换而得到, 这样可 得到 A 的正交三角分解。QR 分解串行算法如下： 算法 18.11 单处理器上矩阵的 QR 分解串行算法算法 输入：矩阵 An×n，单位矩阵 Q 输出：矩阵 Qn×n，矩阵 Rn×n − = 1 1 1 1 s c c s Qij
B (I)for j=l to n do for i=j+I to n do (i)sq=sqrt(a*al Halii*aliD C[//q (iifor k-l to n do qk=c*叫Uk]+s*au[i,] gkk]=c *ql, k]+s*qli, k ailk=-s*al, k]+c*ali, k qik]=-s*ql, k]+c*qli, k (infor k=l to n do al, k=ailk qU, k=qilk a[i, k]=ai[k] qli k]=qi[) end fo end for (2)R=A (3 MATRIX TRANSPOSITION(Q)/*对Q实施算法18.1的矩阵转置* 算法18.11要对m(n-1)2个下三角元素进行消去,每消去一个元素要对两行元素进行旋 转变换,而对一个元素进行旋转变换又需要4次乘法和2次加法,所以若取一次乘法或一次 加法运算时间为一个单位时间,则该算法的计算时间为m(n-1)/2*2n*6=6r2(n-1)=Om3) 16.2矩阵QR分解的并行算法 由于QR分解中消去a时,同时要改变第i行及第j行两行的元素,而在LU分解中, 仅利用主行变更第j行的元素。因此QR分解并行计算中对数据的划分与分布与LU分 解就不一样。设处理器个数为p,对矩阵A按行划分为p块,每块含有连续的m行向量, m=「n/p,这些行块依次记为AA,…Ap,分别存放在标号为0,1…p1的处理器中 在0号处理器中,计算开始时,以第0行数据作为主行,依次和第0,1,…;m-1行数据做 旋转变换,计算完毕将第0行数据发送给1号处理器,以使1号处理器将收到的第0行数据 作为主行和自己的m行数据依次做旋转变换;在此同时,0号处理器进一步以第1行数据作 为主行,依次和第2,3,…m1行数据做旋转变换,计算完成后将第1行数据发送给1号处理 器:如此循环下去。一直到以第m-1行数据作为主行参与旋转变换为止。 在1号处理器中,首先以收到的0号处理器的第0行数据作为主行和自己的m行数据 做旋转变换,计算完将该主行数据发送给2号处理器:然后以收到的0号处理器的第1行数 据作为主行和自己的m行数据做旋转变换,再将该主行数据发送给2号处理器:如此循环 下去,一直到以收到的0号处理器的第m-1行数据作为主行和自己的m行数据做旋转变换 并将第m-1行数据发送给2号处理器为止。然后,1号处理器以自己的第0行数据作为主行 依次和第1,2…,胙-1行数据做旋转变换,计算完毕将第0行数据发送给2号处理器:接着以
Begin (1)for j=1 to n do for i=j+1 to n do (i)sq=sqrt(a[j,j]*a[j,j]+a[i,j]*a[i,j]) c= a[j,j]/sq s= a[i,j]/sq (ii)for k=1 to n do aj[k]= c*a[j,k] + s*a[i,k] qj[k]=c*q[j,k] + s*q[i,k] ai[k]= −s*a[j,k] + c*a[i,k] qi[k]= −s*q[j,k] + c*q[i,k] end for (iii)for k=1 to n do a[j,k]=aj[k] q[j,k]=qj[k] a[i,k]=ai[k] q[i,k]=qi[k] end for end for end for (2)R=A (3)MATRIX_TRANSPOSITION(Q) /* 对 Q 实施算法 18.1 的矩阵转置*/ End 算法 18.11 要对 n(n-1)/2 个下三角元素进行消去，每消去一个元素要对两行元素进行旋 转变换，而对一个元素进行旋转变换又需要 4 次乘法和 2 次加法，所以若取一次乘法或一次 加法运算时间为一个单位时间，则该算法的计算时间为 n(n-1)/2*2n*6=6n 2 (n-1)=O(n 3 )。 1.6.2 矩阵 QR 分解的并行算法 由于 QR 分解中消去 aij 时，同时要改变第 i 行及第 j 行两行的元素，而在 LU 分解中， 仅利用主行 i(i<j)变更第 j 行的元素。因此 QR 分解并行计算中对数据的划分与分布与 LU 分 解就不一样。设处理器个数为 p，对矩阵 A 按行划分为 p 块，每块含有连续的 m 行向量， m = n / p ，这些行块依次记为 A0,A1, …,Ap-1，分别存放在标号为 0,1,…,p-1 的处理器中。 在 0 号处理器中，计算开始时，以第 0 行数据作为主行，依次和第 0,1,…,m-1 行数据做 旋转变换，计算完毕将第 0 行数据发送给 1 号处理器，以使 1 号处理器将收到的第 0 行数据 作为主行和自己的 m 行数据依次做旋转变换；在此同时，0 号处理器进一步以第 1 行数据作 为主行，依次和第 2,3,…,m-1 行数据做旋转变换，计算完成后将第 1 行数据发送给 1 号处理 器；如此循环下去。一直到以第 m-1 行数据作为主行参与旋转变换为止。 在 1 号处理器中，首先以收到的 0 号处理器的第 0 行数据作为主行和自己的 m 行数据 做旋转变换，计算完将该主行数据发送给 2 号处理器；然后以收到的 0 号处理器的第 1 行数 据作为主行和自己的 m 行数据做旋转变换，再将该主行数据发送给 2 号处理器；如此循环 下去，一直到以收到的 0 号处理器的第 m-1 行数据作为主行和自己的 m 行数据做旋转变换 并将第 m-1 行数据发送给 2 号处理器为止。然后，1 号处理器以自己的第 0 行数据作为主行， 依次和第 1, 2,…,m-1 行数据做旋转变换，计算完毕将第 0 行数据发送给 2 号处理器；接着以