Transcript E:WHU教学2010Lecture2008Lec3_Gauss Elimination

线性代数方程组的数值解法(1)
Gauss 消去法
武汉大学数学与统计学院 向华
a11x1  a12 x2    a1n xn  b1
Axb
a21 x1  a22 x2    a2 n xn  b2


an1 x1  an 2 x2    ann xn  bn
(Demos in Matlab: airfoil in 2D)
线性代数方程组的数值解法
直接法：Gauss 消去法，SuperLU
迭代法：定常迭代（Jacobi, GS, SOR, SSOR）
Krylov 子空间方法（CG, MINRES , GMRES, QMR, BiCGStab）
Direct
A = LU
Nonsymmetric
Symmetric
positive
definite
Iterative
y’ = Ay
Pivoting
LU
GMRES,
QMR, …
Cholesky
Conjugate
gradient
More General
More Robust
More Robust
Less Storage
The Landscape of Ax=b Solvers
刘徽 （约220-280）
Gauss（1777-1855）
Gaussian elimination, which first appeared in the text Nine Chapters on the Mathematical Art
written in 200 BC, was used by Gauss in his work which studied the orbit of the asteroid Pallas.
Using observations of Pallas taken between 1803 and 1809, Gauss obtained a system of six linear
equations in six unknowns. Gauss gave a systematic method for solving such equations which is
precisely Gaussian elimination on the coefficient matrix. (The MacTutor History of Mathematics,
http://www-history.mcs.st-andrews.ac.uk/history/index.html)
一个两千年前的例子
今有上禾三秉，中禾二秉，下禾
一秉，实三十九斗；上禾二秉，
中禾三秉，下禾一秉，实三十四
斗；上禾一秉，中禾二秉，下禾
三秉，实二十六斗。问上、中、
下禾实一秉各几何？答曰：上禾
一秉九斗四分斗之一。中禾一秉
四斗四分斗之一。下禾一秉二斗
四分斗之三。-------《九章算术》
3 x  2 y  z  39
2 x  3 y  z  34
x  2 y  3 z  26
一个两千年前的例子(2)
3 x  2 y  z  39

2 x  3 y  z  34
 x  2 y  3 z  26

3 x  2 y  z  39

5
1
y


3
3 z 8

8
4
y

3
3 z  13

3 x  2 y  z  39

5
1
y


3
3 z 8

33
12
z

5
5

z

 3 2 1 39 


 2 3 1 34 
 1 2 3 26


3

0
0

2 1 39 

5
1
8
3
3

8
4
13
3
3

3

0
0

2
1
5
3
1
3
12
5
0
39 

8
33 
5 
11
17
37
, y , x
4
4
4
Basic idea: Add multiples of each row to later rows to make A upper triangular
Axb
 a11

 a21
 

a
 n1
a12
a22

an 2
 a1n   x1   b1 
   
 a2 n   x2   b2 
 




  








 ann   xn   bn 
Solving linear equations is not trivial.
Forsythe (1952)
A(1)
 a11(1)
 (1)
 a21

 
 a (1)
 n1
a12(1)  a1(1n) 

(1)
(1)
a22  a2 n 
A


 
(1) 
an(12)  ann

li1 
ai(11)
(i  2,3,  , n)
a11(1)
aij( 2 )  aij(1)  li1a1(1j) (i, j  2,3,  , n)
A( 2 )
 a11(1)

 0

 
 0

a

 0
 0

 

 0
(1)
11
A( 3)
a12(1)  a1(1n) 

( 2)
( 2)
a22  a2 n 

 
( 2) 
an( 22)  ann

(1)
12
( 2)
22
a
a
0

0
(1)
13
( 2)
23
( 3)
33
a
a
a

an(33)
 a 

 a 
 a 
 
( 3) 
 ann

(1)
1n
( 2)
2n
( 3)
3n
li 2 
ai(22 )
(i  3,4,  , n)
( 2)
a22
aij(3)  aij( 2 )  li 2 a2( 2j) (i, j  3,4,  , n)
A( k )
A( k 1)
 a11(1)




 0

 0
 

 0
 a11(1)




 0

 0
 

 0
a1(,1k)1
a1(1k)


ak( k1,1k)1 ak( k1,1k)

0
akk( k )

0
ak( k)1,k



0
ank( k )


a1(,1k)1

ak( k1,1k)1
ak( k,k)1
ak( k)1,k 1

an( k,k)1

a1(,1k)1
a1(1k)
a1(,1k)1




ak( k1,1k)1

0

0


0
ak( k1,1k)
akk( k )
0

0






ak( k1,1k)1 
ak( k,k)1 
ak( k1,1k)1 

an( k,k11)

a1(1n) 

 
ak( k1,1n) 
(k )

akn

ak( k)1,n 
 

(k ) 
ann 
a1(1n) 

 

ak( k1,1n) 
akn( k ) 

ak( k1,1n) 
 

( k 1) 
ann 
aik( k )
lik  ( k ) (i  k  1, , n)
akk
aij( k 1)  aij( k 1)  lik akj( k )
(i, j  k  1, , n)
a







(1)
11
(1)
12
( 2)
22
(1)
13
( 2)
23
( 3)
33
a
a

a
a
a



(1)
x

a   1  b1 
    ( 2) 
a   x2   b2 


( 3) 


a   x3    b3 
      
 (n) 

(n)  
ann   xn   bn 
(1)
1n
( 2)
2n
( 3)
3n
Gauss 消去过程图示
After k=1
After k=2
After k=3
After k=n-1
用矩阵变换表达消去过程
A
(1)
A( 2 )
A( 3)
 a11(1)
 (1)
 a21

 
 a (1)
 n1
a12(1)  a1(1n) 

(1)
(1)
a22  a2 n 

 
(1) 
an(12)  ann

 a11(1)

 0

 
 0

a12(1)  a1(1n) 

( 2)
( 2)
a22  a2 n 

 
( 2) 
an( 22)  ann

 a11(1)

 0
 0

 

 0
a12(1)
a13(1)
( 2)
a22
0
( 2)
a23
( 3)
a33

0

an(33)
 a1(1n) 

( 2)
 a2 n 
 a3(3n) 
 
( 3) 
 ann

A
(1)
 a11(1)
 A  
 c1
def
r1T 
~ 
A1 
 1
 def
  I  l1e1T
L1   c1
  a(1) I n1 
 11

A
( 2)
a
 L1 A  
 0

(1)
(1)
11
 a11(1)
 def 
r
~  
c1r1T
 a (1)  A1  
11

 0
T
1
( 2)
a22
1


 def
  I  l2 e2T
L2  
1


  ac(22 ) I n  2 
22


A(3)  L2 A( 2)
 a (1)
 11


 0

( 2)
a22
0

r1T

T

r2
T
~ 
 ca2(r22)  A2 
22

c2
r1T 

T
r2 
~
A2 
用矩阵变换表达消去过程(2)
A( n )
 a11(1)


 Ln 1  L2 L1 A(1)  



1 1
1
2
1
n 1
a12(1)  a1(1n) 

( 2)
( 2 ) def
a22  a2 n 
U

  
(n) 
ann

def
A  L L  L U  LU
利用Gauss 变换矩阵的性质：(1) Li 1  I  li eiT   I  li eiT
1



(2) Li L j  I  li eiT I  l j eTj  I  li eiT  l j eTj (i  j )
单位下三角
形
1 1
1
2
1
n 1
L  L L L
 I  l e  l e  l e
T
1 1
T
2 2
T
n 1 n 1
 1
 (1)
 a21
 a11(1)
 a (1)
  31
(1)
 a11
 
 a (1)
 n(11)
a
 11
1
( 2)
a32
( 2)
a22

an( 22)
( 2)
a22
1

an(33)
( 3)
a33







 

 1 

用Gauss消去法求解 A x=b
--- LU分解 A = L U
--- 求解 L y = b
--- 求解 U x = y
(cost = 2/3 n3 flops)
(cost = n2 flops)
(cost = n2 flops)
算法实现： Gauss Elimination Algorithm

版本一
… 对第k列,消去对角线以下元素
… (通过每行加上第k行的倍数)
for i = k+1 to n
… 对第k行以下的每一行i
for j = k to n
… 第k行的倍数加到第 i 行
A(i,j) = A(i,j) - (A(i,k)/A(k,k)) * A(k,j)
for k = 1 to n-1

版本二: 在内循环中去掉常量 A(i,k)/A(k,k) 的计算
for k = 1 to n-1
for i = k+1 to n
m = A(i,k)/A(k,k)
for j = k to n
A(i,j) = A(i,j) - m * A(k,j)
Gauss Elimination Algorithm (2)

上一版本
for k = 1 to n-1
for i = k+1 to n
m = A(i,k)/A(k,k)
for j = k to n
A(i,j) = A(i,j) - m * A(k,j)

版本三: 第k列对角线以下为0,无需计算
for k = 1 to n-1
for i = k+1 to n
m = A(i,k)/A(k,k)
for j = k +1 to n
A(i,j) = A(i,j) - m * A(k,j)
Gauss Elimination Algorithm (3)

上一版本
for k = 1 to n-1
for i = k+1 to n
m = A(i,k)/A(k,k)
for j = k +1 to n
A(i,j) = A(i,j) - m * A(k,j)

版本四: 将乘子 m 存储在对角线以下备用
for k = 1 to n-1
for i = k+1 to n
A(i,k) = A(i,k)/A(k,k)
for j = k +1 to n
A(i,j) = A(i,j) - A(i,k) * A(k,j)
Gauss Elimination Algorithm (4)

上一版本
for k = 1 to n-1
for i = k+1 to n
A(i,k) = A(i,k)/A(k,k)
for j = k+1 to n
A(i,j) = A(i,j) - A(i,k) * A(k,j)

版本五: Split loop
for k = 1 to n-1
for i = k+1 to n
A(i,k) = A(i,k)/A(k,k)
for i = k+1 to n
for j = k+1 to n
A(i,j) = A(i,j) - A(i,k) * A(k,j)
Gauss Elimination Algorithm (5)

上一版本
for k = 1 to n-1
for i = k+1 to n
A(i,k) = A(i,k)/A(k,k)
for i = k+1 to n
for j = k+1 to n
A(i,j) = A(i,j) - A(i,k) * A(k,j)

版本六: 用矩阵运算
for k = 1 to n-1
A(k+1:n,k) = A(k+1:n,k) / A(k,k)
… BLAS 1 (scale a vector)
A(k+1:n,k+1:n) = A(k+1:n , k+1:n ) - A(k+1:n , k) * A(k , k+1:n)
… BLAS 2 (rank-1 update)
What we haven’t told you
 定理: 主元A(k)(k,k)不为0的充要条件是顺序主子矩阵非奇异
 定理: 分解的存在性和唯一性
 选主元策略(当主元A(k)(k,k)为0或很小时)
 向后误差分析
 并行技术
 块算法
 Sparse LU, Band LU
 最新进展(F.Gustavson & S.Toledo, Recursive Algorithm)
 还可用于矩阵求逆,求行列式,秩
BLAS 3 (Blocked) GEPP
for ib = 1 to n-1 step b … Process matrix b columns at a time
end = ib + b-1
… Point to end of block of b columns
apply BLAS2 version of GEPP to get A(ib:n , ib:end) = P’ * L’ * U’
… let LL denote the strict lower triangular part of A(ib:end , ib:end) + I
A(ib:end , end+1:n) = LL-1 * A(ib:end , end+1:n)
… update next b rows of U
A(end+1:n , end+1:n ) = A(end+1:n , end+1:n )
- A(end+1:n , ib:end) * A(ib:end , end+1:n)
… apply delayed updates with single matrix-multiply
… with inner dimension b
Distributed GE with a 2D Block Cyclic Layout
green = green - blue * pink
Matrix multiply of
Matlab 中的相应函数
inv
lu
\
Linpack 中对应的函数
sgea.f
sgefa.f
C, Fortran, Matlab 代码
function x = lsolve(A, b)
% x = lsolve(A, b) returns the solution to the equation Ax = b,
% where A is an n-by-b matrix and b is a column vector of
% length n (or a matrix with several such columns).
% Gaussian elimination with partial pivoting
[n, n] = size(A);
for k = 1 : n-1
% find index of largest element below diagonal in column k
max = k;
for i = k+1 : n
if abs(A(i, k)) > abs(A(max, k))
max = i;
end
end
% swap with row k
A([k max], :) = A([max k], :);
b([k max]) = b([max k]);
% zero out entries of A and b using pivot A(k, k)
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
b(k+1:n)=b(k+1:n)-A(k+1:n,k)*b(k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
%for i = k+1 : n
%alpha = A(i, k) / A(k, k);
%b(i) = b(i) - alpha * b(k);
%A(i, :) = A(i, :) - alpha * A(k, :);
%end
end
% back substitution
x = zeros(size(b));
for i = n : -1 : 1
j = i+1 : n;
x(i) = (b(i) - A(i, j) * x(j)) / A(i, i);
end
/* Computer Soft/c2-1.c Gauss Elimination */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define TRUE 1
/*
a[i][j] : matrix element, a(i,j)
n : order of matrix
eps : machine epsilon
det : determinant
*/
void main()
{
int i, j, _i, _r;
static n = 3;
static float a_init[10][11] = {{1, 2, 3, 6},
{2, 2, 3, 7},
{3, 3, 3, 9}};
static double a[10][11];
void gauss();
/*static int _aini = 1; */
printf( "\nComputer Soft/C2-1 Gauss Elimination \n\n" );
printf( "Augmented matrix\n" );
for( i = 1; i <= n; i++ ){
for( j = 1; j <= n+1; j++ ) {
a[i][j]=a_init[i-1][j-1]; printf( " %13.5e", a[i][j] );
}
printf( "\n" );
}
gauss( n, a );
printf( " Solution\n" );
printf( "-----------------------------------------\n" );
printf( "
i
x(i)\n" );
printf( "-----------------------------------------\n" );
for( i = 1; i <= n; i++ ) printf( " %5d %16.6e\n", i, a[i][n+1] );
printf( "-----------------------------------------\n\n" );
exit(0);
}
void gauss(n, a)
int n; double a[][11];
{
int i, j, jc, jr, k, kc, nv, pv;
double det, eps, ep1, eps2, r, temp, tm, va;
det = 1;
/* Initialization of determinant */
for( i = 1; i <= (n - 1); i++ ){
pv = i;
for( j = i + 1; j <= n; j++ ){
if( fabs( a[pv][i] ) < fabs( a[j][i] ) ) pv = j;
}
if( pv != i ){
for( jc = 1; jc <= (n + 1); jc++ ){
tm = a[i][jc]; a[i][jc] = a[pv][jc]; a[pv][jc] = tm;
}
det = -det;
}
if( a[i][i] == 0 ){
/* Singular matrix */
printf( "Matrix is singular.\n" ); exit(0);
}
for( jr = i + 1; jr <= n; jr++ ){ /* Elimination of below-diag */
if( a[jr][i] != 0 ){
r = a[jr][i]/a[i][i];
for( kc = i + 1; kc <= (n + 1); kc++ ){
temp = a[jr][kc];
a[jr][kc] = a[jr][kc] - r*a[i][kc];
if( fabs( a[jr][kc] ) < eps2*temp ) a[jr][kc] = 0.0;
/*
If the result of subtraction is smaller than
*
2 times machine epsilon times the original
*
value, it is set to zero. */
}
}
}
}
for( i = 1; i <= n; i++ ) {
det = det*a[i][i]; /* Determinant */
}
if( det == 0 ){
printf( "Matrix is singular.\n" ); exit(0);
}
else{ /* Backward substitution starts. */
a[n][n+1] = a[n][n+1]/a[n][n];
for( nv = n - 1; nv >= 1; nv-- ){
va = a[nv][n+1];
for( k = nv + 1; k <= n; k++ ) {va = va
- a[nv][k]*a[k][n+1];}
a[nv][n+1] = va/a[nv][nv];
}
printf( " Determinant = %g \n", det );
return;
}
}
C
C PAGE 220-223: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
C
C FILE: GAUSS.FOR
C
C GAUSSIAN ELIMINATION WITH SCALED PARTIAL PIVOTING (GAUSS,SOLVE,TSTGAUS)
C
DIMENSION A1(4,4),A2(4,4),A3(4,4),B1(4),B2(4),B3(4)
DIMENSION L(4),S(4),X(4)
DATA ((A1(I,J),I=1,4),J=1,4)/3.,1.,6.,0.,4.,5.,3.,0.,3.,-1.,7.,
A 0.,0.,0.,0.,0./
DATA (B1(I),I=1,4)/16.,-12.,102.,0./
DATA ((A2(I,J),I=1,4),J=1,4)/3.,2.,1.,0.,2.,-3.,4.,0.,-5.,1.,-1.,
A 0.,0.,0.,0.,0./
DATA (B2(I),I=1,4)/4.,8.,3.,0./
DATA ((A3(I,J),I=1,4),J=1,4)/1.,3.,5.,4.,-1.,2.,8.,2.,2.,1.,6.,
A 5.,1.,4.,3.,3./
DATA (B3(I),I=1,4)/5.,8.,10.,12./
C
CALL TSTGAUS(3,A1,4,L,S,B1,X)
CALL TSTGAUS(3,A2,4,L,S,B2,X)
CALL TSTGAUS(4,A3,4,L,S,B3,X)
END
SUBROUTINE TSTGAUS(N,A,IA,L,S,B,X)
DIMENSION A(IA,N),B(N),X(N),S(N),L(N)
PRINT 10,((A(I,J),J=1,N),I=1,N)
PRINT 10,(B(I),I=1,N)
CALL GAUSS(N,A,IA,L,S)
CALL SOLVE(N,A,IA,L,B,X)
PRINT 10,(X(I),I=1,N)
RETURN
10 FORMAT(5X,3(F10.5,2X))
END
SUBROUTINE GAUSS(N,A,IA,L,S)
DIMENSION A(IA,N),L(N),S(N)
DO 3 I = 1,N
L(I) = I
SMAX = 0.0
DO 2 J = 1,N
SMAX = AMAX1(SMAX,ABS(A(I,J)))
2 CONTINUE
S(I) = SMAX
3 CONTINUE
DO 7 K = 1,N-1
RMAX = 0.0
DO 4 I = K,N
R = ABS(A(L(I),K))/S(L(I))
IF(R .LE. RMAX) GO TO 4
J=I
RMAX = R
4 CONTINUE
LK = L(J)
L(J) = L(K)
L(K) = LK
DO 6 I = K+1,N
XMULT = A(L(I),K)/A(LK,K)
DO 5 J = K+1,N
A(L(I),J) = A(L(I),J) - XMULT*A(LK,J)
5
CONTINUE
A(L(I),K) = XMULT
6 CONTINUE
7 CONTINUE
RETURN
END
SUBROUTINE SOLVE(N,A,IA,L,B,X)
DIMENSION A(IA,N),L(N),B(N),X(N)
DO 3 K = 1,N-1
DO 2 I = K+1,N
B(L(I)) = B(L(I)) - A(L(I),K)*B(L(K))
2 CONTINUE
3 CONTINUE
X(N) = B(L(N))/A(L(N),N)
DO 5 I = N-1,1,-1