7 线性方程组的直接解法
在求解线性方程组(System of Linear Equations)的算法中，有两类最基本的算法，一类是直接法，即以消去为基础的解法。如果不考虑误差的影响，从理论上讲，它可以在固定步数内求得方程组的准确解。另一类是迭代解法，它是一个逐步求得近似解的过程，这种方法便于编制解题程序，但存在着迭代是否收敛及收敛速度快慢的问题。在迭代过程中，由于极限过程一般不可能进行到底，因此只能得到满足一定精度要求的近似解。本章我们主要介绍几种直接法，迭代法将在下一章讨论。
高斯消去法解线性方程组
在直接方法中最主要的是高斯消去法(Gaussian Elimination)。它分为消元与回代两个过程，消元过程将方程组化为一个等价的三角方程组，而回代过程则是解这个三角方程组。
高斯消去及其串行算法
对于方程组Ax=b，其中A为n阶非奇异阵，其各阶主子行列式不为零，x，b为n维向量。将向量b看成是A的最后一列，此时A就成了一个n×(n+1)的方程组的增广矩阵(Augmented Matrix)，消去过程实质上是对增广矩阵A进行线性变换，使之三角化。高斯消去法按k=1,2,…,n的顺序，逐次以第k行作为主行进行消去变换，以消去第k列的主元素以下的元素。消去过程分为两步，首先计算：
akj=akj/akk , j=k+1, …,n
bk=bk/akk
这一步称为归一化(Normalization)。它的作用是将主对角线上的元素变成1，同时第行上的所有元素与常数向量中的bk都要除以。由于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；
进行回代求
在归?化的过程中，要用akk作除数，当∣akk∣很小时，会损失精度，并且可能会导致商太大而使计算产生溢出。如果系数A虽为非奇异，但不能满足各阶主子式全不为零的条件，就会出现主元素aii为零的情况，导致消去过程无法继续进行。为了避免这种情形，在每次归一化之前，可增加一个选主元(Pivot)的过程，将绝对值较大的元素交换到主对角线的位置上。根据选取主元的范围不同，选主元的方法主要有列主元与全主元两种。
列主元(Column Pivot)方法的基本思想是当变换到第步时，从第列的akk以下（包括akk）的各元素中选出绝对值最大者。然后通过行交换将它交换到akk的位置上。但列主元不能保证所选的akk是同一行中的绝对值最大者，因此采用了列主元虽然变换过程不会中断，但计算还是不稳定的。
全主元方法的基本思想是当变换到第步时，从右下角（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]
(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
在全主元高斯消去法中，由于每次选择主元素的数据交换情况无法预计，因此我们不考虑选主元的时间而仅考虑一般高斯消去法的计算时间复杂度。若取一次乘法和加法运算时间或一次除法运算时间为一个单位时间，则消去过程的时间复杂度为，回代过程的时间复杂度为，算法19.1的时间复杂度为(n3+3n2+2n)/3=О(n3)。
并行高斯消去算法
高斯消去法是利用主行i对其余各行j，(j>i)作初等行变换，各行计算之间没有数据相关关系，因此可以对矩阵A按行划分。考虑到在计算过程中处理器之间的负载均衡，对A采用行交叉划分。设处理器个数为p,矩阵A的阶数为n，，对矩阵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的处理器的第行作为主行，并将它播送给所有的处理器。则编号在my_rank前面的处理器（包括my_rank本身）利用主行对其第i+1,…, m-1行数据和子向量B做行变换。编号在my_rank后面的处理器利用主行对其第i,…,m-1行数据和子向量B做行变换。
回代过程的并行是按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[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]
end for
end if
end for
end for
/*回代过程*/
(2)for i=0 to m-1 do
sum[i]=0.0
end for
(3)for i= m-1 downto 0 do
for j= p-1 downto 0 do
if (my_rank=j) then /*主行所在的处理器*/
(i)x[i*p+j]=(b[i]-sum[i])/a[i,i*p+j]
(ii)将x[i*p+j]广播到所有处理器中
(iii)for k =0 to i-1 do
/*计算有关x[i*p+j]的内积项并累加*/
sum[k]=sum[k]+ a[k,i*p+j]* x[i*p+j]
end for
else /*非主行所在的处理器*/
(iv)接收广播来的x[i*p+j]的值
(v)if (my_rank>j) then
for k =0 to i-1 do
/*计算有关x[i*p+j]的内积项并累加*/
sum[k]=sum[k]+ a[k,i*p+j]* x[i*p+j]
end for
end if
(vi)if (my_rank<j) then
for k =0 to i do
/*计算有关x[i*p+j]的内积项并累加*/
sum[k]=sum[k]+ a[k,i*p+j]* x[i*p+j]
end for
end if
end if
end for
end for
End
消去过程中，参与计算的行向量数在减少，同时参与计算的行向量长度也在减少。设第i次消去参与变换的行向量数为(m-i)，行向量长度为(n-v)，其中v=ip+p/2若取一次乘法和加法运算时间或一次除法运算时间为一个单位时间，则上述算法所需的计算时间为T1= (m-i)*p*(n-v)=(3n2-pn+4mn2)/12，其间共有n行数据作为主行被广播，通信时间为n[(ts+ (n+1)tw) logp；回代过程中，由于0号处理器对对应项进行乘积求和的计算量最大，因此以0号处理器为对象分析。由于0号处理器计算解向量x(0),x(p),…,x((m-1)*p)的时间分别为mp,(m-1)p, …,p，因此其回代过程的计算时间为T2=mp(m+1)/2，解向量x的n个元素被播送给所有处理器的通信时间为n(ts+tw)logp。若不考虑选主元的时间而仅考虑一般高斯消去法的计算时间，则此并行计算时间为Tp= (3n2-pn+4mn2)/12+mp (m+1)/2+nlogp[2ts+ (n+2)tw]。
MPI源程序请参见章末附录。
约当消去法解线性方程组
约当消去及其串行算法
约当消去法(Jordan Elimination)是一种无回代过程的直接解法，它直接将系数矩阵A变换为单位矩阵，经变换后的常数向量b即是方程组的解。这种消去法的基本过程与高斯消去法相同，只是在消去过程中，不但将主对角线以下的元素消成0，而且将主对角线以上的元素也同时消成0。一般约当消去法的计算过程是按k=1,2, …,n的顺序，逐次以第k行作为主行进行消去，以消去第k列除主元素以外的所有元素a1k,a2k, …,ank。记A(0)=A，用约当消去法由A(0)出发通过变换得到方阵序列{ A(k)}，A(k)=[aij(k)]，（k=0,1,2,…,n），每次由A(k－1)到A(k)的过程分为三步：
⑴
⑵
⑶
若取一次乘法和加法运算时间或一次除法运算时间为一个单位时间，则由于在第i次消去中，参与变换的行向量数为n，行向量长度为i，所以下述算法19.3的约当消去法的时间复杂度为=О(n3)。
算法19.3 单处理器上约当消去法求解线性方程组的算法
输入：系数矩阵An×n，常数向量bn×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
endfor
(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) forj=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)for i=1 to n do
if (I ≠ k) then
(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]
(iii)a[i,k]=0
end if
end for
end for
(3)for k=1 to n do
for i=1 to n do
if (shift[i]=k) then 输出x[k]的值b[i] end if
end for
end for
End
约当消去法的并行算法
约当消去法采用与高斯消去法相同的数据划分和选主元的方法。在并行消去过程中，首先对主行元素作归一化操作akj=akj/akk (j=k+1, …,n)，bk= bk/akk ，然后将主行广播给所有处理器，各处理器利用接收到的主行元素对其部分行向量做行变换。若以编号为my_rank的处理器的第i行作为主行，在归一化操作之后，将它广播给所有处理器，则编号不为my_rank的处理器利用主行对其第0,1, …, m-1行数据和子向量做变换，编号为my_rank的处理器利用主行对其除i行以外的数据和子向量做变换（第i个子向量除外）。具体算法框架描述如下：
算法19.4 约当消去法求解线性方程组的并行算法
输入：系数矩阵An×n，常数向量b n×1
输出：解向量xn×1
Begin
对所有处理器my_rank(my_rank=0,…, p-1)同时执行如下的算法:
/*消去过程*/
for i=0 to m-1 do
for j=0 to p-1 do
(1)if (my_rank<j) then /*对于主行前面的块*/
(1.1)v=i*p+j /* v为主元素的行号*/
/*确定本处理器所存的未消元部分的最大元及其位置存于数组lmax中*/
(1.2)lmax(0)=a[i+1,v]
(1.3)for k=i+1 to m-1 do
for t=v to n-1 do
if (│a[k,t]│>lmax[0])
lmax[0]= a[k,t], lmax[1]=k ,
lmax[2]=t ,lmax[3]=my_rank
end if
end for
end for
end if
(2)if (my_rank ≥ j) then
(2.1)v=i*p+j
(2.2)lmax[0]=a[i,v]
(2.3)for k=i to m-1 do
for t=v to n-1 do
if (│a[k,t]│>lmax[0])
lmax[0]= a[k,t], lmax[1]=k,
lmax[2]=t , lmax[3]=my_rank
end if
end for
end for
end if
(3)用Allgather操作将每个处理器中的lmax元素广播到其它所有处理器中
(4)/*确定最大元及其位置*/
maxvalue=getpivort(max),maxrow=getrow(max)
maxcolumn=getcolumn(max),maxrank= getrank(max)
(5)/*列交换*/
if (maxcolumn ≠ v) then
(5.1)for t=0 to m-1 do
交换a[t,v]和a[t,maxcolumn]
end for
(5.2)交换shift[v]和shift[maxcolumn]
end if
(6)/*行交换*/
if (my_rank=j) then
if (maxcolumn ≠ a[i,v]) then
(6.1)if ((maxrank=j) and (i≠maxrow)) then
innerexchangerow( ) /*处理器内部换行*/
end if
(6.2)if (maxrank≠j) then
outerexchangerow( ) /*处理器之间换行*/
end if
end if
end if
(7)if (my_rank=j) then /*主行所在的处理器*/
(7.1)for k=v+1 to n-1 do
a[i,k]= a[i,k]/a[i,v]
end for
(7.2) b[i]=b[i]/a[i,v]，a[i,v]=1
(7.3)for k=v+1 to n-1 do /*将主行元素赋予数组f */
f[k]= a[i,k]
end for
(7.4) f[n]=b[i]
(7.5)广播主行到所有处理器
(7.6)for k=0 to m-1 do /*处理存于该处理器中的非主行元素*/
if (k ≠ i) then
(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 if
end for
else /*非主行所在的处理器*/
(7.7)接收广播来的主行元素存于数组f中
(7.8)for k=0 to m do /*处理非主行元素*/
(i)for w=v+1 to n 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
end for
end for
End
上述算法的计算过程中，参与计算的行向量数为n，行向量长度为(n-v)，其中v=ip+p/2。若取一次乘法和加法运算时间或一次除法运算时间为一个单位时间，则其所需的计算时间为T1= n*(n-v)= n(n-ip-p/2)= mn2-n2/2- n2 (m -1)/2；另外，由于其间共有n行数据依次作为主行被广播，通信时间为n[ts+ (n+1)tw] logp ，所以该算法总共需要的并行计算时间为Tp= mn2-n2/2- n2 (m -1)/2+n[ts+(n+1)tw] log。
MPI源程序请参见所附光盘。
小结
线性代数方程组的求解在科学和工程计算中应用非常广泛，这是因为很多科学和工程问题大多可以化为线性代数方程组来求解。本章主要讨论线性方程组的直接解法，使读者能够了解与掌握利用矩阵变换技巧逐步消元从而求解方程组的基本方法。这种方法可以预先估计运算量，并可以得到问题的准确解，但由于实际计算过程中总存在舍入误差，因此得到的结果并非绝对精确，并且还存在着计算过程的稳定性问题。另外，文献[1]展示了三角形方程组求解器可在多计算机上有效地实现，[2]讨论了共享存储和分布存储结构的并行机上稠密线性方程组的并行算法，[3]是一本很好的综述性专著，它全面地讨论了向量和并行机上线性方程组的直接法和迭代法的并行求解方法，[4]中对各类线性方程组的直接解法有更详尽的讲解与分析，读者可以追踪进一步阅读。
参考文献
Romine C H, Ortega J M. Parallel Solution of Triangular Systems of Equations. Parallel Computing, 1988, 6:109-114
Gallivan K A, Plemmons R J, Sameh A H. Parallel Algorithms for Dense Linear Algebra Computations. SIAM Rev., 1990,32(1):54-135
Ortega J M. Introduction to Parallel and Vector Solution of Linear Systems. Plenum,1988
陈国良 编著．并行算法的设计与分析（修订版）．高等教育出版社，2002.11
附录 全主元高斯消去法并行算法的MPI源程序
源程序 gauss.c
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#include "math.h"
#define a(x,y) a[x*M+y]
#define b(x) b[x]
/*A为M*M的系数矩阵*/
#define A(x,y) A[x*M+y]
#define B(x) B[x]
#define floatsize sizeof(float)
#define intsize sizeof(int)
int M;
int N;
int m;
float *A;
float *B;
int my_rank;
int p;
int l;
MPI_Status status;
void fatal(char *message)
{
printf("%s\n",message);
exit(1);
}
void Environment_Finalize(float *a,float *b,float
*x,float *f)
{
free(a);
free(b);
free(x);
free(f);
}
int main(int argc, char **argv)
{
int i,j,t,k,my_rank,group_size;
int i1,i2;
int v,w;
float temp;
int tem;
float *sum;
float *f;
float lmax;
float *a;
float *b;
float *x;
int *shift;
FILE *fdA,*fdB;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,
&group_size);
MPI_Comm_rank(MPI_COMM_WORLD,
&my_rank);
p=group_size;
if (my_rank==0)
{
fdA=fopen("dataIn.txt","r");
fscanf(fdA,"%d %d", &M, &N);
if (M != N-1)
{
printf("the input is wrong\n");
exit(1);
}
A=(float *)malloc(floatsize*M*M);
B=(float *)malloc(floatsize*M);
for(i = 0; i < M; i++)
{
for(j = 0; j < M; j++)
fscanf(fdA,"%f", A+i*M+j);
fscanf(fdA,"%f", B+i);
}
fclose(fdA);
}
/*0号处理器将M广播给所有处理器*/
MPI_Bcast(&M,1,MPI_INT,0,
MPI_COMM_WORLD);
m=M/p;
if (M%p!=0) m++;
/*各处理器为主行元素建立发送和接收缓冲区(M+1) */
f=(float*)malloc(sizeof(float)*(M+1));
/*分配至各处理器的子矩阵大小为m*M*/
a=(float*)malloc(sizeof(float)*m*M);
/*分配至各处理器的子向量大小为m*/
b=(float*)malloc(sizeof(float)*m);
sum=(float*)malloc(sizeof(float)*m);
x=(float*)malloc(sizeof(float)*M);
shift=(int*)malloc(sizeof(int)*M);
if (a==NULL||b==NULL||f==NULL||
sum==NULL||x==NULL||shift==NULL)
fatal("allocate error\n");
for(i=0;i<M;i++)
shift[i]=i;
/*0号处理器采用行交叉划分将矩阵A划分为
大小为m*M的p块子矩阵,将B划分为大小为m的p块子向量，依次发送给1至p-1号处理器*/
if (my_rank==0)
{
for(i=0;i<m;i++)
for(j=0;j<M;j++)
a(i,j)=A(i*p,j);
for(i=0;i<m;i++)
b(i)=B(i*p);
for(i=0;i<M;i++)
if ((i%p)!=0)
{
i1=i%p;
i2=i/p+1;
MPI_Send(&A(i,0),M,
MPI_FLOAT,i1,i2,
MPI_COMM_WORLD);
MPI_Send(&B(i),1,
MPI_FLOAT,i1,i2,
MPI_COMM_WORLD);
}
}
else
{
for(i=0;i<m;i++)
{
MPI_Recv(&a(i,0),M,MPI_FLOAT,
0,i+1,MPI_COMM_WORLD,
&status);
MPI_Recv(&b(i),1,MPI_FLOAT,
0,i+1,MPI_COMM_WORLD,
&status);
}
}
/*消去*/
for(i=0;i<m;i++)
for(j=0;j<p;j++)
{
/*j号处理器负责广播主行元素*/
if (my_rank==j)
{
/*主元素在原系数矩阵A中
的行号和列号为v*/
v=i*p+j;
lmax=a(i,v);
l=v;
/*在同行的元素中找最大元，
并确定最大元所在的列l*/
for(k=v+1;k<M;k++)
if (fabs(a(i,k))>lmax)
{
lmax=a(i,k);
l=k;
}
/*列交换*/
if (l!=v)
{
for(t=0;t<m;t++)
{
temp=a(t,v);
a(t,v)=a(t,l);
a(t,l)=temp;
}
tem=shift[v];
shift[v]=shift[l];
shift[l]=tem;
}
/*归一化*/
for(k=v+1;k<M;k++)
a(i,k)=a(i,k)/a(i,v);
b(i)=b(i)/a(i,v);
a(i,v)=1;
for(k=v+1;k<M;k++)
f[k]=a(i,k);
f[M]=b(i);
/*发送归一化后的主行*/
MPI_Bcast(&f[0],M+1,
MPI_FLOAT,my_rank,MPI_COMM_WORLD);
/*发送主行中主元素所在的
列号*/
MPI_Bcast(&l,1,MPI_INT,
my_rank,
MPI_COMM_WORLD);
}
else
{
v=i*p+j;
/*接收归一化后的主行*/
MPI_Bcast(&f[0],M+1,
MPI_FLOAT, MPI_COMM_WORLD);
/*接收主行中主元素所在的
列号*/
MPI_Bcast(&l,1,MPI_INT, j,
MPI_COMM_WORLD);
/*列交换*/
if (l!=v)
{
for(t=0;t<m;t++)
{
temp=a(t,v);
a(t,v)=a(t,l);
a(t,l)=temp;
}
tem=shift[v];
shift[v]=shift[l];
shift[l]=tem;
}
}
/*编号小于j的处理器（包括j本身）
利用主行对其第i+1,…,m-1行数据做行变换*/
if (my_rank<=j)
for(k=i+1;k<m;k++)
{
for(w=v+1;w<M;w++)
a(k,w)=a(k,w)
-f[w]* a(k,v);
b(k)=b(k)-f[M]*a(k,v);
}
/*编号大于j的处理器利用主行对
其第i,…,m-1行数据做行变换*/
if (my_rank>j)
for(k=i;k<m;k++)
{
for(w=v+1;w<M;w++)
a(k,w)=a(k,w)
-f[w]*a(k,v);
b(k)=b(k)-f[M]*a(k,v);
}
}
for(i=0;i<m;i++)
sum[i]=0.0;
/*回代过程*/
for(i=m-1;i>=0;i--)
for(j=p-1;j>=0;j--)
if (my_rank==j)
{
x[i*p+j]=(b(i)-sum[i])/
a(i,i*p+j);
MPI_Bcast(&x[i*p+j],1,
MPI_FLOAT,my_rank,
MPI_COMM_WORLD);
for(k=0;k<i;k++)
sum[k]=sum[k]+
a(k,i*p+j)*
x[i*p+j];
}
else
{
MPI_Bcast(&x[i*p+j],1,
MPI_FLOAT,j,
MPI_COMM_WORLD);
if (my_rank>j)
for(k=0;k<i;k++)
sum[k]=sum[k]
+a(k,i*p+j)
*x[i*p+j];
if (my_rank<j)
for(k=0;k<=i;k++)
sum[k]=sum[k]
+a(k,i*p+j)
*x[i*p+j];
}
if (my_rank!=0)
for(i=0;i<m;i++)
MPI_Send(&x[i*p+my_rank],1,
MPI_FLOAT,0,i,
MPI_COMM_WORLD);
else
/*0号处理器从其余各处理器中接收子
解向量x*/
for(i=1;i<p;i++)
for(j=0;j<m;j++)
MPI_Recv(&x[j*p+i],1,
MPI_FLOAT, i,j,
MPI_COMM_WORLD,
&status);
if (my_rank==0)
{
printf("Input of file \"dataIn.txt\"\n");
printf("%d\t%d\n", M, N);
for(i=0;i<M;i++)
{
for(j=0;j<M;j++)
printf("%f\t",A(i,j));
printf("%f\n",B(i));
}
printf("\nOutput of solution\n");
for(k=0;k<M;k++)
{
for(i=0;i<M;i++)
{
if (shift[i]==k)
printf("x[%d]=%f\n",
k,x[i]);
}
}
}
MPI_Finalize();
Environment_Finalize(a,b,x,f);
return(0);
}
2. 运行实例
编译：mpicc –o gauss gauss.cc
运行：可以使用命令 mpirun –np ProcessSize gauss来运行该程序，其中ProcessSize是所使用的处理器个数,这里取为5。本实例中使用了
mpirun –np 5 gauss
运行结果：
Input of file "dataIn.txt"
4 5
1.000000 4.000000 -2.000000 3.000000 6.000000
2.000000 2.000000 0.000000 4.000000 2.000000
3.000000 0.000000 -1.000000 2.000000 1.000000
1.000000 2.000000 2.000000 -3.000000 8.000000
Output of solution
x[0]=1.000000
x[1]=2.000000
x[2]=0.000000
x[3]=-1.000000
说明：该运行实例中，A为4×4的矩阵，B是长度为4的向量，它们的值存放于文档“dataIn.txt”中，其中前4列为矩阵A，最后一列为向量B，最后输出线性方程组AX=B的解向量X。
