7线性方程组的直接解法 在求解线性方程组( System of Linear Equations)的算法中,有两类最基本的算法,一类是 直接法,即以消去为基础的解法。如果不考虑误差的影响,从理论上讲,它可以在固定步数 内求得方程组的准确解。另一类是迭代解法,它是一个逐步求得近似解的过程,这种方法便 于编制解题程序,但存在着迭代是否收敛及收敛速度快慢的问题。在迭代过程中,由于极限 过程一般不可能进行到底,因此只能得到满足一定精度要求的近似解。本章我们主要介绍几 种直接法,迭代法将在下一章讨论。 1.1高斯消去法解线性方程组 在直接方法中最主要的是高斯消去法( Gaussian elimination)。它分为消元与回代两个过 程,消元过程将方程组化为一个等价的三角方程组,而回代过程则是解这个三角方程组。 1911高斯消去及其串行算法 对于方程组Ax=b,其中A为n阶非奇异阵,其各阶主子行列式不为零,x,b为n维向 量。将向量b看成是A的最后一列,此时A就成了一个n×(n+1)的方程组的增广矩阵 ( Augmented matrⅸx),消去过程实质上是对增广矩阵A进行线性变换,使之三角化。高斯消 去法按k=1,2…;n的顺序,逐次以第k行作为主行进行消去变换,以消去第k列的主元素以 下的元素ak+1k,ak+2k,…,amk。消去过程分为两步,首先计算: k=akak,j=k+1,…,n 这一步称为归一化( Normalization)。它的作用是将主对角线上的元素变成1,同时第k行 上的所有元素与常数向量中的bk都要除以ak。由于A的各阶主子式非零,可以保证在消去 过程中所有主元素aA皆为非零。然后计算 a= ajj-aik aki,i/=k+1,…,n b=b-akbk,i=k+1,…n 这一步称为消元,它的作用是将主对角线ak以下的元素消成0,其它元素与向量B中 的元素也应作相应的变换。 在回代过程中,按下式依次解出xm,xm1,…x1: ①直接解出xn,即xn=bn/anm ②进行回代求x=b-∑ 在归-化的过程中,要用ak作除数,当|ak|很小时,会损失精度,并且可能会导致 商太大而使计算产生溢出。如果系数A虽为非奇异,但不能满足各阶主子式全不为零的条 件,就会出现主元素an为零的情况,导致消去过程无法继续进行。为了避免这种情形,在 每次归一化之前,可增加一个选主元( Pivot)的过程,将绝对值较大的元素交换到主对角线的 置上。根据选取主元的范围不同,选主元的方法主要有列主元与全主元两种。 列主元 Column pivot方法的基本思想是当变换到第k步时,从第k列的aM以下(包括
1. 7 线性方程组的直接解法 在求解线性方程组(System of Linear Equations)的算法中，有两类最基本的算法，一类是 直接法，即以消去为基础的解法。如果不考虑误差的影响，从理论上讲，它可以在固定步数 内求得方程组的准确解。另一类是迭代解法，它是一个逐步求得近似解的过程，这种方法便 于编制解题程序，但存在着迭代是否收敛及收敛速度快慢的问题。在迭代过程中，由于极限 过程一般不可能进行到底，因此只能得到满足一定精度要求的近似解。本章我们主要介绍几 种直接法，迭代法将在下一章讨论。 1.1 高斯消去法解线性方程组 在直接方法中最主要的是高斯消去法(Gaussian Elimination)。它分为消元与回代两个过 程，消元过程将方程组化为一个等价的三角方程组，而回代过程则是解这个三角方程组。 19.1.1 高斯消去及其串行算法 对于方程组 Ax=b，其中 A 为 n 阶非奇异阵，其各阶主子行列式不为零，x，b 为 n 维向 量。将向量 b 看成是 A 的最后一列，此时 A 就成了一个 n×(n+1)的方程组的增广矩阵 (Augmented Matrix)，消去过程实质上是对增广矩阵 A 进行线性变换，使之三角化。高斯消 去法按 k=1,2,…,n 的顺序，逐次以第 k 行作为主行进行消去变换，以消去第 k 列的主元素以 下的元素 ak k ak k ank , , , +1 +2 。消去过程分为两步，首先计算： akj=akj/akk , j=k+1, …,n bk=bk/akk 这一步称为归一化(Normalization)。它的作用是将主对角线上的元素变成 1，同时第 k 行 上的所有元素与常数向量中的 bk都要除以 akk 。由于 A 的各阶主子式非零，可以保证在消去 过程中所有主元素 akk皆为非零。然后计算： aij=aij-aik akj , i,j=k+1, …,n bi= bi -aik bk , i =k+1, …,n 这一步称为消元，它的作用是将主对角线 akk以下的元素消成 0，其它元素与向量 B 中 的元素也应作相应的变换。 在回代过程中，按下式依次解出 xn,xn-1, …,x1： ① 直接解出 xn，即 xn=bn/ann； ② 进行回代求 = + = − = − n j i xi bi aij xj i n 1 , 1, , 2,1 在归−化的过程中，要用 akk 作除数，当∣akk∣很小时，会损失精度，并且可能会导致 商太大而使计算产生溢出。如果系数 A 虽为非奇异，但不能满足各阶主子式全不为零的条 件，就会出现主元素 aii 为零的情况，导致消去过程无法继续进行。为了避免这种情形，在 每次归一化之前，可增加一个选主元(Pivot)的过程，将绝对值较大的元素交换到主对角线的 位置上。根据选取主元的范围不同，选主元的方法主要有列主元与全主元两种。 列主元(Column Pivot)方法的基本思想是当变换到第 k 步时，从第 k 列的 akk以下（包括
ak)的各元素中选出绝对值最大者。然后通过行交换将它交换到ak的位置上。但列主元不 能保证所选的ak是同一行中的绝对值最大者,因此采用了列主元虽然变换过程不会中断 但计算还是不稳定的 全主元方法的基本思想是当变换到第k步时,从右下角(n-k+1)阶子阵中选取绝对值 最大的元素,然后通过行变换和列变换将它交换到ak的位置上。对系数矩阵做行交换不影 响计算结果,但列交换会导致最后结果中未知数的次序混乱。即在交换第i列与第j列后 最后结果中x与x的次序也被交换了。因此在使用全主元的高斯消去法时,必须在选主元 的过程中记住所进行的一切列变换,以便在最后结果中恢复。全主元的高斯消去串行算法如 下,其中a我们用/来表示: 算法19.1单处理器采用全主元高斯消去法的消去过程算法 输入:系数矩阵Anxn,常数向量bnx 输出:解向量xn (I)for i=l to n do shifrin end for (2.1) (2.2)for i=k to n do if( a[i,]>d) then d a[i, jsj,iend if end for (23)ifUs≠k)then (ifor i=l to n de 交换q[k]和a[图 end for (i)交换shjf{和 shifrin 2.4)if(l≠k)then (ifor=k to n do 交换叫k和a[l end for (in)交换bk]和b end if (2.5)forj=k+l to n de ak小=a[k小a[kk (2.6)b=b[]ak],akk]=1 (2.7) fori=k+I to n do a[i小=a[i小-a[k]*a[k nd for (i)b=b[小-a[]*b{k]
akk）的各元素中选出绝对值最大者。然后通过行交换将它交换到 akk的位置上。但列主元不 能保证所选的 akk是同一行中的绝对值最大者，因此采用了列主元虽然变换过程不会中断， 但计算还是不稳定的。 全主元方法的基本思想是当变换到第 k 步时，从右下角（n-k+1）阶子阵中选取绝对值 最大的元素，然后通过行变换和列变换将它交换到 akk的位置上。对系数矩阵做行交换不影 响计算结果，但列交换会导致最后结果中未知数的次序混乱。即在交换第 i 列与第 j 列后， 最后结果中 xi 与 xj 的次序也被交换了。因此在使用全主元的高斯消去法时，必须在选主元 的过程中记住所进行的一切列变换，以便在最后结果中恢复。全主元的高斯消去串行算法如 下，其中 aij 我们用 a[i,j]来表示： 算法 19.1 单处理器采用全主元高斯消去法的消去过程算法 输入：系数矩阵 An×n，常数向量 b n×1 输出：解向量 xn×1 Begin (1)for i=1 to n do shift[i]=i end for (2)for k=1 to n do (2.1) d=0 (2.2)for i=k to n do for j=k to n do if (│a[i,j] │>d) then d=│a[i,j] │， js=j， l=i end if end for end for (2.3) if (js ≠ k) then (i)for i=1 to n do 交换 a[i,k]和 a[i,js] end for (ii)交换 shift[k]和 shift[js] end if (2.4) if (l ≠ k) then (i)forj=k to n do 交换 a[k,j]和 a[l,j] end for (ii)交换 b[k]和 b[l] end if (2.5) for j=k+1 to n do a[k,j]= a[k,j]/a[k,k] end for (2.6) b[k]=b[k]/a[k,k]，a[k,k]=1 (2.7) fori=k+1 to n do (i)forj=k+1 to n do a[i,j]=a[i,j]- a[i,k]* a[k,j] end for (ii)b[i]=b[i]-a[i,k]* b[k]
(il)a[,1 nd for (3)fori= n downto l do/采用全主元高斯消去法的回代过程 forj=汁+ I to n do l]=a[*x可 end for end for (4)for k-l to n do for F=l to n do if(shjl=k)then输出x{k]的值x{ i end if end fo end for 在全主元高斯消去法中,由于每次选择主元素的数据交换情况无法预计,因此我们不考 虑选主元的时间而仅考虑一般高斯消去法的计算时间复杂度。若取一次乘法和加法运算时间 或一次除法运算时间为一个单位时间,则消去过程的时间复杂度为22,回代过程的时间 复杂度为Σi,算法19.1的时间复杂度为(m3+3m2+2n)3=0m2 1912并行高斯消去算法 高斯消去法是利用主行i对其余各行j,>)作初等行变换,各行计算之间没有数据相 关关系,因此可以对矩阵A按行划分。考虑到在计算过程中处理器之间的负载均衡,对A 采用行交叉划分。设处理器个数为p矩阵A的阶数为n,m=「m/p,对矩阵A行交叉划分后, 编号为i(=0,1,…,p-1)的处理器含有A的第,计p,…+(m1)p行和向量B的第i, 计+p,…,计+(m-1)一共m个元素 消去过程的并行是依次以第0,1,…;n-1行作为主行进行消去计算,由于对行的交叉划分 与分布,这实际上是由各处理器轮流选出主行。在每次消去计算前,各处理器并行求其局部 存储器中右下角阶子阵的最大元。若以编号为 my rank的处理器的第i行作为主行,则编号 在 my rank后面的处理器(包括 my rank本身)求其局部存储器中第i行至第m1行元素的 最大元,并记录其行号、列号及所在处理器编号;编号在 my rank前面的处理器求其局部 存储器中第计+1行至第m-1行元素的最大元,并记录其行号、列号及所在处理器编号。然后 通过扩展收集操作将局部存储器中的最大元按处理器编号连接起来并广播给所有处理器,各 处理器以此求得整个矩阵右下角阶子阵的最大元 maxvalue及其所在行号、列号和处理器编 号。若 maxvalue的列号不是原主元素ak的列号,则交换第k列与 maxvalue所在列的两列 数据:若 maxvalue的处理器编号不是原主元素ak的处理器编号,则在处理器间的进行行交 换;若 maxvalue的处理器编号是原主元素ak的处理器编号,但行号不是原主元素ak的行 号,则在处理器内部进行行交换。在消去计算中,首先对主行元素作归一化操作ay=a/ak b=bk,然后将主行广播给所有处理器,各处理器利用接收到的主行元素对其部分行向量 做行变换。若以编号为 my rank的处理器的第i行作为主行,并将它播送给所有的处理器 则编号在 my_ rank前面的处理器(包括 my rank本身)利用主行对其第汁1…,m1行数据 和子向量B做行变换。编号在 my rank后面的处理器利用主行对其第i…,m1行数据和子向 量B做行变换
(iii)a[i,k]=0 end for end for (3)for i=n downto 1 do /*采用全主元高斯消去法的回代过程*/ forj=i+1 to n do b[i]= a[i,j]* x[i] end for end for (4)for k=1 to n do for i=1 to n do if (shift[i]=k) then 输出 x[k]的值 x[i] end if end for end for End 在全主元高斯消去法中，由于每次选择主元素的数据交换情况无法预计，因此我们不考 虑选主元的时间而仅考虑一般高斯消去法的计算时间复杂度。若取一次乘法和加法运算时间 或一次除法运算时间为一个单位时间，则消去过程的时间复杂度为 = n i i 1 2 ，回代过程的时间 复杂度为 = n i i 1 ，算法 19.1 的时间复杂度为(n 3+3n 2+2n)/3=О(n 3 )。 19.1.2 并行高斯消去算法 高斯消去法是利用主行 i 对其余各行 j，(j>i)作初等行变换，各行计算之间没有数据相 关关系，因此可以对矩阵 A 按行划分。考虑到在计算过程中处理器之间的负载均衡，对 A 采用行交叉划分。设处理器个数为 p,矩阵 A 的阶数为 n，m = n / p ，对矩阵 A 行交叉划分后， 编号为 i(i=0,1,…, p -1)的处理器含有 A 的第 i, i+p,…,i+(m-1)p 行和向量 B 的第 i, i+p,…,i+(m-1)p 一共 m 个元素。 消去过程的并行是依次以第 0,1,…,n-1 行作为主行进行消去计算，由于对行的交叉划分 与分布，这实际上是由各处理器轮流选出主行。在每次消去计算前，各处理器并行求其局部 存储器中右下角阶子阵的最大元。若以编号为 my_rank 的处理器的第 i 行作为主行，则编号 在 my_rank 后面的处理器（包括 my_rank 本身）求其局部存储器中第 i 行至第 m-1 行元素的 最大元，并记录其行号、列号及所在处理器编号；编号在 my_rank 前面的处理器求其局部 存储器中第 i+1 行至第 m-1 行元素的最大元，并记录其行号、列号及所在处理器编号。然后 通过扩展收集操作将局部存储器中的最大元按处理器编号连接起来并广播给所有处理器，各 处理器以此求得整个矩阵右下角阶子阵的最大元 maxvalue 及其所在行号、列号和处理器编 号。若 maxvalue 的列号不是原主元素 akk的列号，则交换第 k 列与 maxvalue 所在列的两列 数据；若 maxvalue 的处理器编号不是原主元素 akk的处理器编号，则在处理器间的进行行交 换；若 maxvalue 的处理器编号是原主元素 akk的处理器编号，但行号不是原主元素 akk的行 号，则在处理器内部进行行交换。在消去计算中，首先对主行元素作归一化操作 akj=akj/akk , bk=bk/akk，然后将主行广播给所有处理器，各处理器利用接收到的主行元素对其部分行向量 做行变换。若以编号为 my_rank 的处理器的第 i 行作为主行，并将它播送给所有的处理器。 则编号在 my_rank 前面的处理器（包括 my_rank 本身）利用主行对其第 i+1,…, m-1 行数据 和子向量 B 做行变换。编号在 my_rank 后面的处理器利用主行对其第 i,…,m-1 行数据和子向 量 B 做行变换
回代过程的并行是按xn,xn1,…,x1的顺序由各处理器依次计算x(ip+ my rank),一旦 x(itp+ my rank)被计算出来就立即广播给所有处理器,用于与对应项相乘并做求和计算。具 体算法框架描述如下 算法192全主元高斯消去法过程的并行算法 输入:系数矩阵An×n,常数向量bn×1 输出:解向量xnx1 Begin 对所有处理器 my rank( my rank=0…p-1)同时执行如下的算法 /*消去过程* (I)for i=0 to m-I do for产=0 to p-l do if( my rank<)then/*对于主行前面的块* ()y=p+/*v为主元素的行号* /*确定本处理器所存未消元部分的最大元及其位置存于数组lmax中* (ii)Imax(o=ali+l, v (infor k=i+l to m-I do for t=y to n-l do if( a[k,>Imax[0])then Imax[0]=a[k, n), Imax[l=k, Imax[2]=t, Imax[ 3]=my rank end if end for (1.2f( my rank≥then/*对于主行前面的块* (1)y=p/*y为主元素的行号* *确定本处理器所存未消元部分的最大元及其位置存于数组lmax中* (ii)Imax=a[i, vI infor k=i to m-I do for y to n-l do f(alk.n]I>ImaxO)then Imax[O=a[k, n), max[l=k end fo (1.3)用A! gather操作将每一个处理器中的lmax元素广播到其它所有处理器中 (1.4)/确定最大元及其位置* maxvalue=getpivort(max), maxrow=getrow(max maxcolumn=getcolumn(max), maxrank-getrank(max) (1.5)/*列交换* if( maxcolun≠v)the (ifor t=0 to m de
回代过程的并行是按 xn,xn-1,…,x1 的顺序由各处理器依次计算 x(i*p+my_rank)，一旦 x(i*p+my_rank)被计算出来就立即广播给所有处理器，用于与对应项相乘并做求和计算。具 体算法框架描述如下： 算法 19.2 全主元高斯消去法过程的并行算法 输入：系数矩阵 An×n，常数向量 bn×1 输出：解向量 xn×1 Begin 对所有处理器 my_rank(my_rank=0,…,p-1)同时执行如下的算法: /*消去过程*/ (1)for i=0 to m-1 do for j=0 to p-1 do (1.1)if (my_rank<j) then /*对于主行前面的块*/ (i)v=i*p+j /* v 为主元素的行号*/ /*确定本处理器所存未消元部分的最大元及其位置存于数组 lmax 中*/ (ii)lmax(0)=a[i+1,v] (iii)for k=i+1 to m-1 do for t=v to n-1 do if (│a[k,t]│>lmax[0]) then lmax[0]= a[k,t],lmax[1]=k, lmax[2]=t,lmax[3]=my_rank end if end for end for end if (1.2)if (my_rank ≥ j) then /*对于主行前面的块*/ (i)v=i*p+j /* v 为主元素的行号*/ /*确定本处理器所存未消元部分的最大元及其位置存于数组 lmax 中*/ (ii)lmax[0]=a[i,v] (iii)for k=i to m-1 do for t=v to n-1 do if (│a[k,t]│>lmax[0]) then lmax[0]= a[k,t],lmax[1]=k , lmax[2]=t,lmax[3]=my_rank end if end for end for end if (1.3)用 Allgather 操作将每一个处理器中的 lmax 元素广播到其它所有处理器中 (1.4) /*确定最大元及其位置*/ maxvalue=getpivort(max),maxrow=getrow(max) maxcolumn=getcolumn(max),maxrank=getrank(max) (1.5) /*列交换*/ if (maxcolumn≠v) then (i)for t=0 to m do
交换a,y与a[t, mixcolumn end for (i)交换shif与shif end if (1.6)/*行交换* if(mixcolumn fa[i, v))then (iif([maxrank-i] and [maxon))then nnerexchangerow()/*处理器内部换行 nd if (i)if( maxrank≠then outerexchangerow()/*处理器之间换行” d if end if (1.7if( my rank=j)then/*主行所在的处理器* 对主行作归一化运算* (ifor k-y+l to n-1 do a[i, k]=a[i, kla[i, v] (ib(d=b[ai, vla[i, v]=1 /将主行元素赋予数组f* (infor k-v+I to n-l do fk=ali, k (iv)n=b[ (v)广播主行到所有处理器 else/*非主行所在的处理器* 接收广播来的主行元素存于数组f中 (1.8if( my rank≤then for k=i+l to m-l d (ifor w=v+l to n-l do a[k, w]=alk, w]-w]*alk, v] nd fo (iib(k]=b(k]-(n]*alk, rl end for ifor w=v+l to n-l d alk, w]=alk, w]-fw]* alk, vl nd fe (iib(k]=b(k] n]*akk, vl
交换 a[t,v]与 a[t,maxcolumn] end for (ii)交换 shift[v]与 shift[maxcolumn] end if (1.6) /*行交换*/ if (my_rank=j) then if (maxcolumn ≠ a[i,v]) then (i)if ([maxrank=j] and [i≠maxrow]) then innerexchangerow( ) /*处理器内部换行*/ end if (ii)if (maxrank ≠ j) then outerexchangerow( ) /*处理器之间换行*/ end if end if end if (1.7)if (my_rank=j) then /*主行所在的处理器*/ /*对主行作归一化运算*/ (i)for k=v+1 to n-1 do a[i,k]= a[i,k]/a[i,v] end for (ii)b[i]=b[i]/ a[i,v],a[i,v]=1 /*将主行元素赋予数组 f */ (iii)for k=v+1 to n-1 do f[k]= a[i,k] end for (iv)f[n]=b[i] (v)广播主行到所有处理器 else /*非主行所在的处理器*/ 接收广播来的主行元素存于数组 f 中 end if (1.8)if (my_rank ≤ j) then for k=i+1 to m-1 do (i)for w=v+1 to n-1 do a[k,w]= a[k,w]-f[w]* a[k,v] end for (ii)b[k]=b[k]-f[n]* a[k,v] end for end if (1.9)if (my_rank>j) then for k=i to m-1 do (i)for w=v+1 to n-1 do a[k,w]= a[k,w]-f[w]* a[k,v] end for (ii)b[k]=b[k]-f[n]* a[k,v]