1 / 28

# 线性代数方程组的数值解法 (1) Gauss 消去法 - PowerPoint PPT Presentation

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

## PowerPoint Slideshow about '线性代数方程组的数值解法 (1) Gauss 消去法' - zorita-knowles

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript

Direct

A = LU

Iterative

y’ = Ay

More General

Non-

symmetric

Symmetric

positive

definite

More Robust

More Robust

Less Storage

Krylov 子空间方法（CG, MINRES , GMRES, QMR, BiCGStab）

The Landscape of Ax=b Solvers

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)

• Basic idea: Add multiples of each row to later rows to make A upper triangular

Gauss 消去过程图示

After k=1 After k=2 After k=3 After k=n-1

• --- LU分解 A = L U (cost = 2/3 n3 flops)
• --- 求解 L y = b (cost = n2 flops)
• --- 求解 U x = y (cost = n2 flops)

• 版本一

for k = 1 to n-1 … 对第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)

• 版本二: 在内循环中去掉常量 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)
• 还可用于矩阵求逆,求行列式,秩

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;

eps = 1.0; ep1 = 1.0 ; /* eps = Machine epsilon */

while( ep1 > 0 ){

eps = eps/2.0; ep1 = eps*0.98 + 1; ep1 = ep1 - 1;

}

eps = eps*2; eps2 = eps*2;

printf( " Machine epsilon=%g \n", eps );

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-diagonal. */

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 is calculated. */

}

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

SUM = B(L(I))

DO 4 J = I+1,N

SUM = SUM - A(L(I),J)*X(J)

4 CONTINUE

X(I) = SUM/A(L(I),I)

5 CONTINUE

RETURN

END