WORD完美格式
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?1k,ak?2k,?,ank。消去过程分为两步，首先计算：
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；
② 进行回代求xi?bi??aijxj,j?i?1ni?n?1,?,2,1
在归?化的过程中，要用akk作除数，当∣akk∣很小时，会损失精度，并且可能会导致商太大而使计算产生溢出。如果系数A虽为非奇异，但不能满足各阶主子式全不为零的条件，就会出现主元素aii为零的情况，导致消去过程无法继续进行。为了避免这种情形，在每次归一化之前，可增加一个选主元(Pivot)的过程，将绝对值较大的元素交换到主对角线的位置上。根据选取主元的范围不同，选主元的方法主要有列主元与全主元两种。
列主元(Column Pivot)方法的基本思想是当变换到第k步时，从第k列的akk以下（包括
..整理分享..
WORD完美格式
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 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]
..整理分享..
end if WORD完美格式
(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?1n2?ii?1，回代过程的时间
19.1的时间复杂度为(n3+3n2+2n)/3=О(n3)。
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做行变换。
..整理分享..